Team:MIT/Model

Model

Goal and Motivation

Theoretically, when our functional BCAA strain is grown in vitro, for instance in a broth with some arbitrary BCAA concentration, we expect that the BCAA concentration decreases predictably with respect to time given by some rate law. Our goal is to derive this rate law from our gene circuit. If we have a reasonable rate law for BCAA concentrations in vitro, we hope that it can be extended to the human gut microbiome. In addition, knowing a rate law will help with calculating rate constants from empirical data. Lastly, knowing the rate law would help us optimize our design to maximize the rate of BCAA disappearance. Below we derive our rate law by breaking it down into separate models that we combine (using algebraic substitution) to get the rate law for BCAA disappearance.

Multicompartment Modeling

First, let us establish how mass action kinetics work in multi-compartment systems. We know that mass must be conserved, so the number of moles disappearing from one compartment is the number of moles appearing in another compartment, and when accounting for volume: $$\frac{dA_{1}}{dt} = - \frac{dA_{2}}{dt}$$ $$\frac{dA_{1}}{dt} \cdot \frac{Vol_1}{Vol_1} = - \frac{dA_{2}}{dt} \cdot \frac{Vol_2}{Vol_2}$$ $$\frac{dA_{1}}{dt} \cdot \frac{Vol_1}{Vol_1} = - \frac{dA_{2}}{dt} \cdot \frac{Vol_2}{Vol_2}$$ $$\frac{d[A_{1}]}{dt} \cdot Vol_1 = - \frac{d[A_{2}]}{dt} \cdot Vol_2$$ $$\frac{d[A_{1}]}{dt} = - \frac{d[A_{2}]}{dt} \cdot \frac{Vol_2}{Vol_1}$$ This essentially models moles of a molecule moving from a compartment to another. However, the above equations must be adjusted stoichiometrically if it is a chemical reaction as opposed to moles moving between compartments. Also, for relating rate laws to multi-compartment modeling, the volume of the compartment is what matters: $$\frac{dA_{1}}{dt} = k[B_2][C_2]$$ $$\frac{dA_{1}}{dt} \cdot \frac{Vol_1}{Vol_1} = k[B_2][C_2]$$ $$\frac{d[A_{1}]}{dt} \cdot Vol_1 = k[B_2][C_2]$$ In the following sections, we need to be aware of the fact that essentially we have three compartments; the volume of broth in a flask (which we model as constant), and the total volume of activated B. subtilis cells, and the total volume of inactive B. subtilis cells. Though the cells are dispersed in solution, they can be modelled as a mass of cells with some volume.

Recombinase Modeling

