Constraint-Based Modeling: Diatom

Phaeodactylum tricornutum: Genome Scale Metabolic Model

Genome scale metabolic models (GSMMs) allow us to integrate synthetic genetic pathways for products of interest directly into the organism’s metabolic framework in silico. By doing this, we can predict how the introduction of these pathways can impact core cellular functions, resource allocation, and energy balance. Moreover, the GSMM can also simulate different environmental parameters like external CO₂ levels, nutrient availability, and light exposure. With this, we can predict how product synthesis will perform under Martian conditions. 

Thus, GSMMs enable us to identify potential metabolic bottlenecks - whether due to limited precursor availability, cofactor depletion, or energy constraints - that could limit the production of our products of interest. It also provides insights into how to bypass these bottlenecks through targeted genetic modifications or optimization of environmental conditions.

Finally, the predictions of the GSMM can be extended to optimizing the metabolic networks of the two species for the production of our products of interest. We can simulate various genetic and metabolic edits in silico, selecting the best candidates for experimental verification.

As a proof-of-concept, we simulated the metabolic network of our diatom, Phaeodactylum, in order to optimize it for the production of acetaminophen, a generic drug that can be immensely useful to humans who settle on Mars. In sync with our experimental work, we also simulated our bacterium, Pseudomonas, for the production of 4(S)-limonene, a useful platform chemical.  

Our model builds upon an established GSMM of Phaeodactylum tricornutum ^[5]. We have modified this to incorporate acetaminophen production. We began by focusing on the naturally-occuring anthranilate pathway, analyzing how the integration of two novel genes influences the entire metabolic system. 

The addition of the 4abh gene enables the production of 4-aminophenol from anthranilate, which, in the presence of the nhoA gene, leads to the synthesis of acetaminophen. 

4abh and nhoA do not catalyze any further reactions in the network. However, 4-aminophenol doesn't act alone - it is involved in various side reactions ^[6] driven by other enzymes naturally present in the diatom. Using Protein BLAST, we identified hits for enzymes in the diatom’s proteome that can catalyze these side reactions. By encoding these additional reactions into our enhanced model, we captured the full complexity of the metabolic network, providing a comprehensive system to predict and optimize acetaminophen production in P. tricornutum. 

Following the construction of our engineered genome scale metabolic model, we performed various analyses with different computational tools.

Flux Balance Analysis

Flux balance analysis (FBA) is a mathematical approach used to analyze the flow of metabolites through a metabolic network ^[12]. It is widely used in systems biology to predict the growth rate of an organism, the production of metabolites, and the response to environmental changes. FBA models the steady-state distribution of fluxes in a metabolic network, which are the rates at which metabolites are produced or consumed by reactions. 

This requires a stoichiometric matrix S, a matrix representing the stoichiometry of the metabolic network. The rows correspond to metabolites, and the columns correspond to reactions. The entry Sij  indicates the stoichiometric coefficient of the metabolite i in the reaction j. We then define the flux vector v to represent the fluxes / rates of each reaction in the network. 

We define the objective function Z as a linear combination of all reactions in the metabolic network. This defines a biological goal, such as maximizing biomass production or the ATP yield. This function is often chosen to reflect the growth or production objectives of the organism. Mathematically, it is expressed as Z = cTv , where c is a vector that specifies the contribution of each reaction to the objective function.

To perform flux balance analysis, we impose the constraint that the fluxes through the metabolic network are at steady state. Mathematically, this is represented by Sv = 0. Additionally, we can also impose linear constraints on the fluxes through particular reactions in the network based on prior knowledge of thermodynamics, kinetics, or other biological considerations. We represent such constraints with a < v < b, where a and b are vectors representing the bounds of the reactions.

Flux balance analysis solves the linear programming problem given by 

$$\max Z$$ $$ \text{such that } \text{S} \cdot \text{v} = 0 \text{ and } \text{a}\le \text{v} \le \text{b} $$

