Introduction

We engineered a metabolic pathway for naringenin production in E. coli Nissle 1917 in order to produce the probiotics for the prevention part of our project.

We saw an opportunity to expedite the engineering process by using mathematical modeling to pick promoters to use in our engineered pathway.

Derivation

Creating a model that would be able to accurately estimate the amount of naringenin produced by the pathway is an infeasible task before doing any practical experiments. However, we are able to write down a simple model with which we could study the speed of the reactions and that would help us decide on the strength of promoters that should be used.

To this end, we use the staple of modeling in synthetic biology: Michaelis–Menten kinetics. That is we model the following type of enzymatic reaction: $$[E] + [S] \leftrightarrows [ES] \rightarrow [E] + [P],$$ with differential equations: $$\frac{d [P]}{dt} = k_{cat}[E]\frac{[P]}{K_m + [P]},$$ $$\frac{d [S]}{dt} = -k_{cat}[E]\frac{[P]}{K_m + [P]},$$ here \([E]\), \([S]\), \([P]\) are the concentrations of enzyme, substrate and product respectively (and \([x]\) is going to denote the concentration of species \(x\) in all that follows), \(k_{cat}\) is a constant called the turnover number and \(K_m\) is a constant that is called Michaelis constant.

In order to model the concentration of mRNA and enzymes, we use the following differential equations: $$\frac{d[mRNA]}{dt} = \alpha_{mRNA} - \beta_{mRNA}[mRNA],$$ $$\frac{d[Enzyme]}{dt} = \alpha_{enzyme}[mRNA] - \beta_{enzyme}[Enzyme],$$ here \(\beta\)’s denote the decay rates, \(\alpha_{mRNA}\) denotes the transcription rate and \(\alpha_{enzyme}\) denotes the translation rate.

Our team measured the strengths of candidate promoters relative to each other. In other words, we measured how many times a specific promoter is stronger or weaker than the promoter that was used as the positive control, as can be seen from Table 1.

We would like this measurement to be reflected in our model. Thus, we denote some base transcription rate (specified later) as \(\zeta\) and write: $$\alpha_{mRNA} = \gamma_{mRNA}\zeta.$$

Now, our goal is to model the pathway depicted in Fig. 1.

The pathway can be expressed by the following chemical reactions: \begin{equation} \emptyset \rightarrow mRNA(TAL) \rightarrow \emptyset, \end{equation} \begin{equation} \emptyset \rightarrow mRNA(4CL) \rightarrow \emptyset, \end{equation} \begin{equation} \emptyset \rightarrow mRNA(CHS) \rightarrow \emptyset, \end{equation} \begin{equation} \emptyset \rightarrow mRNA(CHI) \rightarrow \emptyset, \end{equation} \begin{equation} mRNA(TAL) \rightarrow mRNA(TAL) + TAL, \end{equation} \begin{equation} TAL \rightarrow \emptyset, \end{equation} \begin{equation} mRNA(4CL) \rightarrow mRNA(4CL) + 4CL, \end{equation} \begin{equation} 4CL \rightarrow \emptyset, \end{equation} \begin{equation} mRNA(CHS) \rightarrow mRNA(CHS) + CHS, \end{equation} \begin{equation} CHS \rightarrow \emptyset, \end{equation} \begin{equation} mRNA(CHI) \rightarrow mRNA(CHI) + CHI, \end{equation} \begin{equation} CHI \rightarrow \emptyset, \end{equation} $$TYR + TAL \rightarrow CACID + TAL,$$ $$CACID + 4CL + CoA \rightarrow CCoA + 4CL,$$ $$CCoA + CHS + 3 \times MalCoA \rightarrow NCHAL + CHS + 4 \times CoA,$$ $$NCHAL + CHI \rightarrow NAR + CHI,$$ $$NAR \rightarrow \emptyset.$$

If we assume that there is an infinite (or alternatively very large) amount of tyrosine, CoA and Mal-CoA (if we wished to model the amount of naringenin produced, then assumption that the concentration of Mal-CoA is infinite would be incorrect as this seems to be the major bottleneck of the pathway. However here we only wish to study the reaction speeds, thus we believe that the assumption is valid for this purpose), we can model these reactions by the following system of differential equations: \begin{equation} \frac{d(TAL)}{dt} = \gamma_{TAL}\zeta - \beta_{m(TAL)}(TAL), \end{equation} \begin{equation} \frac{d(4CL)}{dt} = \gamma_{4CL}\zeta - \beta_{m(4CL)}(4CL), \end{equation} \begin{equation} \frac{d(CHS)}{dt} = \gamma_{CHS}\zeta - \beta_{m(CHS)}(CHS), \end{equation} \begin{equation} \frac{d(CHI)}{dt} = \gamma_{CHI}\zeta - \beta_{m(CHI)}(CHI), \end{equation} \begin{equation} \frac{d[TAL]}{dt} = \alpha_{TAL}(TAL) - \beta_{TAL}[TAL], \end{equation} \begin{equation} \frac{d[4CL]}{dt} = \alpha_{4CL}(4CL) - \beta_{4CL}[4CL], \end{equation} \begin{equation} \frac{d[CHS]}{dt} = \alpha_{CHS}(CHS) - \beta_{CHS}[CHS], \end{equation} \begin{equation} \frac{d[CHI]}{dt} = \alpha_{CHI}(CHI) - \beta_{CHI}[CHI], \end{equation} $$\frac{d[CACID]}{dt} = k_{TAL}[TAL] - k_{4CL}[4CL]\frac{[CACID]}{K_{4CL} + [CACID]},$$ $$\frac{d[CCoA]}{dt} = k_{4CL}[4CL]\frac{[CACID]}{K_{4CL} + [CACID]} - k_{CHS}[CHS]\frac{[CCoA]}{K_{CHS} + [CCoA]},$$ $$\frac{d[NCHAL]}{dt} = k_{CHS}[CHS]\frac{[CCoA]}{K_{CHS} + [CCoA]} - k_{CHI}[CHI]\frac{[NCHAL]}{K_{CHI} + [NCHAL]},$$ $$\frac{d[NCHAL]}{dt} = k_{CHI}[CHI]\frac{[NCHAL]}{K_{CHI} + [NCHAL]} - \beta_{NAR}[NAR],$$ here \((x)\) denotes \([mRNA(x)]\), small \(k\)’s denote the appropriate turnover numbers and big \(K\)’s denote the appropriate Michaelis constants.

