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. In our recombinase diagram, 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.