In our modeling, we used FBA to simulate the metabolic network of diatoms, specifically optimizing the flow of metabolites for acetaminophen production. By modeling steady-state conditions, FBA allowed us to predict outcomes like the maximum potential yield of acetaminophen and the optimal growth rate of diatoms under Martian environmental conditions. 

To verify if our model was accurate, we wanted to check the effect of nitrate in the growth media. It is well known that nitrate is a very important factor in the biomass of the diatom ^[13]. Our FBA simulation showed that the biomass flux shows a linear relationship with the nitrate content before it maxes out at around nitrate content of 4.55. 

We then optimized the model for biomass and acetaminophen production.

We found that the synthesis of glutamate from leucine and oxoglutarate is important for biomass but not acetaminophen production. We see the methyl oxopentanoate mitochondrial reaction is important for both acetaminophen and biomass production. Sodium transport also plays an important role in biomass optimization but has a lower flux when acetaminophen is to be optimized. Isoleucine transport was also found to be important for acetaminophen production, suggesting that it may play an important role in the biosynthesis of our product. The flux through cytochrome b6 was found to be higher when acetaminophen was optimized than when just the biomass was. Interestingly, the flux through the hydroquinone side reaction was higher when biomass was optimized, indicating its role in diatom growth. 

Flux Variance Analysis

We employed flux variance analysis (FVA)  ^[14] to explore the metabolic flexibility of the diatoms under optimal conditions. FVA allowed us to examine the range of possible fluxes through each pathway, revealing which metabolic reactions were essential and rigid and which displayed flexibility. This analysis was pivotal in uncovering potential bottlenecks in acetaminophen synthesis, as well as identifying alternative metabolic routes that could be leveraged to enhance production.

Parsimonious Flux Balance Analysis

Parsimonious Flux Balance Analysis (pFBA) is an extension of traditional FBA that introduces an additional optimization step to minimize the total flux through a metabolic network while still maximizing the same objective ^[15]. The goal of pFBA is to find the most efficient way to reach the objective, with the assumption that biological systems tend to minimize resource usage or metabolic costs. Mathematically, this is achieved by restating the linear programming problem as

$$\max Z$$ $$ \min \sum_{i=1}^{n} |v_{i}|$$ $$ \text{such that } \text{S } \cdot \text{v} = 0 \text{ and } \text{a} \le \text{v} \le \text{b} $$

We used pFBA to  refine the metabolic network of diatoms, allowing us to predict the most resource-efficient pathways for both biomass production and acetaminophen synthesis. By minimizing the overall flux through the network, pFBA generated biologically feasible predictions, ensuring that the diatoms used the least amount of resources and energy while maximizing metabolite yield. This efficiency is particularly critical in the resource-scarce Martian environment, where every molecule counts. pFBA enabled us to identify the most streamlined and energy-efficient pathways, guiding us toward an optimized metabolic design that balances survival and productivity in harsh conditions.

Metabolic Optimization and Modeling Analysis

In our modeling, Metabolic Optimization and Modeling Analysis (MOMA) was used to simulate how the diatoms respond to genetic modifications, particularly the integration of synthetic pathways for acetaminophen production. By assuming that the diatom’s metabolism would adapt swiftly and with minimal disruption to its existing flux distribution, MOMA allowed us to predict the immediate metabolic adjustments following these genetic changes ^[16]. 

This analysis was crucial for evaluating whether introducing the genes for acetaminophen synthesis would adversely affect the overall growth of the diatom. This approach ensured that our modifications optimized acetaminophen production while preserving the diatom's core metabolic functions

Regulatory On / Off Mechanism

We used Regulatory On / Off Mechanism (ROOM) ^[17] analysis for the fine-tuning of the diatom’s metabolic network in response to genetic modifications for acetaminophen production. By predicting the minimal set of reactions that need to be activated or deactivated, ROOM allowed us to strategically adjust the metabolic network to maximize acetaminophen yield while preserving essential functions like photosynthesis. This method was particularly valuable for optimizing the balance between these competing metabolic demands, ensuring that our diatoms could efficiently produce acetaminophen without compromising their core metabolic processes. 

