Flux Balance Analysis

Recent advances in computational biology led to emergence of databases with whole-genome metabolic networks, which for simplicity we will call models. They describe the flow of metabolites within a particular organism. Applying constraint-based methods (see this figure) to a particular model allows one to make quantitative predictions about cell phenotype while eliminating many complex parameters. Of particular interest to our project is flux balance analysis (FBA)^[1], which allows us to predict optimal steady-state biomass yield, which is tightly correlated with cell growth rate and is the most likely observed phenotype^[1,2]. One major advantage of FBA is that it does not require knowledge of enzymatic parameters. This figure^[3] elucidates the mathematics of FBA.

One of the primary questions for the project is the following. Will the synthetic P. denitrificans strain containing the device grow better than a wild type strain in the presence of ammonia? It is important to know the answer to this question because the To answer this question, we performed comparative analysis of the two Paracoccus strains using FBA on the corresponding whole-genome metabolic model. The analysis pipeline involved a slew of open-source computational tools, which we will describe below.

First, we reconstructed a metabolic model of Paracoccus denitrificans strain DSM 413 on a complete medium using ModelSEED^[4]. A complete medium is such that any nutrient, including ammonia, is available for uptake. Thus, the set of reactions included in the model is the largest of all possible sets. Although the largest, this set is incomplete because certain essential reactions that result in cell growth may be missing. In the next step, the model is gapfilled with the reactions necessary for measurable cell growth.

To quantify the difference between wild-type and synthetic strains, we must include several new reactions, metabolites and genes (e.g. oxygenation of ammonia by the AMO enzyme complex) into the model generated by ModelSEED. Such functionality is not available in ModelSEED. To implement this, we turned to COBRApy: Constraint-Based Reconstruction and Analysis^[5] package written in Python. However, COBRApy does not natively work with ModelSEED models. To overcome this compatibility issue, we used Mackinac package^[6] to convert the ModelSEED model into COBRApy-compatible format. Using COBRApy, we then added the new reactions into the model and performed FBA to compare the biomass yields, and hence the growth rates, of the two Paracoccus strains. The script is available here.

First, the FBA was run on the gapfilled unmodified (wild-type) model, which initially contained 1550 reactions and 1556 metabolites. The optimal biomass yield was found to be \( 224.3248 (\text{g dry weight}\cdot\text{h})^{-1} \). Next, we added the nitrification reactions, whose list is given below. \(\ce{Q}\) and \(\ce{QH_2}\) represent ubiquinone and ubiquinol, respectively. \( \text{UqO} \) is the ubiquinone oxidoreductase enzyme which catalyzes the last reaction, and is already present in the proteome of P. denitrificans. \[ \ce{NH_3 + QH_2 + O_2 ->[\text{AMO}] H_2O + Q + NH_2OH} \] \[ \ce{NH_3 + NAD + H_2O ->[\text{AMO}] 2H^+ + NADH + NH_2OH} \] \[ \ce{NH_2OH + O_2 ->[\text{HAO}] NO_2^- + H^+ + H_2O} \] \[ \ce{NH_2OH + 2Q + H_2O ->[\text{HAO}] NO_2^- + 2QH_2} \] \[ \ce{QH_2 ->[\text{UqO}] 2H^+ + Q} \] With the new model containing 1555 reactions and 1559 metabolites (with metabolites hydroxylamine, Q and QH2 added), the optimal biomass yield of the modified model was found to be \( \boxed{228.6980~(\text{g dry weight}\cdot\text{h})^{-1}} \).

Exchange Fluxes and Flux Variability Analysis