We reiterate the need for a recombinase switch: constitutive expression of BCKDH, a large bulky proteins, will certainly harm cell fitness and growth. Therefore, our gene circuit should only be activated before the consumption of the probiotic. We have a constitutive promoter with the gene tetR. tetR is a repressor for the tet promoter of the Cre gene. In this repressed state, the Cre-Lox system does nothing since Cre cannot be synthesized. Once we add a tetracycline derivative (e.g DOX/ATC), the system is activated because DOX/ATC will bind to tetR proteins such that they can no longer function as a repressor, and once Cre synthesis is no longer be repressed, synthesis of Cre begins. Cre "scans" the B. subtilis circular chromosome for the LOX site, cutting out the the terminator and endogenous promoter in a recombination reaction. $$\text{DOX} + \text{TetR} \xrightarrow{k_5} \text{TetR}_{inactive}$$ $$\frac{d[\text{tetR}]}{dt} = \gamma_2 - \varsigma_2 \cdot [\text{tetR}] - k_5[DOX][TetR]$$ $$\frac{d[Cre_{mRNA}]}{dt} = k_6 \cdot (1 - \frac{[tetR]^n}{K_d^n + [tetR]^n}) - \varsigma_3 [Cre_{mRNA}]$$ $$\frac{d[Cre]}{dt} = k_7 \cdot [Cre_{mRNA}] - \varsigma_4 \cdot [Cre]$$ Notice that, when DOX is not present, there is an abundance of TetR, such that and almost no Cre can be synthesized. We proceed to calculate the proportion of the bacterial population that is in the activated state versus the inactivated state (the activated state being bacteria which now constitutively express BCKDH and other gene circuit due to excision of a terminator). In a work by Ringrose et. al., an ODE model was developed that describes Cre recombination reaction. Essentially, there is a multistage process where Cre monomers bind to LoxP sites, forming metastable intermediates that procced to the end product of excising the DNA sequence. Interestingly, the recombination fraction as a function of time asymptotically approaches 60%, and increasing Cre concentration does not increase the maximum recombination portion. As a result of a large number of intermediate states, there are many terms in Ringrose's ODE model, so the 2017 HKSUT developed a simplified model that is much less computationally expensive. Let \(S\) be the proportion of DNA sequences without excised sites, and \(P\) be the proportion of DNA sequences with excised sites. Then the recombination reaction can be characterized by: $$\text{S} \xrightarrow{k_8} \text{P}$$ $$\text{P} \xrightarrow{k_9} \text{S} $$ $$\frac{dP}{dt} = k_8 S - k_9 P$$ $$\frac{dS}{dt} = k_9 P - k_8 S$$ However, the rate of the recombination reaction depends on the concentration of Cre. In the 2017 HKSUT recombinase model, the concentration of Cre was constant, so which is why \(k_5\) and \(k_6\) are constants. For our model, we rewrite the equation taking into account the concentration of Cre: $$\text{Cre} + P_{inactive} \xrightarrow{k_{10}} P_{active} $$ $$\text{Cre} + P_{active} \xrightarrow{k_{11}} P_{inactive} $$ $$\frac{dP_{active}}{dt} = k_{10} [Cre] P_{inactive} - k_{11} [Cre] P_{active}$$ $$\frac{dP_{inactive}}{dt} = k_{11} [Cre] P_{active} - k_{10} [Cre] P_{inactive}$$ Notice that increasing the concentration of Cre does not change the total portion of DNA excised at equilibrium, but only the rate that it reaches equilibrium.

Gene Expression and Protein Synthesis

In B. subtilis cells that are activated, since all of the genes of interest are constitutively expressed, we can model their rates by the following ODES, where ilvK encodes for transaminase, bCAP encodes for BCAA transport, and BCKDH encodes for the enzyme complex that breaks down BCKAs: $$\frac{d[\text{ilvK}]}{dt} = \gamma_{1} - \varsigma_{1} \cdot [\text{ilvK}]$$ $$\frac{d[\text{bCAP}]}{dt} = \gamma_{1} - \varsigma_{1} \cdot [\text{bCAP}]$$ $$\frac{d[\text{BCKDH}]}{dt} = \gamma_{1} - \varsigma_{1} \cdot [\text{BCKDH}]$$ \(\gamma\) represents the respective basal translation rates, and \(\varsigma\) accounts for dilution of the proteins due to cellular division and degradation. The \(\gamma\) and \(\varsigma\) for each of the proteins will be slightly different even though they have the same promoter, but we make a simplifying assumption that they equal. In inactive B. subtilis cells, the concentration of \([ilvK], [bCAP], [BCKHD]\) are relatively low and constant.

Enzyme Kinetics

