Model

Pex-5 Structure Model

Nootkatone Metabolic Model

Overview

In the following we present our model of the Nootkatone biosynthesis pathway, to give you an insight into its behavior and dynamics. We start with an oversimplified luminar model to give you a sense for the behavior of the enzymes in the pathway. Then we continue with a model introducing a function penalizing high concentrations of toxic products. As the toxicity is one of two main culprits of Nootkatone production, we further show the production inside a peroxisome, as we assume that intermediates inside the peroxisome cannot pass the membrane and thus have no toxic effect on the cells. Then we discuss possible ways to overcome the shortage of the cellular FPP pool.

Basic System

The basic reactions of the Nootkatone pathway are the following.

However during research we found that using the p450-BM3 enzyme will simplify and enhance Nootkatone production, giving the following reaction pathway.

We assumed Michaelis-Menten kinetics for each reaction, with the last step being reversible.

Michaelis-Menten kinetics

$$\frac{dP}{dt} = \frac{V_{Max} \cdot c_{S}}{K_{M} + c_{S}}$$

Reversible Michaelist-Menten kinetics

$$\frac{dP}{dt} = \frac{\frac{V_{M+} \cdot c_{S}}{K_{M+}} - \frac{V_{M-} \cdot c_{P}}{K_{M-}}}{1 + \frac{c_{S}}{K_{M+}} + \frac{c_{P}}{K_{M-}}}$$

We further assumed a permanent FPP production, which we fitted in the Penalty model.

This gives us the following system of differential equations.

$$\frac{dFPP}{dt} = v_{\text{FPP in}} - \frac{V_{Max,ValS} \cdot c_{FPP}}{K_{M, ValS} + c_{FPP}} $$

$$\frac{dValencene}{dt} = \frac{V_{Max,ValS} \cdot c_{FPP}}{K_{M, ValS} + c_{FPP}} -\frac{V_{Max,p450\_BM3} \cdot c_{Valencene}}{K_{M, p450\_BM3} + c_{Valencene}}$$

$$\frac{dNootkatol}{dt} = \frac{V_{Max,p450\_BM3} \cdot c_{Valencene}}{K_{M, p450\_BM3} + c_{Valencene}} - \frac{\frac{V_{M,ADH+} \cdot c_{Nootkatol}}{K_{M,ADH+}} - \frac{V_{M,ADH-} \cdot c_{Nootkatone}}{K_{M,ADH-}}}{1 + \frac{c_{Nootkatol}}{K_{M,ADH+}} + \frac{c_{Nootkatone}}{K_{M,ADH-}}}$$

$$\frac{dNootkatone}{dt} = \frac{\frac{V_{M,ADH+} \cdot c_{Nootkatol}}{K_{M,ADH+}} - \frac{V_{M,ADH-} \cdot c_{Nootkatone}}{K_{M,ADH-}}}{1 + \frac{c_{Nootkatol}}{K_{M,ADH+}} + \frac{c_{Nootkatone}}{K_{M,ADH-}}}$$

Parameters

Our crude first assumption was a constant enzyme concentration of $1 \ \text{µM}$. We mostly used kinetic parameters from the Brenda database. If S. cerevisiae was not available as the host we used the nearest relative, if the reaction was not available we used the kinetics for the respective co-substrate as an upper boundary. Another assumption we made is a five-fold reduction in the speed of the reversible reaction of the ADH-21, based on the work by Sebastian Schulz (2015), that the forward reaction is favored.

Simple model

In order to simulate the system of differential equations we set up a few python files, one containing a model class which consists of the parameters, the models and the integration function, one containing all the kinetics used and one which we used to actually run the simulations. We wanted to show some example code on how to set up such a project, while the code size was still managable, in order to help new iGEM teams that might need help with their modeling.

Kinetics File

		  
		def ConstantInflux(substrate,speed):
		    return(speed)

		def MichaelisMenten(
		        substrate,enzymeConcentration,kcat,km):
		    return(
		            (enzymeConcentration*kcat*substrate)
		            /(km+substrate))

		def MichaelisMentenReversible(
		        substrate,product,enzymeConcentration,
		        kcats,kcatp,kms,kmp):
		    return (
		            (
		                    ((enzymeConcentration*kcats*substrate)/kms)
		                    -((enzymeConcentration*kcatp*product)/kmp)
		                    )
		                    /(1+(substrate/kms)+(product/kmp)
		                    ))
		  
		