We compared the fluxes of the reactions of interest across pFBA, MOMA, and ROOM, as they share a similar mathematical constraint structure where we look at minimizing the fluxes through the reactions. We observed only minor differences in the fluxes of pathways activated in MOMA and ROOM. The combined plot highlights these key flux variations while optimizing for both biomass and acetaminophen production. The small variability across the different methods suggests minor yet noticeable differences in the reactions that see some fluxes flowing through them.

Flux Scanning Based on Enforced Objective Flux

In our attempt to understand the intricacies of flow of metabolites, we used the Flux Scanning Based on Enforced Objective Flux (FSEOF) algorithm ^[18]. This analysis enabled us to pinpoint specific reactions within the diatom’s metabolic network that are crucial for maximizing acetaminophen yield. By systematically increasing the flux towards acetaminophen production, FSEOF meticulously scanned the entire metabolic network, revealing which reactions needed to be upregulated to achieve our production goals. This approach provided us with a targeted strategy to enhance the synthesis of acetaminophen, ensuring that every adjustment was focused on amplifying output while maintaining the overall efficiency of the diatom’s metabolism. We found the top 10 overexpression candidates, which are shown below: 

Escher

In order to visualize the fluxes through the diatom’s metabolic network, we turned to the tool Escher ^[19]. However, this is currently built to work only for E. coli and S. cerevisiae models. We successfully adapted Escher to visualize the metabolic pathways in diatoms, marking the first time Escher has been extended beyond its current intended use. 

We successfully integrated data from pFBA, MOMA, and ROOM to construct detailed metabolic maps. Our focus was on the pathways involved in synthesizing three essential amino acids: tryptophan, phenylalanine, and tyrosine. Given that chorismate, our target metabolite, is a precursor for these amino acids, our visualizations provide critical insights into how the introduction of new reactions impacts these pathways. As expected, we found that most of the flux was directed through the chorismate pathway when acetaminophen production was optimized. However, when biomass was optimized, the flux through the side reactions were found to increase.

Constraint-Based Modeling: Bacteria

Using COBRApy, we successfully introduced limonene synthase, enabling the biosynthesis of limonene from its precursor, geranyl diphosphate, into the GSMM for Pseudomonas ^[20]. By optimizing this metabolic pathway, we devised strategies to enhance limonene production.

Flux Balance Analysis

While optimizing limonene production using FBA, we observed that the flux for geranyl diphosphate originates primarily from a synthesis reaction catalyzed by dimethylallyltranstransferase (DMATT). The reaction is Dimethylallyl diphosphate + Isopentenyl diphosphate → Geranyl diphosphate + Diphosphate.

In the case of biomass optimization, the same pathway is activated, but the flux is significantly lower, indicating that although geranyl diphosphate plays a role in biomass formation, its impact is minimal, and the flux can remain low while still being essential.

When analyzing the key pathways and reactions activated during the optimization of limonene production and biomass formation, we observed that several common pathways are triggered. However, the ATP synthase reaction is crucial for biomass production, while it plays a much smaller role in limonene production, as its flux is significantly lower. Additionally, certain reactions, such as the NADH dehydrogenase reaction, and the diffusion of oxygen into the periplasm, show little to no activation, with minimal flux. Overall, although most of the same pathways are activated, their impact on the production of limonene or biomass varies.

Parsimonious Flux Balance Analysis

To minimize resource usage, we performed pFBA. Our simulations indicated that the flux through various key reactions can be reduced while still maintaining optimal biomass and limonene production. Interestingly, most of these reactions were related to glucose metabolism - glucose transport via diffusion or proton symport, and the glucose dehydrogenase reactions. 

Flux Scanning Based on Enforced Objective Flux

Overexpressing the genes associated with the above pathways leads to increased production of limonene. In this case, we focused on optimizing the production of both biomass and limonene. We identified the top 10 pathways and reactions that have the greatest impact on their fluxes.

Dynamic Modeling

Mechanisms for Silicate Solubilization

