Difference between revisions of "Team:Vilnius-Lithuania/Model"

Revision as of 15:50, 21 October 2021

Motivation

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 performed mathematical analysis in order to expedite the engineering process by deciding what promoters to use based on the model for the pathway we derived.

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.

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.

Fig. 1. Naringenin synthesis pathway.

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:

Table 1. Turnover numbers ($k_{cat}$)

Enzyme	Values (1/s)	Average (1/s)	Reference
Tyrosine ammonia-lyase (TAL)	107	119	[1]
	114		[1]
	115		[1]
	139		[1]
4-coumarate-CoA ligase (4CL)	0.2163	0.3354	[2]
	0.2205		[2]
	0.7821		[2]
	0.1225		[2]
Chalcone synthase (CHS)	0.045	0.0575	[3]
	0.178		[4]
	0.11		[4]
	0.085		[4]
	0.05		[4]
	0.0202		[5]
	0.0167		[6]
	0.042		[7]
	0.007		[7]
	0.021		[8]
Chalcone isomerase (CHI)	5	89.5	[9]
	7.8		[9]
	9.6		[9]
	35.2		[9]
	56.9		[9]
	130.3		[9]
	134.7		[9]
	197.7		[9]
	228.2		[9]

Table 2. Michaelis constants ($K_{M}$)

Enzyme	Values (mM)	Average (mM)	Reference
4-coumarate-CoA ligase (4CL)	0.389	0.276	[2]
	0.155		[2]
	0.283		[2]
Chalcone synthase (CHS)	0.0049	0.0049	[7]
Chalcone isomerase (CHI)	0.0024	0.007	[9]
	0.0048		[9]
	0.0048		[9]
	0.0061		[9]
	0.007		[9]
	0.0085		[9]
	0.0086		[9]
	0.0099		[9]
	0.0105		[9]

From Table 1 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:

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

Table 3. Parameter values used in sensitivity analysis

Parameter	Value range
$\gamma_{TAL}$	$0.33 - 3$
$\gamma_{4CL}$	$0.33 - 3$
$\gamma_{CHS}$	$0.33 - 3$
$\gamma_{CHI}$	$0.33 - 3$
$\beta_{NAR}$	$3.6\mathrm{e}{-5} \pm 3.6\mathrm{e}{-6} \: (1/s)$
$k_{TAL}$	$119 \pm 11.9 \: (1/s)$
$k_{4CL}$	$0.3354 \pm 0.034 \: (1/s)$
$k_{CHS}$	$0.0575 \pm 0.006 \: (1/s)$
$k_{CHI}$	$89.5 \pm 8.95 \: (1/s)$
$K_{4CL}$	$0.276 \pm 0.028 \: (mM)$
$K_{CHS}$	$0.0049 \pm 0.0005 \: (mM)$
$K_{CHI}$	$0.007 \pm 0.0007 \: (mM)$

The results of sensitivity analysis are presented in Table 4.

Table 4. Results of sensitivity analysis

Parameter	Correlation coefficient
$\gamma_{TAL}$	$0.0242$
$\gamma_{4CL}$	$0.0339$
$\gamma_{CHS}$	$0.9833$
$\gamma_{CHI}$	$0.0008$
$\beta_{NAR}$	$-0.0938$
$k_{TAL}$	$-0.0113$
$k_{4CL}$	$0.0041$
$k_{CHS}$	$0.1042$
$k_{CHI}$	$-0.0009$
$K_{4CL}$	$-0.0161$
$K_{CHS}$	$-0.0199$
$K_{CHI}$	$-0.0096$

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 is the bottleneck of the naringenin synthesis pathway. Thus we decided to use a pSlpA that we noticed by analysing possible promotors for E. coli and L. paracasei. Furthermore, BBa_J23101, also known as one from the Anderson's collection promotor, was compared and picked up for other naringenin pathway protein expression.

References

1.

Zhou, S., Liu, P., Chen, J., Du, G., Li, H., Zhou, J. (2016). Characterization of mutants of a tyrosine ammonia-lyase from Rhodotorula glutinis. Appl. Microbiol. Biotechnol. 100, 10443-10452. To the article.

2.

Gao, S., Yu, H. N., Xu, R. X., Cheng, A. X., & Lou, H. X. (2015). Cloning and functional characterization of a 4-coumarate CoA ligase from liverwort Plagiochasma appendiculatum. Phytochemistry, 111, 48–58. To the article.

3.

Guo, H.-L., Yang, Y.-D., Ma, Y.-D., Liu, W.-B., Feng, J., Luo, Z.-Q., … Ma, L.-Q. (2016). A bifunctional type III polyketide synthase from raspberry (Rubus idaeus L.) with both chalcone synthase and benzalacetone synthase activity. Journal of Plant Biochemistry and Biotechnology, 26(1), 80–90. To the article.

4.

Shen, Y., Li, X., Chai, T., & Wang, H. (2016). Outer-sphere residues influence the catalytic activity of a chalcone synthase from Polygonum cuspidatum. FEBS open bio, 6(6), 610–618. To the article.

