^{1} Civil and Environmental Engineering, University of Wisconsin-Madison,
^{2} Biomedical Engineering, Technical University of Denmark,
^{3} Chemical and Biological Engineering, University of Wisconsin-Madison

A major goal of systems biology is to develop mathematical models that can simulate an organism's biochemical characteristics, or phenotype, based on its biological components (e.g. genes, enzymes, metabolites). A common approach taken by systems biologists to model the phenotype of microorganisms rely on constraint-based reconstruction and analysis (COBRA) methods that simulate the flow of metabolites through an organisms metabolic network to predict system properties, such as rates of biomass growth or chemical production (O'Brien et al., 2015). COBRA methods have been under development since the mid-1990s and have become powerful tools for enabling predictive biology. This project focuses on the application of COBRA methods, specifically, flux balance analysis (FBA) (Orth et al., 2010), to model the biochemical process of *nitrification* mediated by two different groups of microorganisms.

Nitrification is an important step in the global nitrogen cycle, where ammonia (NH_{3}) is first oxidized to nitrite (NO_{2}^{-}) by ammonia oxidizing bacteria (AOB), and subsequently to nitrate (NO_{3}^{-}) by nitrite oxidizing bacteria (NOB). The process was discovered by Sergei Winogradsky in 1890 (Dworkin, 2012) and plays an important role in controlling the availability of nitrogen in soil and marine ecosystems, and removing ammonia from contaminated waters during water and wastewater treatment (see Nitrification primer).

In aerobic environments, the first step of nitrification (ammonia oxidation) it often done by bacteria affiliated with the genus *Nitrosomonas*, whereas the second step (nitrite oxidation) is often done by bacteria affiliated with the genus *Nitrobacter*. The overall chemical transformations of nitrification are shown below:

Step 1. $ \hspace{1cm} 2 \hspace{0.1 cm} NH_3 + 3 \hspace{0.1 cm} O_2 \rightarrow 2 \hspace{0.1 cm} NO_2^- + 2 \hspace{0.1 cm} H_2O + 4 \hspace{0.1 cm} H^+ \hspace{1cm}$(*Nitrosomonas*)

Step 2. $ \hspace{1cm} 2 \hspace{0.1 cm} NO_2^- + O_2 \rightarrow 2 \hspace{0.1 cm} NO_3^- \hspace{4.4cm}$(*Nitrobacter*)