On Earth, soil microbes extract useful and nutritional compounds from the soil for cell plant uptake. This includes insoluble silicate sources that are inaccessible to the plant. They do this through a geochemical cycle called . Bioweathering is the primary source for transforming polymerized silica to monomeric forms ^[21].  

For silicate sources, the mechanisms for bioweathering and solubilization vary with composition. Metal-bound silicates require a combination of pH shifts and ligand interactions for effective release. On the other hand, minerals like quartz, silica, and phytoliths can often be dissolved through simpler processes such as proton exchange. Different bacterial species have different capabilities when it comes to solubilizing minerals. For example, earlier research found that Bacillus and Pseudomonas species were more effective at solubilizing magnesium trisilicate compared to silicate minerals like illite, quartz, and muscovite ^[22].

The three broadly accepted mechanisms for silicate solubilization are:

• lowering the pH of the environment through the production of organic and inorganic acids, 
• producing chelating compounds and ligands, and 
• performing nucleophilic attacks and exchange reactions. 

Among these, the most noticeable mechanism is acidolysis, where the pH around the microorganisms decreases due to their activity ^[23].

The rate of dissolution is influenced by the acidic conditions created by microbial activity, which shifts the dissolution rate away from equilibrium. Two organic acids that have been identified in media with silicate sources are gluconic acid and ketogluconic acid. While a specific marker acid for silicate deterioration remains unclear, our dynamic model aims to establish a positive correlation between acid concentration and silicate solubilization. 

Single Compartment Network Model

The network we constructed exclusively uses the organic acid buildup pathway for solubilization. To minimize noise and unavailable data, we focus on the proton contributions by gluconic acid and ketogluconic acid, the two biggest contributors to the acid pool. 

Our initial model was built considering that the acids are only in the cellular compartment of P. fluorescens and not secreted outside. This was done because we did not know the exact parameters for the acids’ exchange reaction and decided to assume that everything happens inside the cell. 

This model contained a growth media with glucose and other supplements. Inside the cell, the network contained:

• the organic acid production pathways 
• the silicate solubilization pathway

Dual Compartment Network Model

To stay biologically consistent, we also built a dual compartment set up containing intracellular and external spaces. Organic acids get produced inside the cell through acid synthesis pathways and then undergo exchange reactions to solubilize silica outside the cell. 

Broadly, the network can be classified into three pathways:

• the intracellular organic acid production pathways
• the discharge of these acids into external space
• the external solubilization 

Internal Concentrations

• Glucose: We assume a collective 75mM of internal glucose after making estimates about the glucose carrying capacity of various bacteria in different media ^[24].

• Gluconate : Another study conducted to find the absolute intracellular metabolic concentrations of E. coli ^[25] allowed us to make calculated estimates on concentrations of gluconate and gluconolactone in our P. fluorescens  model.

• Ketogluconate: We assumed that gluconate and ketogluconate are at the same concentration.

External Concentrations

While the external sample space of the bacteria is surrounded by soil, we assume that the bacterial population will utilize external metabolites in localized bubbles because of biofilm formation and other physical constraints.

• Glucose: The medium surrounding the bacteria was modeled to contain 2 g/L of glucose or 11 mM of external glucose ^[26].

• Silica: Taking references from the Martian soil simulant being used in our experimental work in the lab, and the amount of soil in the final co-culture setup, we estimate a total of 200 mM of SiO₂ in the soil.

• We assume that the Martian soil's external pH is 7.7 [27].

Reaction Rates

Using experimental data available in literature, we estimated the parameters of the reactions in our network model.

Results

Owing to data limitations in getting exact rate constants, our modeling work aimed at corroborating broad patterns established in literature as opposed to predicting specific rate values. Two of the claims that we set out to establish were:

• To show an inverse relationship of pH and the rate of solubilization: Since the solubilization process depends on organic acid buildup, it makes sense that a greater proton buildup (and the subsequent fall in pH) would result in faster solubilization.
• To show external pH converging to 4 or 5 ^[36].

While both our single and dual-compartment representations showed the inverse relationship between the pH and the rate of silicon solubilization, the relation is starkly seen in the single compartment set-up.