Model File

		  
		import kinetics as k
		import numpy as np
		import scipy.integrate

		class model(object):
		    def __init__(self):
		        self.c_ValS = 1E-6 #[M]
		        self.c_ADH = 1E-6 #[M]
		        self.c_p450 = 1E-6 #[M]

		        self.km_ValS = 1.04E-6  #[M]
		        self.km_p450 = 126E-6 #[M]
		        self.km_ADH = 161E-6 #[M]

		        self.kcat_ValS = 0.0032*3600 #[1/s]
		        self.kcat_p450 = 6*3600 #[1/s]
		        self.kcat_ADH = 2.619*3600 #[1/s]
		        self.k_influxSpeed = 1.03562761862e-05 #[1/s]


		    def FPPInflux(self,FPP):
		        return(
		                k.ConstantInflux(
		                        FPP,
		                        self.k_influxSpeed)
		                )

		    def ValenceneSynthase(self,FPP):
		        return(
		                k.MichaelisMenten(
		                        FPP,
		                        self.c_ValS,
		                        self.kcat_ValS,
		                        self.km_ValS)
		                )
		    def p450(self,Valencene):
		        return(
		                k.MichaelisMenten(
		                        Valencene,
		                        self.c_p450,
		                        self.kcat_p450,
		                        self.km_p450)
		                )
		    def ADH(self,Nootkatol,Nootkatone):
		        return (
		                k.MichaelisMentenReversible(
		                        Nootkatol,
		                        Nootkatone,
		                        self.c_ADH,
		                        self.kcat_ADH,
		                        self.kcat_ADHrev,
		                        self.km_ADH,
		                        self.km_ADH)
		                )

		    def p450_Model(self,time,initvalues):
		        FPP,Val,Nkl,Nkn = initvalues
		        dFPPdt = self.FPPInflux(FPP) - self.ValenceneSynthase(FPP)
		        dValdt = self.ValenceneSynthase(FPP) - self.p450(Val)
		        dNkldt = self.p450(Val) - self.ADH(Nkl,Nkn)
		        dNtkdt = self.ADH(Nkl,Nkn)
		        return [dFPPdt,dValdt,dNkldt,dNtkdt]


		    def Integration(self,time,initvalues,model):
		        Y = np.zeros([len(time),len(initvalues)])
		        Y[0] = np.transpose(initvalues)
		        r = scipy.integrate.ode(model)
		        r.set_integrator('lsoda')
		        r.set_initial_value(initvalues)
		        cnt=1
		        while r.successful() and cnt < len(time):
		            Y[cnt,:] = r.integrate(time[cnt])
		            cnt+=1
		        return Y
		  
		

Simulation file

		  
		from model import model
		import numpy as np
		import matplotlib.pyplot as plt

		m = model()
		time = np.linspace(0,1,1000)
		data = m.Integration(time,[0,0,0,0],m.p450_Model)
		plt.plot(time,data[:,0],label=(
		        r"FPP, $k_{FPP}$ = "
		        +"{:.2e}".format(m.k_influxSpeed)
		        +r" $\frac{1}{h}$"
		        ))
		plt.plot(time,data[:,1],label=(
		        r"Valencene, $c_{ValS}$ = "
		        +"{:.0e}".format(m.c_ValS)
		        +r" $\frac{mol}{L}$"
		        ))
		plt.plot(time,data[:,2],label=(
		        r"Nootkatol, $c_{p450}$ = "
		        +"{:.0e}".format(m.c_p450)
		        +r" $\frac{mol}{L}$"
		        ))
		plt.plot(time,data[:,3],label=(
		        r"Nootkatone, $c_{ADH-21}$ = "
		        +"{:.0e}".format(m.c_ADH)
		        +r" $\frac{mol}{L}$"
		        ))
		plt.yscale("log")
		plt.legend(loc=9, bbox_to_anchor=(0.5, -0.2), ncol=2)
		plt.title("First Simulation")
		plt.xlabel(r"Time in $[h]$")
		plt.ylabel(r"Yield in $\frac{mol}{cell}$")
		plt.show()
		  
		

The results of our first simulation are the following:

What is visible in this first simulation is that the FPP import is faster than Valencene Synthase since FPP is agglomerating. However Valencene synthase, p450 and ADH are all saturated around 20 minutes, since the valencene concentration becomes stationary and Nootkatol and Nootkatone levels keep rising proportionally, due to the back-reaction of ADH.

Bioreactor model

In order to check the validity of our model we took the results of Tamara Wriessnegger 2014 as a point for reference. Wriessnegger achieved $208 \ \frac{\text{mg}}{\text{L}}$ Nootkatone production after 108 h in a in-vitro production in Pichia pastoris. We therefore had to scale up our single cell model to a bioreactor simulation.

We first calculated the number of product molecules per cell.

$$\frac{\text{N}_{\text{Product}}}{\text{cell}} \ [\text{mol}]= c_{Product} [\frac{\text{mol}}{\text{L}}]\cdot V_{\text{Yeast}} \ [\text{L}]$$

Then we calculated the number of cells per Liter bioreactor.

$$ \text{cells per Liter} \ [\frac{\text{1}}{\text{L}}] = \frac{\rho_{\text{Cells in Bioreactor}} \ [\frac{\text{g d.w.}}{\text{L}}] }{m_{\text{Yeast dry weight}} \ [\text{g d.w.}] } $$

With these two equations we could scale up the production of a single cell for each product.

$$\text{Yield}_{\text{Bioreactor}} \ [\frac{\text{mg}}{\text{L}}] = \frac{\text{N}_{\text{Product}}}{\text{cell}} \ [\text{mol}] \cdot \text{cells per Liter} \ [\frac{\text{1}}{\text{L}}] \cdot \text{M}_{\text{Product}} \ [\frac{\text{g}}{\text{mol}}] \cdot 1000$$

We used the following constants for the calculations above.

Bioreactor simulation

Our model predicts a Nootkatone Yield of $154.24 \ \frac{\text{mg}}{\text{L}}$ and a Nootkatol Yield of $31.06 \ \frac{\text{mg}}{\text{L}}$. Wriessnegger produced $208 \ \frac{\text{mg}}{\text{L}}$ Nootkatone and $44 \ \frac{\text{mg}}{\text{L}}$ Nootkatol. This showed us that our assumption of a five-fold reduction in the backwards reaction of alcohol dehydrogenase is valid, as our $\frac{\text{Nootkatone}}{\text{Nootkatol}}$ ratio was 4.96, which is quite close to the 4.72 by Wriessnegger. We then took a look at how differently expressing the enzymes would have an influence on Nootkatone production.

Expression analysis

The first figure shows the Nootkatone yield based on a varied enzyme concentration of all enzymes, while the other enzyme concentrations were kept constant at $1 \ µM$. The second plot shows varied enzyme concentrations for all enzymes at the same time. The smallest concentration with almost maximal Nootkatone production is at about $2 \ µM$ Valencene Synthase, when not over increasing the other enzymes. Therefore we chose to double the expression of Valencene Synthase to ensure that our pathway is saturated.

This increased our yield above the results by Wriessnegger. But our model still lacks the important factor of product toxicity at high concentration, which is why we did not fit the FPP import to this, but to the next model.

Toxicity

According to Carole Gavira 2013 both Nootkatone and Nootkatol exhibit toxicity at high levels. At $100 \ \frac{\text{mg}}{\text{L}}$ Nootkatol the cell viability was decreased to 55%, at $100 \ \frac{\text{mg}}{\text{L}}$ Nootkatone to 80%. For both substances all cells were dead around $1000 \ \frac{\text{mg}}{\text{L}}$. We used hill kinetics to fit this behaviour.

$$\frac{K_{Tox}^n}{K_{Tox}^n + [S]^n}$$

We fitted the data of Gavira using scipy.optimize.brent and multiple factors n. This is due to the fact that n is only valid in integer values and scipy does not support mixed-integer optimizations. We further chose to only implement Nootkatone toxicity, even if it is the less toxic compound, as it is concentrated five times higher than Nootkatol in the cell. This also helped to keep our model simple.

Even though the error produced by functions with a higher n was decreasing we chose to use n = 3 to prevent overfitting, since n = 3 already fit both points quite well.