This model is overly complicated for our purposes. We can reduce it by noting that the reactions \((1) - (12)\) happen on a faster time scale then the rest. Therefore, we can assume that the reactions \((1) - (12)\) are in the steady state for the entirety of the process. With this assumption we have additional conditions: \begin{equation} \frac{d(TAL)}{dt} = 0, \end{equation} \begin{equation} \frac{d(4CL)}{dt} = 0, \end{equation} \begin{equation} \frac{d(CHS)}{dt} = 0, \end{equation} \begin{equation} \frac{d(CHI)}{dt} = 0, \end{equation} \begin{equation} \frac{d[TAL]}{dt} = 0, \end{equation} \begin{equation} \frac{d[4CL]}{dt} = 0, \end{equation} \begin{equation} \frac{d[CHS]}{dt} = 0, \end{equation} \begin{equation} \frac{d[CHI]}{dt} = 0. \end{equation}

By combining \((13)-(16)\) with \((21)-(24)\) we get $$(x) = \frac{\gamma\zeta}{\beta_{mRNA}},$$ and then by combining \((17)-(20)\) with \((25)-(28)\) we get $$[x] = \frac{\alpha\gamma\zeta}{\beta_{mRNA}\beta_{enzyme}}.$$

We can additionally assume that translation rates and decay rates of mRNA and enzyme are similar for different species. Then by taking the base transcription rate \(\zeta\) such that $$\frac{\alpha\zeta}{\beta_{mRNA}\beta_{enzyme}}$$ is equal to 1 we can reduce the original model to a simpler model with less equations: $$\frac{d[CACID]}{dt} = k_{TAL}\gamma_{TAL} - k_{4CL}\gamma_{4CL}\frac{[CACID]}{K_{4CL} + [CACID]},$$ $$\frac{d[CCoA]}{dt} = k_{4CL}\gamma_{4CL}\frac{[CACID]}{K_{4CL} + [CACID]} - k_{CHS}\gamma_{CHS}\frac{[CCoA]}{K_{CHS} + [CCoA]},$$ $$\frac{d[NCHAL]}{dt} = k_{CHS}\gamma_{CHS}\frac{[CCoA]}{K_{CHS} + [CCoA]} - k_{CHI}\gamma_{CHI}\frac{[NCHAL]}{K_{CHI} + [NCHAL]},$$ $$\frac{d[NAR]}{dt} = k_{CHI}\gamma_{CHI}\frac{[NCHAL]}{K_{CHI} + [NCHAL]} - \beta_{NAR}[NAR].$$

Analysis

We see that in the steady state we have $$[NAR] = \frac{k_{TAL}\gamma_{TAL}}{\beta_{NAR}}.$$ This makes intuitive sense - the more substrate one puts in, the more product one expects to get. However, the steady-state might take an exorbitant amount of time to reach depending on the parameters. Thus, we decided to study the system after simulating it for 16 hours (taking the initial concentrations of all proteins in the pathway to be 0) as these are the timescales that the performance of the engineered pathway would be measured in.

Next, we researched the literature to compile probable values for turnover numbers and Michaelis constants. We came up with the following figures:

From Table 2 we see that the reaction producing naringenin chalcone seems to be around 10 times slower than the second slowest one in the pathway. This makes sense since this is a sequential reaction involving 4 molecules. Seeing this, we hypothesized that this reaction is the major bottleneck of the pathway. That is, the only parameters that have a major impact on the output of the model are \(k_{CHS}\) and \(\gamma_{CHS}\).

We validated this hypothesis by performing a simple sensitivity analysis as follows:

1. Generate 10000 samples of parameter values by uniformly sampling from the intervals detailed in Table 4. The average value for \(\beta_{NAR}\) was derived from [10].
2. Simulate the model with generated random parameters for 16 hours and save the concentration of naringenin.
3. Compute the correlation coefficients between the parameters and concentration of naringenin.

The results of sensitivity analysis are presented in Table 5.

The sensitivity analysis confirmed our hypothesis. We note that it also showed that another important parameter is the decay rate of naringenin.

Conclusion

We derived a simple mathematical model for the naringenin pathway that our team wanted to implement in vivo. By performing sensitivity analysis, we determined that the reaction which turns Coumaryl-CoA to naringenin chalcone synthase is the bottleneck of the naringenin synthesis process.

In the end we decided to use the pSlpA promoter to produce the enzyme responsible for the bottleneck reaction as this promoter was the strongest one from the one's that we measured (as can be seen from Table 1). For the rest of the enzymes we could use a weaker promoter according to the model. Thus, we chose BBa_J23101 sequence, since it demonstrated adequate performance in our measurements.