# Model

## Overview

We aim to develop mathematical models that act as a bridge between the theoretical and the physical realization of our biological work, to further optimize our experiments and understand our results.

We demonstrated that our E. coli can convert tyrosine - the precursor of p-Cresol, into p-Coumaric acid, a beneficial byproduct by diverting the original fermentation pathway. For that, we simulated the interactions between E. coli Nissle and Clostridium difficile (C. difficile) under co-culture conditions.

Our Goals:

1. Better understand and improve our CreSolve therapeutic drug’s effects on reducing p-Cresol production.

2. Explore the effectiveness of bacteriocin on reducing Clostridium difficile population.

## Co-culture Model

Oh My Gut’s CreSolve is engineered as a living therapeutic and so, it is vital to simulate our therapeutic drug’s efficiency in reducing excess tyrosine. However, due to time constraints, we were not able to perform a co-culture experiment to further validate our engineered E. coli Nissle. Therefore, to prove the feasibility of CreSolve, we built a series of co-culture model under two different conditions:

1. C. difficile and CreSolve are cultured under the same initial OD value of 0.06 (1 : 1).

2. CreSolve cultured in a C. difficile-dominant environment (1:400).

The following equations are used in the co-culture models:

Table 1. Co-culture model equations
Equations Description
$\frac{\mathit{d}\left[\mathit{TAL}\right]}{\mathit{d}t}={k}_{1}⁢\left[\mathit{Nissle}\right]$ We assumed that the rate of change of TAL (tyrosine ammonia-lyase) concentration is proportional to E. coli Nissle cell density. k1 stands for a coefficient.
$\frac{\mathit{d}{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}{\mathit{d}t}=-\frac{{k}_{\mathit{cat1}}⁢{\left[\mathit{TryP}\right]}_{\mathit{outside}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}{{K}_{\mathit{m1}}+{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}$ Transportation of tyrosine into the cell through tyrosine transporter according to Michaelis-Menten equation. [TryP]outside stands for the tyrosine transporter concentration; kcat1 is the estimated turnover number of tyrosine transporter; Km1 is the tyrosine concentration at which tyrosine transportation rate is at half-maximum.
$\frac{\mathit{d}{\left[\mathit{Tyr}\right]}_{\mathit{Nissle}}}{\mathit{d}t}=\frac{{k}_{\mathit{cat1}}⁢{\left[\mathit{TryP}\right]}_{\mathit{cell}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}{{K}_{\mathit{m1}}+{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}-\frac{{k}_{\mathit{cat2}}⁢{\left[\mathit{TAL}\right]}_{\mathit{Nissle}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{Nissle}}}{{K}_{\mathit{m2}}+{\left[\mathit{Tyr}\right]}_{\mathit{Nissle}}}$ Change of tyrosine concentration in E. coli Nissle over time based on the transportation of tyrosine into the cell and the conversion of tyrosine into p-Coumaric acid by TAL. Both tyrosine transportation and p-Coumaric production are based on Michaelis-Menten equation. kcat2 is the estimated turnover number of TAL; Km2 is the tyrosine concentration at which tyrosine fermentation rate by TAL is at half-maximum.
$\frac{\mathit{d}{\left[\mathit{PCA}\right]}_{\mathit{tot}}}{\mathit{d}t}=\frac{{k}_{\mathit{cat2}}⁢{\left[\mathit{TAL}\right]}_{\mathit{tot}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{Nissle}}}{{K}_{\mathit{m2}}+{\left[\mathit{Tyr}\right]}_{\mathit{Nissle}}}$ p-Coumaric production rate after tyrosine is imported through tyrosine transporter and converted by TAL. Based on Michaelis-Menten equation.
$\frac{\mathit{d}\left[\mathit{Nissle}\right]}{\mathit{d}t}={\mu }_{1}\cdot \left[\mathit{Nissle}\right]\frac{1-\left(\left[\mathit{Nissle}\right]+{\alpha }_{1}\left[\mathit{CD}\right]\right)}{{K}_{1}}$ Rate of E. coli Nissle cell density increase based on the competitive Lotka-Volterra equation. μ1 stands for the specific growth rate of E. coli Nissle, and K1 stands for carrying capacity. α1 represents the effect C. difficile has on E. coli Nissle.
$\frac{\mathit{d}\left[\mathit{PHPAD}\right]}{\mathit{d}t}={k}_{2}⁢\left[\mathit{CD}\right]$ We assumed that the rate of change of p-HPA decarboxylase (p-Hydroxyphenylacetate decarboxylase) concentration is proportional to C. difficile cell density. k2 stands for a coefficient.
$\frac{\mathit{d}{\left[\mathit{PHPAD}\right]}_{\mathit{CD}}}{\mathit{d}t}=\frac{{k}_{\mathit{cat1}}⁢{\left[\mathit{TryP}\right]}_{\mathit{cell}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}{{K}_{\mathit{m1}}+{\left[\mathit{Tyr}\right]}_{\mathit{outside}}}-\frac{{k}_{\mathit{cat3}}⁢{\left[\mathit{PHPAD}\right]}_{\mathit{CD}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{CD}}}{{K}_{\mathit{m3}}+{\left[\mathit{Tyr}\right]}_{\mathit{CD}}}$ Change of tyrosine concentration in C. difficile over time based on the transportation of tyrosine into the cell and the conversion of tyrosine into p-Cresol by p-HPA decarboxylase. Both tyrosine transportation and p-Cresol production are based on Michaelis-Menten equation. kcat3 is the estimated turnover number of p-HPA decarboxylase; Km3 is the tyrosine concentration at which tyrosine fermentation rate by p-HPA decarboxylase is at half-maximum.
$\frac{\mathit{d}{\left[\mathit{PC}\right]}_{\mathit{tot}}}{\mathit{d}t}=\frac{{k}_{\mathit{cat3}}⁢{\left[\mathit{PHPAD}\right]}_{\mathit{tot}}⁢{\left[\mathit{Tyr}\right]}_{\mathit{CD}}}{{K}_{\mathit{m3}}+{\left[\mathit{Tyr}\right]}_{\mathit{CD}}}$ p-Cresol production rate after tyrosine is imported through tyrosine transporter and converted by p-HPA decarboxylase. Based on Michaelis-Menten equation.
$\frac{\mathit{d}\left[\mathit{CD}\right]}{\mathit{d}t}={\mu }_{2}\cdot \left[\mathit{CD}\right]\frac{1-\left(\left[\mathit{CD}\right]+{\alpha }_{2}\left[\mathit{Nissle}\right]\right)}{{K}_{2}}$ Rate of C. difficile cell density increase based on the competitive Lotka-Volterra equation. μ2 stands for the specific growth rate of C. difficile, and K2 stands for carrying capacity. α2 represents the effect E. coli Nissle has on C. difficile.

The equations we used in the co-culture models are mostly based on Michaelis-Menten equation and competitive Lotka-Volterra equation. The Michaelis-Menten equation is commonly used in biochemistry to model enzyme kinetics. Under the assumption that enzyme concentration is much less than the substrate concentration, the rate of product formation is shown in the equation below:

$v = d P d t = V max ⁢ S K M + S = k cat ⁢ E 0 S K M + S$

$v$ is the rate of enzymatic reactions (rate of product formation, $\frac{\mathit{d}\left[\mathit{P}\right]}{\mathit{d}t}$). [S] is the concentration of a substrate. ${V}_{\mathrm{max}}$ represents the maximum reaction rate achieved by the system. The value of Michaelis constant ${K}_{m}$ is the substrate concentration at which the reaction rate reaches 1/2 ${V}_{\mathrm{max}}$. ${k}_{\mathrm{cat}}$, the turnover number, is the maximum number of substrate molecules converted to product per enzyme molecule per second.

Though TryP is not an enzyme, the kinetics of membrane transport can also be described by Michaelis-Menten equation.

The competitive Lotka-Volterra equation is similar to the logistic equations, which is a commonly used exponential population model, but with an additional term to account for the species' interactions:

$d x 1 d t = r 1 ⁢ x 1 1 - x 1 + α 12 ⁢ x 2 K 1$

$d x 2 d t = r 2 ⁢ x 2 1 - x 2 + α 21 ⁢ x 1 K 2$

$x$ is the population size of the two species. $r$ is the specific growth rate. $K$ is the carrying capacity, representing the maximum population size of the species that the environment can sustain indefinitely. The interaction between the two species is described by α12 and α21, which stands for the effect species 2 has on the population of species 1 and the effect species 1 has on the population of species 2, respectively.

These are the parameters we used:

Table 2. Co-culture model parameters
Parameters Description
k1= 0.01 Coefficient for the proportional relationship between the changing rate of TAL concentration and E. coli Nissle cell density.
k2= 0.01 Coefficient for the proportional relationship between the changing rate of p-HPA decarboxylase concentration and C. difficile cell density.
[TryP]outside=6.64
*10-4
TyrP transport protein accessible to growth medium.
[TryP]cell=0.151 TyrP transport protein expression in a cell.
kcat1=50 (s-1) The estimated turnover number of tyrosine transporter.
kcat2=0.015 (s-1) The estimated turnover number of TAL.
kcat3=2.7 (s-1) The estimated turnover number of p-HPA decarboxylase.
Km1=0.57 (μM) The tyrosine concentration at which tyrosine transportation rate is at half-maximum.
Km2=15.5 (μM) The tyrosine concentration at which tyrosine fermentation rate by TAL is at half-maximum.
Km3=2800 (μM) The tyrosine concentration at which tyrosine fermentation rate by p-HPA decarboxylase is at half-maximum.
α1=0.1 The effect C. difficile has on E. coli Nissle.
α2=0.1 The effect E. coli Nissle has on C. difficile.
μ1=0.82 The specific growth rate of E. coli Nissle.
μ2=0.13 The specific growth rate of C. difficile.
K1=0.62 The carrying capacity of E. coli Nissle.
K2=1.75 The carrying capacity of C. difficile.
[Tyr]Nissle=29 (μM) The initial tyrosine concentration in a growing E. coli cell.
[Tyr]CD=50 (μM) The initial tyrosine concentration in a growing C. difficile cell. We estimated the number to be larger than [Tyr]Nissle.
[TAL]cell The TAL concentration in an E. coli Nissle (CreSolve) cell.
[TAL]tot=0.0044
*[TAL]cell
The total TAL concentration both inside and outside E. coli Nissle (CreSolve). Vcell/Vtot=0.0044.
[PHPAD]cell The p-HPA decarboxylase concentration in a C. difficile cell.
The total p-HPA decarboxylase concentration both inside and outside a C. difficile cell. Vcell/Vtot=0.0044.
[Tyr]outside=1000(μM) The initial tyrosine concentration in the environment.

1. The models are used to simulate the in vitro co-culture experiment.
2. The expression of TAL is directly proportional to E. coli biomass.
3. Tyrosine can only be fermented to either p-Cresol or p-Coumaric acid.
4. p-Cresol can only be produced by C. difficile from tyrosine; p-Coumaric acid can only be produced by E. coli Nissle (CreSolve).
5. The tyrosine transporter in C. difficile and E. coli Nissle are the same (same function, same concentration in each cell, same transportation efficiency).
6. The tyrosine concentration in C. difficile ([Tyr]CD) is larger than E. coli Nissle ([Tyr]Nissle).
7. In the tyrosine fermentation pathway of C. difficile (Fig. 2), the rate-determining stage is the conversation of p-HPA into p-Cresol. The km3 and kcat3 values in equations 7 and 8 are determined by the enzyme p-HPA decarboxylase.
8. C. difficile and E. coli Nissle are grown under anaerobic conditions in BHI medium.
9. The p-Cresol level is not high enough to inhibit the growth of E. coli Nissle.
10. The effects of C. difficile and E. coli Nissle on each other (α1, α2) are the same.
11. CFU counting by OD value: 109 CFU/mL/OD600 nm.

### Co-culture model with the same OD value (1 : 1)

Under the first condition, C. difficile and E. coli Nissle are cultured under the same initial OD value of 0.06 (1 : 1), where the same initial amount of C. difficile and engineered E. coli Nissle (CreSolve) are used to simulate growth, tyrosine fermentation, p-Coumaric acid and p-Cresol productions. We started this co-culture model with an initial OD ratio of 1:1, to observe the interaction between C. difficile and engineered E. coli Nissle without additional factors. The growth model of both CreSolve and C. difficile co-culture, and C. difficile monoculture are then generated (Fig. 3). The results indicated that the addition of CreSolve slowed the growth of C. difficile. In addition, during co-culture, the overall tyrosine fermentation rate will decrease at a faster rate over time when compared with C. difficile monoculture alone (Fig. 4). From Fig. 5, the presence of CreSolve indeed lowered the production of p-Cresol by C. difficile. The tyrosine concentration decreases as expected (Fig. 4), and the decrease of p-Cresol is shown (Fig. 5). With the modeling results, we then proceed with a disparity population ratio between the two species to further explore CreSolve’s effectiveness. Fig. 3. Simulated growth rate of CreSolve and C. difficile at 1:1. Fig. 4. Tyrosine consumption by C. difficile monoculture compared with co-culture. Fig. 5. p-Cresol production rate by C. difficile monoculture compared with co-culture.

### Co-culture model with CreSolve cultured in a C. difficile-dominant environment (1:400)

In the second condition, we aim to know whether CreSolve is effective in a C. difficile-dominant environment, we assumed the ratio of 400 for C. difficile and 1 for the engineered E. coli Nissle (CreSolve). From Fig. 6, we can see that even though the initial CreSolve population is far less than C. difficile, it still has a significant effect on C. difficile growth. The tyrosine concentration decreases faster as expected (Fig. 7), and significant decrease of p-Cresol is shown (Fig. 8). With this model, we are more certain that CreSolve is effective in decreasing the p-Cresol level of CKD patients. Fig. 6. Simulated growth rate of C. difficile and CreSolve at 400:1. Fig. 7. Tyrosine consumption by C. difficile monoculture compared with co-culture (400:1). Fig. 8. p-Cresol production rate by C. difficile monoculture compared with co-culture (400:1).

## Bacteriocin Model

Since the Clostridium cluster is one of the major p-Cresol producers, besides changing the tyrosine fermentation pathway, we also aim to reduce the population of Clostridium spp. in the future. To show the killing efficiency of bacteriocin CBM-B, we modeled the relationship between C. difficile survival rate and CBM-B concentration.

As our model of the net growth rate of C. difficile population $\psi$, when exposed to CBM-B concentration $a$, we assumed the following relationship:

$ψ a = ψ max - μ a$

${\psi }_{\mathit{max}}$ is the growth rate of the C. difficile population in the absence of CBM-B, and $\mu \left(a\right)$ is the survival rate of the C. difficile population exposed to CBM-B concentration $a$, which we assume to be a Hill function:

$μ a = E max ⁢ a / EC 50 κ 1 + a / EC 50 κ$

The reason why we choose the Hill equation for this model is that according to our research, applications of bacteriocins are being tested to assess their application as a narrow-spectrum antibiotic, and pharmacodynamic models of antibiotics are often based on the Hill equation. Therefore, we consider it is suitable to use the Hill equation for our CBM-B effectivity model.

The Hill equation is commonly used in biochemistry, which refers to the binding of ligands to macromolecules. In pharmacology, it is extensively used to describe the drug concentration–effect relationship. In the equation above, ${E}_{\mathrm{max}}$ designates the maximum CBM-B-mediated survival rate, ${\mathrm{EC}}_{50}$is the CBM-B concentration at which the survival rate is at half of its maximum, ${E}_{\mathrm{max}}/2$, and $\kappa$ denotes the Hill coefficient, which determines the steepness of the sigmoid relationship between $\mu$ and $a$.

Hill coefficient describes the cooperativity of ligand binding in the following way:

• $\kappa >0$ (Positively cooperative binding): Once one ligand molecule is bound to the enzyme, its affinity for other ligand molecules increases.

• $\kappa <0$ (Negatively cooperative binding): Once one ligand molecule is bound to the enzyme, its affinity for other ligand molecules decreases

• $\kappa =0$ (Noncooperative binding): The affinity of the enzyme for a ligand molecule is not dependent on whether or not other ligand molecules are already bound.

CBM-B’s mechanism of reducing C. difficile population is far more complicated than just a single ligand-receptor interaction. Therefore, the κ in this model can be seen as a coefficient reflecting the extent of cooperativity among multiple ligand binding sites, which eventually causes the reduction of C. difficile and inhibits the colony’s growth. This model explored how bacteriocin will definitely increase our CreSolve’s values.

### Bacteriocin dynamics

From the graph below, when ${\mathrm{EC}}_{50}$ value is fixed, the slope becomes steeper as the $\kappa$ value increases (Fig. 10). Similar trend can be found with fixed $\kappa$ value and varied ${\mathrm{EC}}_{50}$ value (Fig. 11). Larger $\kappa$ value and larger ${\mathrm{EC}}_{50}$ value both indicate the increasing interaction strength of bacteriocin on C. difficile.

After fitting with experimental data provided by one of our PIs, Professor Huang, we obtained a $\kappa$ value around 11, which shows that bacteriocin CBM-B reduces C. difficile efficiently. It also proves that CBM-B is a good choice for us when engineering a p-Cresol-reducing E. coli Nissle.

$\kappa$ = 11 is then used to simulate the time and concentration of bacteriocin required to perform effectively (Fig. 12). It shows that when bacteriocin concentration is higher than 1 $\mu g/\mathrm{ml}$, C. difficile survival rate reaches a plateau everytime, which indicates that no matter how high the concentration is, it still takes time for bacteriocin to work. After 6 hours, the C. difficile survival rate nearly reaches 1. Therefore, we speculate that the time bacteriocin CBM-M needs to fully kill C. difficile is 6 hours. Fig. 11. Bacteriocin efficiency at κ= 11 with varying EC50 values. Fig. 12. Bacteriocin efficiency at κ = 11 with varying interaction times.

## References

1. Bomze, I. M. (1983). Lotka-Volterra equation and replicator dynamics: A two-dimensional classification. Biological Cybernetics, 48(3), 201–211. doi: 10.1007/bf00318088
2. Newbold, C. J. (n.d.). Microbial feed additives for ruminants. Biotechnology in Animal Feeds and Animal Feeding, 259–278. doi: 10.1002/9783527615353.ch13
3. Selmer, T., & Andrei, P. I. (2001). p-Hydroxyphenylacetate decarboxylase from Clostridium difficile. European Journal of Biochemistry, 268(5), 1363–1372. doi: 10.1046/j.1432-1327.2001.02001.x
4. Pettit, L. J., Browne, H. P., Yu, L., Smits, W., Fagan, R. P., Barquist, L., … Lawley, T. D. (2014). Functional genomics reveals that Clostridium difficile Spo0A coordinates sporulation, virulence and metabolism. BMC Genomics, 15(1), 160. doi: 10.1186/1471-2164-15-160
5. Berner, M., Krug, D., Bihlmaier, C., Vente, A., Muller, R., & Bechthold, A. (2006). Genes and Enzymes Involved in Caffeic Acid Biosynthesis in the Actinomycete Saccharothrix espanaensis. Journal of Bacteriology, 188(7), 2666–2673. doi: 10.1128/jb.188.7.2666-2673.2006
6. Bennett BD, Kimball EH, Gao M, Osterhout R, Van Dien SJ, Rabinowitz JD. Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. Nat. Chem. Biol. 2009 Aug;5(8):593–9.
7. Liu, D., Wei, Y., Liu, X., Zhou, Y., Jiang, L., Yin, J., … Zhang, Y. (2018). Indoleacetate decarboxylase is a glycyl radical enzyme catalysing the formation of malodorant skatole. Nature Communications, 9(1). doi: 10.1038/s41467-018-06627-x
8. Passmore, I. J., Letertre, M. P. M., Preston, M. D., Bianconi, I., Harrison, M. A., Nasher, F., … Dawson, L. F. (2018). Para-cresol production by Clostridium difficile affects microbial diversity and membrane integrity of Gram-negative bacteria. PLOS Pathogens, 14(9). doi: 10.1371/journal.ppat.1007191 