The ultimate goal of this project was to predict the population dynamics and nutrient cycling (NH_{3}, NO_{2}^{-}, and NO_{3}^{-}) of a *Nitrosomonas*/*Nitrobacter* co-culture performing nitrification in a batch bioreactor, based on each organisms metabolic network stoichiometry and nutrient uptake kinetics (Figure 1). The ability to predict and simulate nitrification is essential for the design of biological wastewater treatment systems, and is a classic problem in environmental engineering. Moreover, being able to predict the nitrification process is critical to managing the global nitrogen cycle, which has been identified as a grand challenge for the 21st century by the National Academy of Engineering (NAE).

Model development was conducted in a stepwise manner to achieve our ultimate project goal. The growth of each organism was first modelled at steady-state to determine their maximum growth rate and nutrient uptake/export rates for NH_{3}, NO_{2}^{-}, and NO_{3}^{-}. A sensitivity and robustness analysis was also conducted to evaluate how each organisms growth was impacted by intracellular metabolite imbalances and changes in nutrient uptake rates. Subsequently, resulting steady-state models were solved over multiple timepoints to simulate the population dynamics and chemical transformations occuring in a bioreactor cycle.

Because several species are known to carry out the nitrification process, we were also interested in predicting the outcome of competiton experiments between different species of *Nitrobacter* involved in the second step of nitrification (NO_{2}^{-} oxidation). This was simulated using the dynamic microbial community model, where each species started at the same initial concentration and their abundance at the end of the cycle was determined.

The network stoichiometry used for *Nitrosomonas* and *Nitrobacter* were adapted from (Louca and Doebeli, 2015); kinetic parameters for nutrient uptake and affinity were adapted from Laanbroek et al., 1994 and Nowka et al., 2015.

Metabolic networks are comprised of compounds (nodes) and biochemical reactions (edges). An organisms metabolic network determines its physiological and biochemical properties, including what nutrients it consumes, how it generates energy (e.g. ATP) and reducing power (electrons in the form of NADH), and how it synthesizes biomass (i.e. proteins, lipids, nucleic acids, carbohydrates). The proportion of reactants to products for a given biochemical reaction is termed stoichiometry. Metabolic networks can be represented mathematically as a matrix of compounds (rows) and reactions (columns), where the matrix entries represent stoichiometric coefficients. Figure 2 provides a conceptual example of how a metabolic network is converted to a mathematical representation.

**Nitrifier metabolism and stoichiometry**

Nitrifying bacteria, such as *Nitrosomonas* and *Nitrobacter*, are "chemolithoautotrophic" microorganisms, meaning they grow using CO_{2} as a carbon source (similar to plants) and use inorganic substrates (of mineral origin) to generate energy and reducing power (electrons) for biosynthesis and cellular maintenance. For nitrifying bacteria, inorganic nitrogen compounds are used as primary substrates: *Nitrosomonas* use NH_{3}, whereas *Nitrobacter* use NO_{2}^{-}. Like humans, nitrifying bacteria also require O_{2} as an electron acceptor for energy production via aerobic respiration.

For this project, simplifed metabolic networks for *Nitrosomonas* and *Nitrobacter* were used to model cell growth (Louca and Doebeli, 2015). Networks were first converted to stoichiometric tables in excel and saved as csv files. Resulting tables were then read into Julia as matrices and stored as NamedArrays. The total number of compounds and reactions for each metabolic network are provided on Table 1.

Model | Reactions | Compounds |
---|---|---|

Nitrosomonas | 30 | 28 |

Nitrobacter | 26 | 24 |

Figure 3 shows the metabolic network for *Nitrosomonas* used in this project. The details of the metabolic network for *Nitrobacter* can be found in the supplementary information.

**Monod kinetics**

The kinetics of substrate uptake for *Nitrosomonas* and *Nitrobacter* were modelled using the Monod equation:

$ \hspace{1.7cm} r^k_{uptake,i} = V_{m,i}^k(\frac{C^k_i}{K_{m,i}^k + C^k_i}) \hspace{0.5cm} (1)$

Where,

$r^k_{uptake,i}$ = Uptake rate of substrate *i* by organism *k*

$V_{m,i}^k$ = maximum uptake rate of substrate *i* for organism *k*

$K_{m,i}^k$ = half-saturation coefficient of substrate *i* for organism *k*

$C^k_i$ = Extracellular concentration of substrate *i* available to organism *k*

The Monod equation is commonly used to model kinetic limitations on substrate uptake by microorganisms. Under non-limiting conditions, substrate uptake occurs at the maximum rate, *V _{m}*. Under substrate-limiting conditons, the substrate uptake rate is based on the organisms affinity for the substrate (

For our modelling purposes, kinetic constraints were only imposed on the primary substrates for each organism; NH_{3} for *Nitrosomonas* and NO_{2}^{-} for *Nitrobacter*. All other nutrients required for growth (e.g. phosphorous, CO_{2}, O_{2}, etc.) were assumed to be non-limiting (see Section 2.B. for further details).

In [1]:

```
#Nitrosomonas stoichiometry and kinetics
using NamedArrays
#read stoichiometry table as matrix
rawNS = readdlm("Lawson_Nitrosomonas.csv",',');
(ms,ns) = size(rawNS)
#range of reactions and metabolites in matrix
ns_reactions = 2:ns
ns_metabolites = 2:ms
#names of reactions and metabolites
reactionsNS = rawNS[1,ns_reactions][:]
metabolitesNS = rawNS[ns_metabolites,1][:]
#collect exchange reactions
EXrxnsNS = []
for i in reactionsNS
if ismatch(r"EX_",i) == true
push!(EXrxnsNS,i)
end
end
#collect intracellular reactions
INTrxnNS = []
for i in reactionsNS
if ismatch(r"EX_|T_",i) == false
push!(INTrxnNS,i)
end
end
#stochiometry matrix
stoicmatrixNS = NamedArray( Array{Float64}(rawNS[ns_metabolites,ns_reactions]), (metabolitesNS,reactionsNS),("metabolitesNS","reactionsNS") )
println(stoicmatrixNS)
#Nitrosomonas europaea kinetics
KdNS = 0.04 #Biomass decay coefficient (days-1)
VmNS = 2.14e-13 #Nitrosomonas max growth rate (mol NH3/cell-d)
KmNS = 1.050 #Half-saturation value for NH3 (mM)
wt_NS = 3e-13 #cell weight (g DW/cell)
Vm_NS = (VmNS*1000)/wt_NS #(mmol/gDW-d)
MW = 1198/1000 #gDW/mmol, calculated from Smatrix ;
```

In [2]:

```
#Nitrobacter stoichiometry and kinetics
using NamedArrays
#read stoichiometry table as matrix
rawNB = readdlm("Lawson_Nitrobacter.csv",',')
(mb,nb) = size(rawNB)
#range of reactions and metabolites in matrix
nb_reactions = 2:nb
nb_metabolites = 2:mb
#names of reactions and metabolites
reactionsNB = rawNB[1,nb_reactions][:]
metabolitesNB = rawNB[nb_metabolites,1][:]
#collect exchange reactions
EXrxnsNB = []
for i in reactionsNB
if ismatch(r"EX_",i) == true
push!(EXrxnsNB,i)
end
end
#collect intracellular reactions
INTrxnNB = []
for i in reactionsNB
if ismatch(r"EX_|T_",i) == false
push!(INTrxnNB,i)
end
end
#stochiometry matrix
stoicmatrixNB = NamedArray( Array{Float64}(rawNB[nb_metabolites,nb_reactions]), (metabolitesNB,reactionsNB),("metabolitesNB","reactionsNB"))
println(stoicmatrixNB)
#Nitrobacter winogradskyi kinetics
KdNB = 0.04 #Biomass decay coefficient (days-1)
VmNB = 8.88e-14 #Nitrobacter max growth rate (mol NO2/cell-d)
KmNB = 0.309 #Half-saturation value for NO2 (mM)
wt_NB = 4e-13 #cell weight (g DW/cell)
Vm_NB = (VmNB*1000)/wt_NB #(mmol/gDW-d)
MW = 1198/1000 #gDW/mmol, calculated from Smatrix;
```

Flux balance analysis (FBA) is a mathematical optimization approach for analyzing the flow of metabolites through a metabolic network (Orth et al., 2010). Here, flux is defined as the molar flow rate of molecules through a given reaction per unit biomass [mmol /g of dry cell weight(DW) /hour]. Both *Nitrosomonas'* and *Nitrobacter's* metabolic network were initially analyzed using FBA to determine their maximum growth rate and corresponding uptake/export rates for NH_{3}, NO_{2}^{-}, and NO_{3}^{-} at steady-state. The decision variables for the FBA problem were the metabolic fluxes through each reaction. The objective function was to maximize flux through the biomass synthesis reaction (i.e. maximize cell growth), subject to the following constraints:

- Mass conservation: the stoichiometry matrix imposes flux balance constraints on the network, such that the total production rate of a given metabolite must equal its total consumption rate (i.e. steady-state constraint).
- Upper and lower flux bounds on each reaction.
- Substrate uptake rate constraints set by experimentally determined kinetic parameters.

Figure 4 provides a conceptual image of the biomass objective function.

The problem is typically written as:

$$ \begin{aligned} \underset{v \in \mathbb{R^n} }{\text{maximize}}\qquad& \mu^k = c^Tv_j^k = v^k_{biomass} \hspace{0.5cm} (2)\\ \text{subject to:}\qquad& \sum S^k_{i,j}v^k_j = b_i && \forall \hspace{0.1cm} i \in I \hspace{0.5cm} (3)\\ & LB^k_j \le v^k_j \le UB^k_j && \forall \hspace{0.1cm} j \in J \hspace{0.5cm} (4)\\ & v^k_{uptake,i} \le V_{m,i}^k(\frac{C^k_i}{K_{m,i}^k + C^k_i}) && \hspace{0.1cm} i = NH_3, NO_2^- \hspace{0.5cm} (5)\\ \end{aligned} $$Where,

*Sets*

$K$ = Set of organisms

$J^k$ = Set of reactions in metabolic network of organism *k*

$I^k$ = Set of metabolites in metabolic network of organism *k*

*Parameters*

$S^k_{i,j}$ = Stoichiometric coefficients of metabolite $i \in I^k$ in the reaction $j \in J^k$

$LB^k_j$ and $UB^k_j$ = Lower and upper flux bound on reaction $j \in J^k$ for community member *k*

$V_{m,i}^k$ = maximum uptake rate of substrate *i* for organism *k*

$K_{m,i}^k$ = half-saturation coefficient of substrate *i* for organism *k*

$C^k_i$ = Extracellular concentration of substrate *i* available to organism *k*

$c^T$ = cost array containing all zeros except for a 1 in the position corresponding to the biomass flux reaction

$b_i$ = the rate of accumulation/depletion of metabolite *i*, assumed here to be zero (i.e. steady-state)

*Variables*

$v^k_j$ = Metabolic flux through reaction *j* in community member *k*

$v^k_{uptake,i}$ = substrate uptake flux by community member *k*

Because the righthand side of Eq. 5 reduces to a constant, the problem stated above can be solved using linear programming. A function for performing FBA is coded in Julia + JuMP below.

**Dual problem and shadow prices**

The sensitivity of FBA to flux imbalances can be captured by the dual problem and its corresponding variables (shadow prices):

Where,

$\lambda_i$ = the shadow prices for the steady-state constraints on metabolite *i*

$q_{1,j}$ = the shadow prices for the lower bound flux constraints on reaction *j*

$q_{2,j}$ = the shadow prices for the upper bound flux constraints on reaction *j*

$\gamma_j$ = the shadow prices for the substrate uptake constaints on reaction *j*

$r_i^k$ = the constant obtained from evaluating the right hand side of Eq. 5

The shadow prices for the metabolites ($\lambda_i$) are most important to us and quantify the "value" of each metabolite to the objective function; that is, how sensitive the objective function is to flux imbalances. A negative $\lambda_i$ means that adding metabolite *i* would increase the biomass growth rate. A positive $\lambda_i$ means that removing metabolite *i* would increase the biomass growth rate. A $\lambda_i$ of zero means that adding or removing metabolite *i* would have no affect on the biomass growth rate.

The dual values are automatically output from the FBA function coded below and are used in subsequent analysis detailed in this report.

**Model assumptions**

Our FBA model is based on two key assumptions:

That is, there is no intracellular accumulation or depletion of metabolites in the network. This assumption is justified becuase the rate of turnover of most metabolic intermediates is high relative to intracellular concentrations. Transient effects of metabolism due to perturbations are thus resolved quickly relative to the rate of cell growth.*Metabolic networks operate under steady-state conditions.*While the use of appropriate objective functions (e.g. minimize ATP usage, minimize ribosome dilution, etc.) is currently an open topic in metabolic engineering, the maximization of growth assumption works well for microbial systems in suspended cell culture. Due to the exponential nature of cell growth, the fastest growing cells will always be dominate in the total cell population. In other words, rapidly growing phenotypes will always outcompete slower variants of the same species.*The objective of a cell is to maximize its growth rate.*

Other assumptions relating to the physical environment were also made. These include ignoring the effects of temperature, pH, and diffusion limitations on the chemical reactions considered in the metabolic network.

In [3]:

```
using JuMP, Gurobi
#solve steady-state FBA model
function fbaNS(u,C) #set uptake rate (value) for reaction (key)
λNS = Dict() #dictionary of metabolite shadow prices
z = Model(solver=GurobiSolver(OutputFlag=0))
@variable(z, flux[reactionsNS]) #flux variables for each reaction
@variable(z, t[key in keys(u)]) #t equals absolute value of flux["r"]
#steady-state metabolite constraints
@constraint(z, stoichiometry[i in metabolitesNS], sum{flux[j]*stoicmatrixNS[i,j], j in reactionsNS} == 0)
@constraint(z, direction[i in INTrxnNS], flux[i] >= 0) #Directionality constraints on reaction fluxes
@constraint(z, flux["EX_NO"] == 0)
@constraint(z, c1[key in keys(u)], t[key] <= u[key]) #uptake rate limited by nutrient availability
@constraint(z, c2[key in keys(u)], t[key] >= flux[key] ) #key sets exchange flux reaction
@constraint(z, c3[key in keys(u)], t[key] >= -flux[key] ) #key sets exchange flux reaction
#If nutrient concentration is 0, uptake flux must be zero
if C[:NH3] <= 0
@constraint(z, flux["EX_NH3"] >= 0 )
end
if C[:NO2] <= 0
@constraint(z, flux["EX_NO2"] >= 0 )
end
if C[:NO3] <= 0
@constraint(z, flux["EX_NO3"] >= 0 )
end
@objective(z, Max, flux["R_biomassNS"]) #maximize biomass growth
solve(z)
#output model solutions
global μ = (getobjectivevalue(z))
global rNH3 = getvalue(flux["EX_NH3"])
global rNO2 = getvalue(flux["EX_NO2"])
global rNO3 = getvalue(flux["EX_NO3"])
global fluxes = getvalue(flux)
for i in metabolitesNS
λNS[i] = getdual(stoichiometry[i])
end
global λNS
return μ, rNH3, rNO2, rNO3, fluxes, λNS
end
#Determine max growth rate
C1=Dict(:NH3=>100 ,:NO2=>0.0 ,:NO3=>0.0 )
rate = Vm_NS
u = Dict("EX_NH3" => rate)
outputNS = Dict(:a => [])
push!(outputNS[:a],fbaNS(u,C1)) #run model
#print model results
println("μ: ",μ, ", rNH3: ", rNH3,", rNO2: ",rNO2, ", rNO3: ",rNO3)
for i in λNS
#println(i) #shadow price for metabolite i
end
println("reaction fluxes: ")
for j in reactionsNS
println("$j: ",fluxes[j]) #fluxes
end
```

In [4]:

```
using JuMP, Gurobi
#solve steady-state FBA model
function fbaNB(u, C) #set uptake rate (u) for reaction (r)
λNB = Dict() #dictionary of metabolite shadow prices
z = Model(solver=GurobiSolver(OutputFlag=0))
@variable(z, flux[reactionsNB]) #flux variables (mmol/mgDW/day)
@variable(z, t[key in keys(u)]) #t equals absolute value of flux["r"]
#steady-state metabolite constraints
@constraint(z, stoichiometry[i in metabolitesNB], sum{flux[j]*stoicmatrixNB[i,j], j in reactionsNB} == 0)
@constraint(z, direction[i in INTrxnNB], flux[i] >= 0) #Directionality constraints on reaction fluxes
@constraint(z, flux["EX_NO"] == 0)
@constraint(z, c1[key in keys(u)], t[key] <= u[key]) #uptake rate limited by nutrient availability
@constraint(z, c2[key in keys(u)], t[key] >= flux[key] ) #n equals exchange flux reaction
@constraint(z, c3[key in keys(u)], t[key] >= -flux[key] ) #n equals exchange flux reaction
#If nutrient concentration is 0, uptake flux must be zero
if C[:NH3] <= 0
@constraint(z, flux["EX_NH3"] >= 0 )
end
if C[:NO2] <= 0
@constraint(z, flux["EX_NO2"] >= 0 )
end
if C[:NO3] <= 0
@constraint(z, flux["EX_NO3"] >= 0 )
end
@objective(z, Max, flux["R_biomassNB"]) #maximize biomass growth
solve(z)
#output model solutions
global μ = (getobjectivevalue(z))
global rNH3 = getvalue(flux["EX_NH3"])
global rNO2 = getvalue(flux["EX_NO2"])
global rNO3 = getvalue(flux["EX_NO3"])
global fluxes = getvalue(flux)
for i in metabolitesNB
λNB[i] = getdual(stoichiometry[i])
end
global λNB
return μ, rNH3, rNO2, rNO3, fluxes , λNB
end
#Determine max growth rate
C2=Dict(:NH3=>0.0,:NO2=>1.0,:NO3=>0.0)
u = Dict("EX_NO2" => Vm_NB)
test = Dict(:b => [])
push!(test[:b],fbaNB(u, C2)) #run model
#print results
println("μ: ",μ, ", rNH3: ", rNH3,", rNO2: ",rNO2, ", rNO3: ",rNO3)
for i in λNB
#println(i) #metabolite shadow prices
end
println("reaction fluxes:")
for j in reactionsNB
println("$j: ",fluxes[j]) #fluxes
end
```

Phenotypic phase plane (PhPP) analysis was conducted to determine the sensitivity of the objective function (biomass growth rate) to simultaneous changes in oxygen and ammonia (two key metabolites for growth). In PhPP analysis, two uptake fluxes are varied simultaneously on an x,y-plane and used to describe changes in the objective function (Edwards et al., 2002). The optimal biomass flux solutions are calculated for all points examined on this plane and the resulting regions or "phases" represent different metabolic network utilizations (i.e. phenotypes).

The PhPP can be defined using the shadow prices from the solutions to the dual FBA problems. Once all points have been evaluted in the two-dimensional plane, the shadow prices for all metabolites are calculated within this region and lines are drawn to demarcate regions of constant shadow prices. Figure 5 shows a conceptual PhPP.

The isoclines or contours of the phase plane represent the combination of uptake fluxes that result in the same objective value. The slope of the isocline can be calculated from the ratio of shadow prices for the metabolites corresponding to each uptake flux as follows:

$ \hspace{1.7cm} \alpha = -\frac{\lambda_A}{\lambda_B} = -\frac{\frac{dZ}{db_A}}{\frac{dZ}{db_B}} = \frac{db_B}{db_A} \hspace{0.5cm} (8)$

Where,

$\alpha$ = the relative change in the objective function for the two examined exchange fluxes,

$\lambda$ = the shadow price of metabolite A or B, and

$\frac{dZ}{db}$ = the change or sensitivity of the objective function (Z) to changes in the availability of metabolites A or B.

**Solution**

PhPP analysis was only conducted for *Nitrosomonas*. Here, the maximal growth rate of *Nitrosomonas* was examined as a function of simultaneously varying the oxygen uptake and ammonia uptake fluxes. A line representing the optimal relationship between the two metabolic fluxes was also determined. This "line of optimality" (LO) was calculated by fixing the ammonia uptake flux along the x-axis and determining the optimal biomass flux using FBA. The code for PhPP is provided below in Julia + JuMP.

In [5]:

```
#Phenotype Phase Plane Analysis (PhPP)
#inputs
nPts = 50 #number of points to plot in each dimension
range1 = 700 #range of reaction 1 to plot
range2 = 700 #range of reaction 2 to plot
rxnID = ["EX_NH3", "EX_O2"]
metID = ["NH3_NH4_e", "O2_c"]
#create empty vector for PhPP results
C3=Dict(:NH3=>100 ,:NO2=>100 ,:NO3=>100 )
b1 = collect(linspace(0, range1, nPts))
b2 = collect(linspace(0, range2, nPts))
b = Dict(rxnID[1] => b1, rxnID[2] => b2 )
PhPP = []
growthRates = []
shadowPrice1 = []
shadowPrice2 = []
#create empty vectors for LO results
LO = []
growthRatesLO = []
rxn2Flux = []
#calculate PhPP points
for i in 1:nPts
for j in 1:nPts
b_sub = Dict()
b_sub[rxnID[1]] = b[rxnID[1]][i] #flux limit on NH4
b_sub[rxnID[2]] = b[rxnID[2]][j] #flux limit on O2
push!(PhPP,fbaNS(b_sub, C3) ) #push values to list
end
end
#collect growthRates
for i in 1:length(PhPP)
push!(growthRates,PhPP[i][1] ) # optimal growth rate
end
growthRates = reshape(growthRates,nPts,nPts);
#collect shadowPrices
for i in 1:length(PhPP)
push!(shadowPrice1,PhPP[i][6][metID[1]] ) #shadow prices for NH4
push!(shadowPrice2,PhPP[i][6][metID[2]] ) #shadow prices for O2
end
shadowPrice1 = reshape(shadowPrice1,nPts,nPts);
shadowPrice2 = reshape(shadowPrice2,nPts,nPts);
#calculate line of optimality
for i in 1:nPts
b_sub = Dict()
b_sub[rxnID[1]] = b[rxnID[1]][i] #flux limit on NH4
push!(LO,fbaNS(b_sub, C3) )
end
#collect LO results
for i in 1:nPts
push!(growthRatesLO, LO[i][1]) # LO growth rate
push!(rxn2Flux, -1*LO[i][5][rxnID[2]])
end
#print results
println("NH3 limited region")
println("λNH3: ",shadowPrice1[35,7],", pt 37,7")
println("λNH3: ",shadowPrice1[42,24],", pt 42,24")
println("λNH3: ",shadowPrice1[14,3],", pt 14,3")
println("λNH3: ",shadowPrice1[21,14],", pt 21,14")
println("λNH3: ",shadowPrice1[21,6],", pt 21,6")
println("λO2: ",shadowPrice2[35,7],", pt 35,7")
println("λO2: ",shadowPrice2[42,24],", pt 42,24")
println("λO2: ",shadowPrice2[14,3],", pt 14,3")
println("λO2: ",shadowPrice2[21,14],", pt 21,14")
println("λO2: ",shadowPrice2[21,6],", pt 21,6")
println("O2 limited region")
println("λNH3: ",shadowPrice1[25,46],", pt 25,46")
println("λNH3: ",shadowPrice1[4,38],", pt 4,38")
println("λNH3: ",shadowPrice1[2,14],", pt 2,14")
println("λNH3: ",shadowPrice1[13,21],", pt 13,21")
println("λNH3: ",shadowPrice1[42,42],", pt 42,42")
println("λO2: ",shadowPrice2[25,46],", pt 25,46")
println("λO2: ",shadowPrice2[4,38],", pt 4,38")
println("λO2: ",shadowPrice2[2,14],", pt 2,14")
println("λO2: ",shadowPrice2[13,21],", pt 13,21")
println("λO2: ",shadowPrice2[42,42],", pt 42,42")
```

Following our steady-state analysis of each organism, both FBA models for *Nitrosomonas* and *Nitrobacter* were combined to simulate the population dynamics and nutrient cycling occuring in a fed-batch bioreactor performing nitrification over one reactor cycle. At time zero, 1 mM of NH_{3} was added to the bioreactor. Subsequently, both FBA models were solved over multiple timepoints, where solution outputs from each point (biomass growth rate and substrate uptake/export rates) were used to update the biomass and extracellular nutrient concentrations (NH_{3}, NO_{2}^{-}, and NO_{3}^{-}) used for the following time step. This formulation has previously been applied to model multi-species microbial communities of *Rhodoferax* and *Geobacter* in anoxic subsurface environments (Zhuang et al., 2011).

Ordinary differential equations for mass conservation around biomass and nutrients in the bioreactor were solved to determine the biomass and nutrient concentrations at each timepoint:

$\hspace{1cm} \frac{dX^k}{dt} = (\mu^k - K_d^k)X^k \hspace{0.5cm} \forall \hspace{0.1cm} k \in K \hspace{0.5cm} (9)$

$ \hspace{1cm} \frac{dC_i}{dt} = \sum r^k_iX^k \hspace{0.5cm} i = NH_3, NO_2^-, NO_3^- \hspace{0.5cm} (10)$

Where,

$C^k_i$ = Extracellular concentration of shared metabolite *i* available to organism *k*

$X^k$ = Biomass concentration of organism *k*

$\mu^k = MW*v^k_{biomass}$ = biomass flux for organism *k* times biomass molecular weight

$r^k_i$ = uptake/export rate of shared metabolite *i* for organism *k*

$K_d^k$ = Biomass decay rate for organism *k*

**Competition experiment**

In the original FBA model, kinetic parameters for *Nitrobacter* were taken from the species *Nitrobacter winogradskyi*. In order to model competition for NO_{2}^{-} between multiple organisms, a new species, *Nitrobacter vulgaris*, was added to the model and simulations were run to determine which organism would proliferate over the other in the bioreactor. For simplicity, the network stoichiometry between *N. vulgaris* and *N. winogradskyi* was assumed to be the same, but the maximum uptake rate and affinity parameters for NO_{2}^{-} were changed to match those determined in previous studies (Nowka et al., 2015). The Julia + JuMP code for implementing the dyanmic microbial community modelling is provided below.

In [6]:

```
#Bioreactor model
using JuMP
#initial conditions
species = [:NS,:NB, :NV] #list of species
nutrients = [:NH3 :NO2 :NO3] #list of extracellular metabolites
EXnut = ["EX_NH3","EX_NO2","EX_NO3"] #exchanges fluxes for extracellular metabolites
X0_NS = float(1e10)*wt_NS #inital Nitrosomonas biomass conc (gDW/L)
X0_NB = float(1e10)*wt_NB #inital Nitrobacter biomass conc (gDW/L)
X0_NV = float(1e10)*wt_NB #inital N. vulgaris biomass conc (gDW/L)
C=Dict(:NH3=>[1.0],:NO2=>[0.0],:NO3=>[0.0]) #Dict of initial extracellular metabolite conc (mM)
X=Dict(:NS =>[X0_NS], :NB => [X0_NB], :NV => [X0_NV]) #Dict of species concentrations
#kinetics for Nitrobacter vulgaris
VmNV = 3.144E-13
KmNV = 0.049
Vm_NV = (VmNV*1000)/wt_NB #(mmol/gDW-d)
#dynamics
Δt = 0.01 #time step in days
n = 100 #number of timesteps
θ = Δt*n #bioreactor residence time (days)
#store FBA model outputs
out = Dict( :NS => [], :NB => [] , :NV => [] )
for j in 1:n
#nutrient status at time t
nut = Dict(:NH3=> C[:NH3][j],:NO2 => C[:NO2][j],:NO3 => C[:NO3][j] )
#step 1A: solve FBA model for Nitrosomonas
if C[:NH3][j] <= 0
uptakeNS = Dict("EX_NH3" => 0 )
push!(out[:NS],fbaNS(uptakeNS,nut))
else
uptakeNS = Dict("EX_NH3" => Vm_NS*( C[:NH3][j]/ (C[:NH3][j] + KmNS )) )
push!(out[:NS],fbaNS(uptakeNS,nut))
end
#step 1B: solve FBA model for Nitrobacter
if C[:NO2][j] <= 0
uptakeNB = Dict("EX_NO2" => 0 )
uptakeNV = Dict("EX_NO2" => 0 )
push!(out[:NB],fbaNB(uptakeNB,nut))
push!(out[:NV],fbaNB(uptakeNV, nut))
else
uptakeNB = Dict("EX_NO2" => Vm_NB*( C[:NO2][j]/(C[:NO2][j] + KmNB )) )
uptakeNV = Dict("EX_NO2" => Vm_NV*( C[:NO2][j]/(C[:NO2][j] + KmNV )) )
push!(out[:NB],fbaNB(uptakeNB, nut))
push!(out[:NV],fbaNB(uptakeNV, nut))
end
#step 2: mass balance around biomass (X)
X[:NS] = [X[:NS]; X[:NS][j] + MW*out[:NS][j][1]*X[:NS][j]*Δt - KdNS*X[:NS][j]*Δt ]
X[:NB] = [X[:NB]; X[:NB][j] + MW*out[:NB][j][1]*X[:NB][j]*Δt - KdNB*X[:NB][j]*Δt ]
X[:NV] = [X[:NV]; X[:NV][j] + MW*out[:NV][j][1]*X[:NV][j]*Δt - KdNB*X[:NV][j]*Δt ]
#step 3: mass balance around nutrients (C)
C[:NH3] = [C[:NH3]; C[:NH3][j] + out[:NS][j][2]*X[:NS][j]*Δt + out[:NB][j][2]*X[:NB][j]*Δt + out[:NV][j][2]*X[:NV][j]*Δt ]
C[:NO2] = [C[:NO2]; C[:NO2][j] + out[:NS][j][3]*X[:NS][j]*Δt + out[:NB][j][3]*X[:NB][j]*Δt + out[:NV][j][3]*X[:NV][j]*Δt ]
C[:NO3] = [C[:NO3]; C[:NO3][j] + out[:NS][j][4]*X[:NS][j]*Δt + out[:NB][j][4]*X[:NB][j]*Δt + out[:NV][j][4]*X[:NV][j]*Δt ]
end
```

Table 2 provides the steady-state flux values (mmol/gDW/day) for biomass growth and nutrient uptake by *Nitrosomonas* and *Nitrobacter*. These results agree well with what we would expect for chemolithoautotrophic growth. For *Nitrosomonas*, O_{2}, CO_{2}, and NH_{3} are all transported into the cell to generate biomass and NO_{2}^{-}. The maximum specific growth rate for *Nitrosomonas* can be calculated by multiplying the molecular weight of a cell (1.198 gDW/mmol) by the biomass flux. This results in a growth rate of 12.9 days^{-1} (1.85 hour doubling time) for *Nitrosomonas*, which is approximately 3 times greater than expected (typ. doubling time is ~7 hours).

For *Nitrobacter*, O_{2}, CO_{2}, and NO_{2}^{-} are all transported into the cell to generate biomass and NO_{3}^{-}, as expected. We also note that when the extracellular concentration of NH_{3} is zero (a key metabolite for biomass synthesis), *Nitrobacter* is still capable of growth because ammonia can be produced by reducing nitrite via the enzyme nitrite reductase (reaction "R_nirBD"). The maximum specific growth rate calculated for *Nitrobacter* was 0.80 days^{-1} (doubling time of 29 hours), which agrees well with values reported in the literature (26 hours, Nowka et al., 2015).

Reaction | Nitrosomonas Flux |
Nitrobacter Flux |
---|---|---|

Biomass | 10.8 | 0.67 |

EX_O2 | -945.1 | -102.9 |

EX_NH3 | -713.3 | 0.0 |

EX_NO2 | 702.4 | -221.9 |

EX_NO3 | 0.0 | 221.3 |

EX_H | 425.7 | 25.7 |

EX_NAD | -21.7 | -1.3 |

EX_CO2 | -54.2 | -3.3 |

EX_H2O | 759.1 | 2.8 |

Figure 6 and 7 show the shadow prices for all metabolites during maximum growth of *Nitrosomonas* and *Nitrospira*. For *Nitrosomonas*, we see biomass growth is most sensitive to the metabolites protein, hydroxylamine (NH_{2}OH), NADH (reducing power), and ammonia (NH_{3}), where *Nitrobacter* biomass growth was most sensitive to protein, NADH, and ammonia. Increasing the intracellular amounts of these metabolites would improve the observed maximum growth rate of these organisms, and are possible targets for metabolic engineering efforts.

In [7]:

```
#plot shadow prices
using PyPlot
dualNS = []
for i in metabolitesNS
push!(dualNS,λNS[i])
end
x = collect(1:length(metabolitesNS));
plot(x, dualNS, "ro", x, zeros(length(metabolitesNS)), "--")
xticks(x, metabolitesNS, rotation="vertical")
xlabel("Metabolites"); ylabel("Shadow price")
title("Figure 6 - Nitrosomonas Shadow Prices")
```

Out[7]:

In [8]:

```
#plot shadow prices
using PyPlot
dualNB = []
for i in metabolitesNB
push!(dualNB,λNB[i])
end
x = collect(1:length(metabolitesNB))
plot(x, dualNB, "ro", x, zeros(length(metabolitesNB)), "--")
xticks(x, metabolitesNB, rotation="vertical")
xlabel("Metabolites"); ylabel("Shadow price")
title("Figure 7 - Nitrobacter Shadow Prices")
```

Out[8]:

Figure 8 shows the x-y phenotypic phase plane for *Nitrosomonas* under varying fluxes for O_{2} and NH_{3} uptake. The line of optimality (red) corresponds to where the optimal oxygen and ammonia uptake fluxes provide the maximum biomass yield. From the slope of the isoclines, we can see that two distinct phenotypic phases or "metabolic network states" exist:

- Phase 1: Ammonia is the single limiting substrate (left portion of phase plane, vertical isoclines).
- Phase 2: Oxygen is the single limiting substrate (right portion of phase plane, horizontal isoclines).

The two phenotypic phases can also be observed from the shadow prices shown on Table 2, where both metabolites have a constant value in a given phase. In the ammonia limiting phase, we see that the objective function (biomass growth) is sensitive to changes in NH_{3}, but not O_{2}, whereas the reverse case holds for points in the oxygen limiting phases.

Metabolite | λ(520:100) | λ(590:340) | λ(200:40) | λ(300:200) | λ(300:85) | λ(350:645) | λ(55:530) | λ(30:200) | λ(185:300) | λ(590:590) |
---|---|---|---|---|---|---|---|---|---|---|

NH_{3} |
-0.01521 | -0.01521 | -0.01521 | -0.01521 | -0.01521 | 0 | 0 | 0 | 0 | 0 |

O_{2} |
0 | 0 | 0 | 0 | 0 | -0.01148 | -0.01148 | -0.01148 | -0.01148 | -0.01148 |

Phase | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 2 |

Because the metabolic network used in this analysis was very small (30 reactions), no alternate phenotypic phases were found where both substrates became limiting. However, we anticipate that for larger metabolic reconstuctions of *Nitrosomonas* (e.g. genome-scale model) multiple phenotypic phases would be observed.

In [9]:

```
#plot PhPP
using PyPlot
function meshgrid{T}(b1::AbstractVector{T}, b2::AbstractVector{T})
m, n = length(b1), length(b2)
b1 = reshape(b1, 1, n)
b2 = reshape(b2, m, 1)
(repmat(b1, m, 1), repmat(b2, 1, n))
end
(X1,Y1) = meshgrid(b1,b2);
#LO line end
a = find(rxn2Flux .<= length(1:range1) )
n = (a[length(a)])
pygui(false)
figure(figsize=(6,5))
contour(growthRates, origin="lower", extent=(0,range1,0,range2), 60, )
plot( b1[1:n], rxn2Flux[1:n], "r-" )
axis("image"); xlabel(rxnID[1]); ylabel(rxnID[2]);
colorbar()
tight_layout();
savefig("PhPP_contour.pdf")
title("Figure 8A - Phenotypic Phase Plane")
```

Out[9]:

In [10]:

```
#plot PhPP
pygui(false)
fig = surf(X1,Y1, growthRates+10, rstride=5, cstride=5,cmap="rainbow")
contour(X1, Y1, growthRates, zdir="z", 50, offset=0, origin="lower" )
#plot( b1[1:n], rxn2Flux[1:n],growthRatesLO[1:n], "r-" )
xlabel("NH3 uptake rate"); ylabel("O2 uptake rate"); zlabel("Growth rate")
tight_layout()
savefig("PhPP.pdf")
title("Figure 8B - 3D Phase Plane")
```

Out[10]:

The ultimate goal of our project was to be able to simulate the population dynamics and nutrient cycling occuring in a nitrifying bioreactor. Figure 9 shows the nutrient profile simulated over one bioreactor batch cycle. As anticipated for nitrification, NH_{3} is initially consumed followed by a transient accumulation of NO_{2}^{-} and finally NO_{3}^{-} production. Complete nitrification is achieved at a hydraulic retention time of approximately 15 hours, when the biomass concentrations are initially at 1x10^{10} cells/L.

From the biomass concentrations shown on Figure 10, we can see that *N. vulgaris* was able to outcome N. winogradskyi for NO_{2}^{-} and proliferate in the system. This is expected becuase both its maximum uptake rate and affinity coefficient for NO_{2}^{-} are better.

A major limitation of the bioreactor model and FBA analyses done here is they do not account for the effects of toxic intermediates. For example, hydroxylamine, nitric oxide, nitrous acid, and ammonia at high concentrations are known to inhibit the growth of nitrifying bacteria (Vadivelu et al., 2006). This could be corrected in the model by adding kinetic terms that explicitly account for toxicity effects.

In [11]:

```
#plot nutrient concentrations
using PyPlot
x = collect(0:Δt:θ)*24
figure(figsize=(12,4))
plot(x, C[:NH3],"-", x, C[:NO2],"-", x, C[:NO3] )
xlabel("time (hrs)"); ylabel("Concentration (mM)")
legend(["NH3", "NO2", "NO3"], loc="right")
title("Figure 9 - Bioreactor nutrient profile");
```

In [12]:

```
#plot biomass growth
x = collect(0:Δt:θ)*24
figure(figsize=(12,4))
plot(x, X[:NS]*1000, "--", x, X[:NB]*1000, "--", x, X[:NV]*1000, "--")
xlabel("time (hrs)"); ylabel("conc. (mg DW/L)")
legend(["NS", "NB", "NV"], loc="lower right");
title("Figure 10 - Biomass growth");
```

The following conclusions can be drawn from this report:

- The growth phenotype of an organism can be predicted from it's metabolic network using linear programming.
- The maximum specific growth rates for
*Nitrosomonas*and*Nitrobacter*agree reasonably well with data from the literature. - The metabolite shadow prices from the dual FBA problem can be used to evalute the sensitivity of an organisms maximum growth rate to metabolite imbalances.
- Metabolite shadow prices can also be used to examine the different phenotypes of an organism via phenotypic phase plane analysis.
- Dynamics FBA modelling and mathamatical optimization can be used to predict environmental processes, such as nitrification, that rely on mixed-species microbial communities.
*N. vulgaris*outcompetes*N. winogradskyi*for NO_{2}-^{}under the conditions simulated in this study.

Future work will include the addition of more cooperative species to complete the nitrogen cycle. In particular, the addition of a denitrifying bacteria, such as *Thiobacillus dentrificans*, *Micrococcus denitrificans*, or *Pseudomonas* species would allow for complete conversion of NH_{3} to nitrogen gas (N_{2}) in our bioreactor system.

Another promising future direction is the application of genome-scale metabolic models. The ~25-30 reaction models used here are a known oversimplication of the actual metabolic network. We predict that incorporating a more complete metabolic reconstruction will reveal alternative phenotypic phases, and allow for a more accurate and comprehensive analysis of microbial metabolic networks.

Finally, we may continue to investigate bi-level dynamic optimzation of the microbial system. Rather than modeling the extracellular enviroment as an ordinary differential equation network, the bioreactor dynamics may also be framed as an optimization problem (see Zomorrodi et al., 2014). This formulation seek to maximize the global biomass concentration subject to constraints generated by "inner" flux balance analysis problems, and represents a nonconvex nonlinear problem.