The breakdown of all three branched chain amino acids follow a similar pathway: $$\text{BCAA}_{in} + \text{oxoglutarate} + \text{ilvK} ^{\xrightarrow[k_{-1}]{k_{1}}}_{\leftarrow} complex ^{\xrightarrow[k_-2]{k_2}}_{1 \leftarrow} \text{BCKA} + \text{ilvK} + \text{glutamate}$$ $$\text{BCKA} + \text{BCKDH} + NAD^{+} + \text{CoA} ^{\xrightarrow[k_{-3}]{k_{3}}}_{\leftarrow} complex_2 \xrightarrow{k_4} \text{butanol-CoA} + \text{BCKDH} + \text{CO}_2 + \text{NADH}$$ Transamination is readily reversible while the reaction catalyzed by the BCKDH complex is irreversible. We do not concern ourselves with the breakdown of butanol-COA since we assume that it is quickly broken down by the B. subtilis cell as part of the cellular respiration pathway. Also, we can assume that the concentration of NAD+, oxoglutarate, and CoA are relatively constant as they are replenished by other pathways in the cell, as these are important reactants of cellular respiration and should be quickly used up. The enzyme kinetics in isolation look like this: $$\frac{d[BCAA_{in}]}{dt} = -k_1[BCAA][ilvK][oxoglutarate] + k_{-1}[complex_1]$$ $$\frac{d[ilvK]}{dt} = -k_1[BCAA][ilvK][oxoglutarate] + k_{-1}[complex_1] + k_2[complex_1]$$ $$\frac{d[complex_1]}{dt} = k_1[BCAA][ilvK][oxoglutarate] - k_{-1}[complex_1] - k_2[complex_1] + k_{-2}[BCKA][ilvK][glutamate]$$ $$\frac{d[BCKA]}{dt} = k_2[complex_1] - k_{-2}[BCKA][ilvK][glutamate] - k_3[BCKA][BCKDH][NAD+][CoA] + k_{-3}[complex_2]$$ $$\frac{d[complex_2]}{dt} = k_3[BCKA][BCKDH][NAD+][CoA] - k_{-3}[complex_2] - k_{4}[complex_2]$$ $$\frac{d[BCKDH]}{dt} = -k_3[BCKA][BCKDH][NAD+][CoA] + k_{-3}[complex_2] + k_4[complex_2]$$ $$\frac{d[oxoglutarate]}{dt} = \frac{d[NAD+]}{dt} = \frac{d[CoA]}{dt} = 0$$

Amino Acid Transport

Transport of BCAA in the cell can be characterized as an enzymatic reaction via a symporter; the symporter encoded by BcaP is a likely a sodium cotransporter. BCAAs cannot diffuse through the cellular membrane because they are large and polar. We assume the reverse reaction is impossible since it is likely that the binding sites for BCAAs and NA is on the portion of the symporter outside of the cell. We assume the NA concentration inside and outside the cell is constant due to the activity of NA+/K+ -ATPase; experimental study has also shown that NA+ is not a rate limiting factor. $$\text{BCAA}_{out} + \text{Na}^{+}_{out} + bCAP ^{\xrightarrow[k_{b_{-1}}]{k_{b_1}}}_{\leftarrow} \text{bCAP}_{complex} \xrightarrow{k_{b_2}} \text{Na}^{+}_{out} + BCAA_{in} + bCAP$$ $$\frac{d[BCAA_{out}]}{dt} = (-k_{b_1}[BCAA_{out}][NA^{+}_{out}][bCAP] + k_{b_{-1}}[bCAP_{complex}]) \cdot \frac{1}{V_{broth}}$$ $$\frac{d[BCAA_{in}]}{dt} = (k_{b_2}[bCAP_{complex}]) \cdot \frac{1}{V_{active}}$$ $$\frac{d[NA^+_{in}]}{dt} = \frac{d[NA^+_{out}]}{dt} = 0$$ $$\frac{d[bCAP]}{dt} = (-k_{b_1}[BCAA_{out}][NA^{+}_{out}][bCAP] + k_{b_{-1}}[BCAP_{complex}] + k_{b_2}[bCAP_{complex}]) \cdot \frac{1}{V_{active}} $$ $$\frac{d[bCAP_{complex}]}{dt} = (k_{b_1}[BCAA_{out}][NA^{+}_{out}][bCAP] - k_{b_{-1}}[bCAP_{complex}] - k_{b_2}[bCAP_{complex}]) \cdot \frac{1}{V_{active}} $$

B. subtilis Population and Volume

When our B. subtilis cells are introduced to an in-vitro broth, the population of cells are not constant. The two most famous formulas for modeling microbial growth are the Monod equation and logistic growth curve, which are classic growth models in many textbooks. For the Monod equation, the population's specific growth rate is modeled as a function of some limiting substrate: $$\frac{dP}{dt} = \mu \cdot P_0$$ $$\mu = \mu_{max} \cdot \frac{[S]}{K_s+[S]}$$ However, for our model, we are not concerned with growth as a function of substrate concentration but rather population of B. subtilis as a function of time, so we use the logistic growth curve. For the logistic growth, the ODE can used to derive the population as a function of time: $$\frac{dP}{dt} = rP\cdot(1-\frac{P}{K}) $$ $$P(t) = \frac{K\cdot P_0}{(K-P_0)e^{rt}+P_0}$$ Where \(K\) is the asymptote, \(P_0\) is the initial population, and \(r\) is a constant. \(K\) and \(r\) depend on the substrate and temperature. The volume of a bacterial cell follows a normal distribution. For high cell counts, the total volume as a function of time is the product of the population and the average volume of a cell, \(v_{cell}\), such that: $$\frac{dV}{dt} = v_{cell} \cdot \frac{dP}{dt}$$