To gain a deeper understanding of how inclusion of new genes affects the metabolism of the cell, it is important to consider the difference in exchange fluxes between the wild type and synthetic strain. Exchange fluxes determine the intake and expulsion of certain extracellular metabolites, such as glucose, H⁺, or certain dipeptides. The set of exchange reactions also determines the composition of the growth medium. The set of fluxes that gives rise to the optimal biomass, however, is not unique. This includes exchange fluxes. The system of stoichiometric equations for FBA is usually underdetermined, which means that there are actually infinitely many solutions to the problem, all of which are optimal up to certain degree of significance. In other words, to investigate the significance of all exchange fluxes with respect to the optimum, we need to move away from analyzing specific numbers determined by FBA and instead consider ranges, or variations of fluxes that allow for the optimal solution.

To obtain the allowed range of fluxes, we employ a technique called Flux Variability Analysis (FVA)^[7]. The main parameter of FVA is the so called fraction of optimum, denoted by \(f\), which specifies the allowed deviation from the optimal biomass. In essence, during FVA we perturb each specified exchange flux in both directions under the constraint that the solution must within the fraction of optimum, thereby obtaining the allowed lower and upper bounds on the flux. If the fraction equals 1 (i.e. no deviation from the optimal biomass yield is allowed, in contrast to 0.95, which corresponds to 95% of the optimum), then any fluxes within the obtained range will reproduce the optimal solution. COBRApy package contains a module of functions which performs FVA on an existing model (see script for details).

The data obtained during FVA with \(f=0.999\) (i.e. a really tight variability constraint) confirms that the addition of nitrification circuit does, in fact, enhance the metabolism of the entire cell. However, the data are insufficient to determine how exactly the additional reactions and metabolites allow for such improvement.

Nitrite Exchange Knock-out and New Project Directions

To visualize metabolic pathways in our model, we used escher. The two figures below are based on the data from COBRApy simulations.

The first figure represents the nitrification-denitrification metabolic pathway of the strain with the device. As one can see in the first figure, all the nitrite produced is expelled from the cell through nitrite transport. This most likely happens due to cytotoxicity of nitrite^[8]. However, with nitrite transport genes knocked out, we can see in the second figure that all nitrite flows directly into the denitrification pathway. This offers a potential idea for a side project: how much will denitrification efficiency be improved due by knocking out nitrite exchange in Paracoccus denitrificans? This idea is worthwhile to explore in particular because knocking out nitrite exchange does not severely decrease the biomass yield (see figure captions).

Results, Summary and Discussion

Analysis of the data obtained from FBA revealed that the biomass yield of the P. denitrificans strain containing our device is greater than the wild-type yield by nearly 2%, which means that our bacterium will use the inserted nitrification pathway in order to enhance its metabolism. As such, our device potentially confers fitness advantage. We predict that our synthetic strain will out-compete the native wild-type strain if both strains are growing in the same environment, including sewage sludge. If our device is used in wastewater treatment facilities, we predict that there will be no need to artificially sustain the new culture.

As was stated earlier, the expected cell phenotype is the one with the largest biomass yield. This is only true under the assumption that the cell lives in an ideal or near-ideal environment with no over-population. Several studies have shown that under different growth environments or culture densities, cells can exhibit non-optimal yield metabolism^[9,10], which means that pathways different from those obtained through FBA will be employed by the cell. We do not know whether wastewater causes a similar shift in metabolism of Paracoccus denitrificans. However, seeing as it is able to thrive in such environment gives enough reason to believe that the growth rate will not lose correlation with the biomass. Furthermore, it should be noted that this model is oblivious to protein structure and enzyme kinetics, and as such is unable to capture the full scale of effects caused by inserting the nitrification construct into the model.

Description

In this section, we present a kinetic model of bioreactors which are able to perform both nitrification and denitrification. In particular, sequential batch reactors and bioreactors containing our synthetic strain of P. denitrificans are described by the model. Literature searching led to the discovery of a denitrification model called Activated Sludge Model with Indirect Coupling of Electrons (ASM-ICE), proposed by Pan et al. in 2013. ^[11] Mathematica notebook is available here.