Our dual-compartment model was able to corroborate the pH conditions of the solubilization process incredibly well. Studies have shown that solubilization of silicates in a growth media by bacteria is generally accompanied by a fall in pH which then plateaus at a value in the ballpark of 4-5 after a week or so. In the case of siliceous Earth as the silica source, this value is 4.19 [36].

Our model shows the pH transiently dipping to values below 3.5 at peak solubilization. After this, the pH of the cell gradually increases to a value around 4. We also see the silicic acid content begin to saturate around the same time the pH starts stabilizing. 

Our dynamic network model for silicate solubilization effectively captures experimental data. The relation between pH and silicate solubilization was further corroborated by our wet lab work.

Community Modeling

We leveraged the MICOM (Metagenome-Scale Modeling to Infer Metabolic Interactions) framework to explore and analyze the intricate metabolic interactions within a co-culture system composed of our diatom and bacteria. MICOM is designed specifically to model and predict microbial interactions by integrating genome scale metabolic models, allowing researchers to study how different organisms co-metabolize substrates and influence each other’s growth dynamics. In general, a community model represents a group of individual models that interact within a shared environment. This environment is represented by a single compartment within which the organisms interact with each other and with entities outside the compartment.

Although the MICOM framework is used to study microbial communities with a large number of species, we suitably modified it to study a co-culture of a prokaryotic and eukaryotic species.

Community Building

As a GSMM for P. fluorescens was not available, we use the model for Pseudomonas putida as a proxy. To represent silicate solubilization and silicic acid uptake in our model, we incorporated a series of reactions that capture the complete process from how it is solubilized by the bacteria to how it can be taken up by the diatom as part of its biomass. We used a single compartment model for the bacterial solubilization, similar to our initial dynamic model.

After running the simulation, we found that the bacteria is indeed taking up the silica and solubilizing it. We also found that a fraction of the solubilized silica is taken up by the diatom for its biomass. The diatom is also seen to produce oxygen for the co-culture. Due to the photosynthetic nature of the diatom, the community takes up carbon dioxide and light for growth.

Community Media

We then used MICOM to determine the growth medium that optimizes the growth of both the bacteria and diatom in the co-culture. The ideal uptake fluxes that maximizes the growth are tabulated below:

The fluxes obtained are indicative of the relative concentrations of different media components required.

Growth Rates and Tradeoff

When we analyzed the relative growth rates of the two organisms in the optimized media, we found that the bacteria grows faster than the diatom. This was in good concurrence with our experimental results.

The tradeoff refers to the balance between optimizing individual species' growth and maintaining community-wide cooperation ^[37]. Each species in a microbial community has its own metabolic demands, and optimizing for maximum growth of one species can limit the resources available to others, leading to reduced overall community stability. MICOM incorporates this tradeoff by using a cooperative optimization approach, which ensures that species grow in a balanced way, sharing nutrients and resources more equitably, rather than one dominating the system. This tradeoff is particularly relevant in our co-culture experiments with Phaeodactylum tricornutum and Pseudomonas fluorescens, where silicate solubilization and resource allocation are critical for both organisms. Balancing the metabolic needs of both species allows for a more accurate simulation of their interactions and helps identify key metabolic exchanges driving growth and silicate uptake and utilization.

Our experimental hypothesis that increase in the bacteria’s growth directly correlates with the diatom’s growth was observed when we analyzed the tradeoff. The growth rates we observed in the optimized media were found to correspond to a tradeoff value of around 0.5, a perfect balance of cooperativity and self-utilization.

Effect of Silica on Growth

We analyzed the community formed between the diatom and bacteria by comparing it to a system without added silica, observing how the diatom's growth rates varied in the presence and absence of solubilized silica. 

We observed that there was a two-fold increase in growth when silica was added to the community and was solubilized by the bacteria. We observed the same results when we tested this experimentally.

Our community model demonstrates high accuracy, providing strong validation for its reliability and predictive power. This solidifies the model as a valuable tool for understanding and simulating diatom-bacteria interactions.