Of note, the total change in population for the inactive and active B. subtilis cells depends on recombinase activity as well. The growth rate for the inactivated B. subtilis cells will follow typical growth curves, while the growth rate for activated B. subtilis cells is likely very slow or negative due to decreased cell fitness. We assume that the population of active B. subtilis cells is non-diving and only changes due to recombinase activity. From our B. subtilis growth assay, we found \( K = 4.5 \cdot 10^7 \) B. sub cells per mL, \(r= -0.019 \), with initial condition \(P_0 = 3.8 \cdot 10^6\) B. sub cells per mL.

Complete Model

Combining all the equations from pervious sections into one system, we arrive at a complete model. The term \( \frac{d[BCAA_{out}]}{dt} \) is the rate of disappearance of BCAAs from the in-vitro environment, i.e. LB broth, that the B. subtilis cells are growing in, which is what we were trying to derive all along.

$$\frac{dP_{active}}{dt} = r_1 \cdot P_{active}\cdot(1-\frac{P_{active}}{K_1}) + (k_7 [Cre]P_{inactive} - k_8[Cre]P_{active} )$$ $$\frac{dP_{inactive}}{dt} = k_8 [Cre]P_{active} - k_7 [Cre]P_{inactive} $$ $$V_{active} = v_{cell} \cdot P_{active}$$ $$V_{inactive} = v_{cell} \cdot P_{inactive}$$ $$V_{total} = V_{active} + V_{inactive}$$ $$\frac{d[\text{tetR}]}{dt} = (\gamma_2 - \varsigma_2 \cdot [\text{tetR}] - k_5[DOX][TetR]) \cdot \frac{1}{V_{total}}$$ $$\frac{d[Cre_{mRNA}]}{dt} = (k_6 \cdot (1 - \frac{[tetR]^n}{K_d^n + [tetR]^n}) - \varsigma_3 [Cre_{mRNA}]) \cdot \frac{1}{V_{total}}$$ $$\frac{d[Cre]}{dt} = (k_7 \cdot [Cre_{mRNA}] - \varsigma_4 \cdot [Cre]) \cdot \frac{1}{V_{total}} $$ $$\frac{d[BCAA_{out}]}{dt} = (-k_{b_1}[BCAA_{out}][NA^{+}_{out}][bCAP] + k_{b_{-1}}[bCAP_{complex}]) \cdot \frac{1}{V_{broth}}$$ $$\frac{d[bCAP]}{dt} = (-k_{b_1}[BCAA_{out}][NA^{+}_{out}][bCAP] + k_{b_{-1}}[bCAP_{complex}] + k_{b_2}[bCAP_{complex}] + (\gamma_1 - \varsigma_1 \cdot [\text{bCAP}])) \cdot \frac{1}{V_{active}}$$ $$\frac{d[bCAP_{complex}]}{dt} = (k_{b_1}[BCAA_{out}][NA^{+}_{out}][bCAP] - k_{b_{-1}}[bCAP_{complex}] - k_{b_2}[bCAP_{complex}]) \cdot \frac{1}{V_{active}}$$ $$\frac{d[BCAA_{in}]}{dt} = (-k_1[BCAA][ilvK][oxoglutarate] + k_{-1}[complex_1] + k_{b_2}[bCAP_{complex}]) \cdot \frac{1}{V_{active}}$$ $$\frac{d[ilvK]}{dt} = (-k_1[BCAA][ilvK][oxoglutarate] + k_{-1}[complex_1] + k_2[complex_1] + (\gamma_1 - \varsigma_1 \cdot [\text{ilvK}])) \cdot \frac{1}{V_{active}} $$ $$\frac{d[complex_1]}{dt} = (k_1[BCAA_{in}][ilvK][oxoglutarate] - k_{-1}[complex_1] - k_2[complex_1] + k_{-2}[BCKA][ilvK][glutamate]) \cdot \frac{1}{V_{active}}$$ $$\frac{d[BCKA]}{dt} = (k_2[complex_1] - k_{-2}[BCKA][ilvK][glutamate] - k_3[BCKA][BCKDH][NAD+][CoA] + k_{-3}[complex_2]) \cdot \frac{1}{V_{active}}$$ $$\frac{d[complex_2]}{dt} = (k_3[BCKA][BCKDH][NAD+][CoA] - k_{-3}[complex_2] - k_{4}[complex_2]) \cdot \frac{1}{V_{active}}$$ $$\frac{d[BCKDH]}{dt} = (-k_3[BCKA][BCKDH][NAD+][CoA] + k_{-3}[complex_2] + k_4[complex_2] + (\gamma_1 - \varsigma_1 \cdot [\text{BCKDH}])) \cdot \frac{1}{V_{active}}$$ $$\frac{d[NA^+_{in}]}{dt}, \frac{d[NA^+_{in}]}{dt}, \frac{d[NA^+_{out}]}{dt}, \frac{d[oxoglutarate]}{dt}, \frac{d[NAD+]}{dt}, \frac{d[CoA]}{dt} = 0$$