In essence, the model is given by a system of coupled nonlinear time-dependent differential equations that produces for the following intermediate concentrations: \(S_{NO_3}\) - nitrate, \(S_{NO_2}\) - nitrite, \(S_{NO}\) - nitric oxide, \(S_{N_2O}\) - nitrous oxide, \(S_{N_2}\) - nitrogen gas. The rate of change of each of these variables is coupled to \(X\) - the heterotrophic biomass, \(S_S\) - carbon source concentration, and one of \(S_{Mox}\) or \(S_{Mred}\) - concentrations of oxidized (reduced) electron carriers. The relationships are summarized by a stoichiometric matrix below.

It is more instructive to interpret the table by looking at columns. Thus, for instance, we have \(\frac{dS_{NO}}{dt}=R3-R4\), where \(R3, R4\) are the expressions defined on the right side. The rates \(R1\) through \(R5\) depend on a large number of kinetic parameters which vary in different scenarios. We have chosen the set of parameters from Case 1 described in the paper as it is the closest to the experiment our wet lab team performed: "Case 1 (Kucera et al., 1983): The branching of the electron flow to individual terminal acceptors NO₃, NO₂ and N₂O was investigated in a pure denitrifying culture P. denitrificans (N.C.1.B. 8944). The culture was cultivated anaerobically to early stationary phase. A closed reactor with magnetic stirrer was used to carry out batch tests, during which N₂ was provided into the reactor to ensure oxygen free environment[...]"

Nitrification Amendment

Unfortunately, the above model does not describe the nitrification process. To that end, our main goal is to update the model in a proper way so that it incorporates both nitrification and denitrification, which is the full ammonia removal process. To do this, we need to include two new dependent variables (i.e. concentrations): \(S_{NH_3}\) - ammonia and \(S_{NH_2OH}\) - hydroxylamine. Furthermore, we would need to incorporate several new reaction rates that couple the concentrations. These rates would be of the same functional form as \(R1\) since each nitrification step consumes oxidized acceptors and returns their reduced counterparts:

Evidently, each such reaction would introduce 3 new kinetic parameters that would need to be estimated. With at least 2 new reactions (conversion of ammonia and hydroxylamine), the number of new parameters would be at least 6. Wolfram Mathematica has the capacity to obtain the best fit parameter estimates of a system of parametric nonlinear differential equations without knowing the closed form of the solutions. However, for more complicated models such as ours, the fitting becomes progressively slower as the dimension of the parameter space increases. This issue is resolved due to specifics of nitrification biochemistry: it was found that the concentration of hydroxylamine during ammonia oxidation is negligible^[12]. As such, we can assume that it gets instantly converted by the Hydroxylamine Oxidoreductase (HAO), thereby allowing us to write a single combined reaction for conversion of ammonia as long as we properly account for electron transfer.

