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 get 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 show an opt-knock flux-balance analysis to show how the second culprit, the cellular FPP pool can be overcome.

Basic System

The basic reactions of the Nootkatone pathway that are introduced by our team 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:

Bioreactor simulation

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.

Our model predicts a Nootkatone Yield of $268.14 \ \frac{\text{mg}}{\text{{L}}$ and a Nootkatol Yield of $54 \ \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 model is already quite accurate and that our assumption of a five-fold reduction in the backwards reaction of alcohol dehydrogenase is valid as well, as our $\frac{\text{Nootkatone}}{\text{Nootkatol}}$ ratio was 4.96, which is close to the 4.72 by Wriessnegger. But our model still lacks an important factor, 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 therefore used hill kinetics to fit this behaviour.

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

We fitted the data by Gavira using scipy.optimize.minimize and got the following results. Note that due to numerical reasons the factor n had to be calculated as a float but then was rounded to an integer.

Although Nootkatol is more toxic to the cell than Nootkatone the later is concentrated five times higher in our model. We therefore limited us to the Nootkatol penalty to keep the model simple. We used the arithmetic mean of the two functions as a penalty function for 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.

$$\frac{K_{\text{Tox Nkn}}^2}{K_{\text{Tox Nkn}}^2 + [\text{Nkn}]^2}$$

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 is closest to the actual biological situation. 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 $4.18E-9 \frac{1}{\text{s}}$, which we then also used in all other models.

Penalty model expression

We wanted to check whether over- or underexpression of any of the enzymes could improve our Nootkatone yield. We therefore set up one analysis in which we held two of the three enzyme concentrations constant while varying the third one and a second one in which all enzyme concentrations were varied at the same time.

Since we fitted the FPP concentration an enzyme concentration of $1 \ \text{µM}$, a further increase in enzyme concentration did not result in a higher yield, as was to be expected. But interestingly down-regulating valencen synthase led to an increase in Nootkatone yield, which was probably due to a submaximal penalty. We further wanted to explore how the system reacts when the FPP import is increased by a factor of 10.

Penalty model increased FPP import dynamics

While the maximal yield is increased the penalty function of the model almost approaches $10E-2$, which means that at the end of the production almost all cells are dead. The increased yield will therefore most probably result as an artifact of a very high valencen production in the beginning. Since only the imoprt function is penalized, all further reactions can still take place.

Penalty model increased FPP import expression

The expression analysis further solidifies the interpretation that the overall yield is linearliy dependent on valencene synthase and FPP import and that the increased yield is an artifact of high early valencene production.

Penalty model increased FPP import and enzyme concentration dynamics

Peroxisome model expression

Peroxisome increased FPP import dynamics

Peroxisome model increased FPP import expression

Peroxisome model increased FPP increased enzyme concentration dynamics