Results

We used our theoretical model to inform our design choices, since from the equations we were able to see how our choices would impact the rate of disappearance of BCAAs based on the rate law. Attempting to numerically solve the ODEs with estimated rate constants produced the below graph:

We see that the rate of disappearance of BCAAs is low at first, but increases once the recombinase begins converting inactive cells to active cells, and once the cells begin synthesizing catabolic enzymes. Of note, our parts are not well characterized in literature, so our estimates for rate constants are error prone. We would like to proceed with experiments that could derive more accurate rate constants, and produce experimental results that could refine the model.

References

  1. de Jong, H., Casagranda, S., Giordano, N., Cinquemani, E., Ropers, D., Geiselmann, J., & Gouzé, J.-L. (2017). Mathematical modelling of microbes: metabolism, gene expression and growth. J. R. Soc. Interface., 14(136), 20170502. 10.1098/rsif.2017.0502
  2. den Hengst, C. D., Groeneveld, M., Kuipers, O. P., & Kok, J. (2006). Identification and Functional Characterization of the Lactococcus lactis CodY-Regulated Branched-Chain Amino Acid Permease BcaP (CtrA). J Bacteriol, 188(9), 3280–3289. 10.1128/jb.188.9.3280-3289.2006
  3. HKUST 2017 iGEM Team. (2017). Team:Hong Kong HKUST/Model/Recombination - 2017.igem.org. Single Cell ODE Model for the Recombination Module. https://2017.igem.org/Team:Hong_Kong_HKUST/Model/Recombination
  4. Hofmeyr, J.-H. S. (2020). Kinetic modelling of compartmentalised reaction networks. 10.31219/osf.io/5as9h
  5. McLellan, M. A., Rosenthal, N. A., & Pinto, A. R. (2017). Cre-loxP-Mediated Recombination: General Principles and Experimental Considerations. Current Protocols in Mouse Biology, 7(1), 1–12. 10.1002/cpmo.22
  6. Ringrose, L., Lounnas, V., Ehrlich, L., Buchholz, F., Wade, R., & Stewart, A. F. (1998). Comparative kinetic analysis of FLP and cre recombinases: mathematical models for DNA binding and recombination. Journal of Molecular Biology, 284(2), 363–384. 10.1006/jmbi.1998.2149
  7. Saier Lab. (n.d.). Branched chain amino acid (Leucine/isoleucine/valine) uptake transporter of 469 aas and 12 TMSs, BcaP or CitA. Retrieved October 21, 2021, from https://www.tcdb.org/search/result.php?tc=2.A.3.3.23
  8. Souba, W. W., & Pacitti, A. J. (1992). Review: How Amino Acids Get Into Cells: Mechanisms, Models, Menus, and Mediators. JPEN J Parenter Enteral Nutr, 16(6), 569–578. 10.1177/0148607192016006569