Sometimes life happens and you get scooped. This happened to us less than a week away from the end of our research. It turns out that Hydroxylamine Oxidoreductase (HAO) converts hydroxylamine to nitric oxide (NO), as opposed to nitrite (NO₂^-) as we had originally thought^[13. Nitrite is a non-enzymatic product of a chemical reaction of NO with oxygen O₂. With this new knowledge, we had to change our course of action and define the combined reaction from the previous section to be \(\ce{NH_3 ->[\text{AMO+HAO}] NO}. \) Now we can finally write down the new equations to be added to the model:

Incorporating nitrification also changes electron carrier concentrations:

We are not yet able to simulate the model, for we do not know the values of the maximal reaction rate \(r_{NH_3,max}\) and the affinity constants \( K^N_{NH_3}, K_{Mox,1}. \)

Parameter Fitting and Sensitivity

Using the data we obtained during the bioreactor experiment (culture: N. europaea, ammonia spike at t=0), we ran the nonlinear model fit algorithm in Mathematica to obtain the best fit parameters. Their values are as follows:

Nonlinear fit performs best when provided with sufficiently good initial guess values. To prevent the algorithm from producing extremely unlikely results, we imposed the following constraints:

Unfortunately, using constrained optimization made it impossible to directly obtain the errors on these values. To amend that, we manually performed parameter sensitivity analysis to gain a better understanding of the errors. The results are summarized in the figure below.

Model Validation

To test the accuracy of the model, we have performed several simulations. In the first simulation, we set \(S_{NH_3}(0) = 34~\text{mmol/L}\) and the concentrations of all denitrification intermediates to zero. Furthermore, assume we are in excess of carbon source and have a large enough heterotrophic biomass. In this case, we expect that denitrification reactions do not produce enough oxidized electron carriers to support ammonia oxidation. The results of the simulation are in accord with these expectations.

In the second simulation, we reproduced Case 1 scenario from the Pan et al. paper and compared the nitrogen concentrations of the original and expanded models. We expect the nitrous oxide produced from ammonia oxidation to undergo conversion to nitrogen gas, the final product of ammonia removal. This is indeed the case, as is illustrated by the figure below.

Applications

Let us apply our model to solve a simple influent timing optimization problem. Suppose between windows of t₀=1.5 hours we have an influx of wastewater with a constant concentration of ammonia equal to 10 mmol NH₃/L and nitrite concentration of 2 mmol NO₃/L.

Zooming in on the first period, we see that nitrogen gas production stalls due to exhaustion of ammonia.

Mathematica can obtain the exact moment in time at which \(S_{NH_3}\) crosses a certain small threshold using an appropriate trigger WhenEvent[]. In our case, exhaustion occurs at \(t\approx 0.83~\text{h}.\) This tells us that to maximize the productivity of the bioreactor, we should receive an influent at or slightly lower than 0.83 hours. We can simulate this scenario just as easily.

And finally let us put nitrogen production rates side by side.

The average rate of nitrogen production is clearly higher in the optimized model, even with a cutoff imposed. Naturally, this example could be extended in many different ways to predict more realistic situations (e.g. randomly distributed influent concentrations), but in the absence of extensive concrete wastewater data measurements, we have decided to stop here.

Results and Discussion

Let us discuss some strengths and weaknesses of this model. On the pro side, the nitrification extension of the original model was proven to be faithful -- i.e., the model behaves as it is supposed to, retaining the full range of benefits of the original one. Secondly, by proxy of being implemented in Mathematica, our model offers a great deal of flexibility in choosing which scenarios can be simulated. Time-dependent and variable-dependent triggers in differential equations are extremely powerful tools in that respect. Only minor tweaks are required to simulate different initial conditions. Or, if one suspects that some kinetic parameters differ from those that we have found earlier, using our code one can fit the model to experimental data to find more accurate estimates.

On the other hand, there are several major drawbacks. First, this model is incomplete: non-enzymatic reaction \(\ce{NO + 1/2 O_2 -> NO_2^-}\) and the action of nitrite oxidoreductase (Nor) given by \(\ce{NO_2^- ->[\text{Nor}] NO_3^-}\) are not part of this model. Although it is possible to expand the model further and incorporate these reactions, the lack of sufficient data and a growing number of technical difficulties prevented us from doing so. Secondly, this model does not account for oxygen O₂ levels in the environment. Incorporating it into the model, given the complexity of its interactions with nearly all of the species present, is beyond the scope of a summer project, and most likely beyond the numeric capacity of Mathematica's algorithms.

The current age is, without a doubt, an age of numerics. Closed-form formulas amenable to rigorous mathematical analysis in mathematical biology are a rarity. Let us also not forget the rise of all the different -omics with their seemingly obscure meta-data that, upon heavy analysis, leads to a completely new understanding of the structure and function of living organisms. It is imperative that future generations of computational biologists embrace that many problems do not have a nice solution, and the only way to study them is to stand on the shoulders of numerical giants (whose names, for some reason, always start with letter M). We can only hope that one day, someone somewhere will stumble upon our code and repurpose it to their own project.

Generic template of the Metabolic Modeling script for use by other teams
Script used in the Metabolic Modeling section
Mathematica script with nitrogen removal kinetics