We used this penalty function only on our FPP import. The idea behind that was that since the whole pathway depends linearly on the FPP influx it can act representative for the cell viability. This turned out to show some problems later in the model but we chose to keep this assumption, as the problems can be addressed by proper interpretation of the data and because we wanted to keep our model simple.

This changes our system in the following way.

$$\frac{dFPP}{dt} = v_{\text{FPP in}} \cdot \frac{K_{\text{Tox Nkn}}^2}{K_{\text{Tox Nkn}}^2 + [\text{Nkn}]^2} - \frac{V_{Max,ValS} \cdot c_{FPP}}{K_{M, ValS} + c_{FPP}} $$

Penalty model

Dynamics

This is the model where we actually fitted our import parameter to the data supplied by Wriessnegger, as we assumed it to be sufficiently close to the actual biological reality. We did this by using scipy.optimize.brent to minimize the error function $\text{abs}(208 - [S]_{\text{Nootkatol}})$. This gave us an import rate of $6.55E-9 \frac{1}{\text{s}}$, which we then also used in all other models, including the ones shown previously.

The Nootkatone yield is approximately at the results of Wriessnegger and we can see that the FPP import declines with increasing Nootkatone production. The penalty function can be interpreted as the fraction of alive cells, which at the final point of the simulation is $52.16 \ \%$. This huge number already gives us an understanding that Nootkatol and Nootkatone production limit themselves in cytosolic production, thus showing us the first culprit of Nootkatone production.

Expression

We wanted to check whether over- or underexpression of any of the enzymes could improve our Nootkatone yield under penalty conditions.

Interestingly the peak Nootkatone yield is now under submaximal Valencene Synthase expression. This is possibly due to the increased toxicity. In order to check this we further increased the FPP import tenfold over the fitted value.

Increased FPP import dynamics

The maximal yield under increased FPP import conditions is increased, but at the end of the simulation only around $11 \ \%$ of the cell population is alive. This is the problem of the model that we talked earlier about. In this case the early FPP import so fast compared to the enzymatic reactions so that FPP accumulates. Since only the import is penalized, but not the other reactions, all that follows is an simulation of dead cells still transforming an large FPP pool. This is obviously biologically wrong and an artifact of the way the penalty is set up. The effects of this error can be seen in the variation of single enzyme concentrations in the plots below.

Increased FPP import expression

Increased FPP import and enzyme concentration dynamics

Even with this submaximal production only $4 \ \%$ of cells are alive at the end of the simulation. So even if the results of the production with a tenfold increased FPP influx have to be interpreted carefully, the general takeaway of the simulation is that given a sufficiently high FPP influx, Nootkatone toxicity is the main remaining culprit of Nootkatone production. This is where our Nootkatone subproject comes in. By transferring the enzymes needed for Nootkatone production into the peroxisome we assume that the products cannot diffuse back into the cell and therefore do not impose any toxicity on the cells. This can be seen by the great increase in Nootkatone production in the models below.

Peroxisome model

Expression

Peroxisome increased FPP import dynamics

Peroxisome model increased FPP import expression

Peroxisome model increased FPP increased enzyme concentration dynamics

Conclusion

As can be seen from our simulations, the two main culprits of Nootkatone production are the toxicity of the products at high concentrations and the FPP influx. The toxicity can be avoided by transfering the pathway into the peroxisome, which is one of the projects of our team. To address the problem of FPP influx our team did an opt-knock FBA analyis to show how S. cerevisiae cells can by engineered to stabily produce more FPP. Note that the values marked with an asterix are assumed to be false due to the model setup as explained in the sections above.

Further thoughts

A false assumption in the peroxisomal models above is that the size of the peroxisome is equal to that of the yeast cell. This is obviously not correct. We wanted to get a feeling how large the peroxisome actually has to be in order for production to be superior to cytosolic production. We therefore took the the results of the cytosolic model with increased FPP intakte and enzyme concentration and then calculated the yield of the corresponding peroxisomal model with varying peroxisome size.

We calculated the volume for equal production

$$V_{\text{Perox. Equal}} = \frac{\text{Yield}_{\text{Cyt.}} \cdot V_{\text{Yeast}}}{\text{Yield}_{\text{Perox.}}}$$

to be $1.24E-14 \ \text{L}$, so approximately a fifth of the cell volume. To further increase the possible yield our team is working on a Pex11-knockout that aims to increase peroxisome size

Optknock Analysis