5.

Stewart, C., Jr, Woods, K., Macias, G., Allan, A. C., Hellens, R. P., & Noel, J. P. (2017). Molecular architectures of benzoic acid-specific type III polyketide synthases. Acta crystallographica. Section D, Structural biology, 73(Pt 12), 1007–1019. To the article.

6.

Abe, I., Watanabe, T., & Noguchi, H. (2004). Enzymatic formation of long-chain polyketide pyrones by plant type III polyketide synthases. Phytochemistry, 65(17), 2447–2453. To the article.

7.

Liu, B., Falkenstein-Paul, H., Schmidt, W., & Beerhues, L. (2003). Benzophenone synthase and chalcone synthase from Hypericum androsaemum cell cultures: cDNA cloning, functional expression, and site-directed mutagenesis of two polyketide synthases. The Plant journal : for cell and molecular biology, 34(6), 847–855. To the article.

8.

Morita, H., Takahashi, Y., Noguchi, H., & Abe, I. (2000). Enzymatic formation of unnatural aromatic polyketides by chalcone synthase. Biochemical and biophysical research communications, 279(1), 190–195. To the article.

9.

Park, S. H., Lee, C. W., Cho, S. M., Lee, H., Park, H., Lee, J., & Lee, J. H. (2018). Crystal structure and enzymatic properties of chalcone isomerase from the Antarctic vascular plant Deschampsia antarctica Desv. PloS one, 13(2), e0192415. To the article.

10.

Kanaze, F. I., Bounartzi, M. I., Georgarakis, M., & Niopas, I. (2006). Pharmacokinetics of the citrus flavanone aglycones hesperetin and naringenin after single oral administration in human subjects. European Journal of Clinical Nutrition, 61(4), 472–477. To the article.

@@ Line 65: / Line 65: @@
                    <p> 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. </p>
-                   <p> 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
+                   <p> 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.</p>
-                      seen from <b>Table 1.</b></p>
                    <p> 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.$$ </p>
                    <p> Now, our goal is to model the pathway depicted in <b>Fig. 1.</b></p>
@@ Line 376: / Line 375: @@
                    We validated this hypothesis by performing a simple sensitivity analysis as follows: </p>
                    <ol>
-                       <li> Generate 10000 samples of parameter values by uniformly sampling from the intervals detailed in <b>Table 4.</b> The average value for \(\beta_{NAR}\) was derived from [10]. </li>
+                       <li> Generate 10000 samples of parameter values by uniformly sampling from the intervals detailed in <b>Table 3.</b> The average value for \(\beta_{NAR}\) was derived from [10]. </li>
                        <li>Simulate the model with generated random parameters for 16 hours and save the concentration of naringenin. </li>
                        <li>Compute the correlation coefficients between the parameters and concentration of naringenin. </li>
@@ Line 441: / Line 440: @@
                        </table>
                    </div>
-                   <p> The results of sensitivity analysis are presented in <b>Table 5.</b></p>
+                   <p> The results of sensitivity analysis are presented in <b>Table 4.</b></p>
                    <div class="table-container">
                        <div class="table-headline"> <b>Table 4.</b> Results of sensitivity analysis </div>

Parameter	Value range
\(\gamma_{TAL}\)	\(0.33 - 3\)
\(\gamma_{4CL}\)	\(0.33 - 3\)
\(\gamma_{CHS}\)	\(0.33 - 3\)
\(\gamma_{CHI}\)	\(0.33 - 3\)
\(\beta_{NAR}\)	\(3.6\mathrm{e}{-5} \pm 3.6\mathrm{e}{-6} \: (1/s)\)
\(k_{TAL}\)	\(119 \pm 11.9 \: (1/s)\)
\(k_{4CL}\)	\(0.3354 \pm 0.034 \: (1/s)\)
\(k_{CHS}\)	\(0.0575 \pm 0.006 \: (1/s)\)
\(k_{CHI}\)	\(89.5 \pm 8.95 \: (1/s)\)
\(K_{4CL}\)	\(0.276 \pm 0.028 \: (mM)\)
\(K_{CHS}\)	\(0.0049 \pm 0.0005 \: (mM)\)
\(K_{CHI}\)	\(0.007 \pm 0.0007 \: (mM)\)

Parameter	Correlation coefficient
\(\gamma_{TAL}\)	\(0.0242\)
\(\gamma_{4CL}\)	\(0.0339\)
\(\gamma_{CHS}\)	\(0.9833\)
\(\gamma_{CHI}\)	\(0.0008\)
\(\beta_{NAR}\)	\(-0.0938\)
\(k_{TAL}\)	\(-0.0113\)
\(k_{4CL}\)	\(0.0041\)
\(k_{CHS}\)	\(0.1042\)
\(k_{CHI}\)	\(-0.0009\)
\(K_{4CL}\)	\(-0.0161\)
\(K_{CHS}\)	\(-0.0199\)
\(K_{CHI}\)	\(-0.0096\)