Approach

---

Modeling is an inherent part of the engineering cycle, even more so in the context of biomanufacturing, a discipline that relies on optimization when it comes to assessing feasibility. For our group, modeling was the driving force behind our work. Its importance was emphasized by the feedback we received during our meetings with stakeholders, which inquired regarding the production rates and costs. Thus, helping us identify which aspects needed to be modeled.

We worked on multiple models and employed a stacking approach, where the output of one model was used as the input of another. This was done to optimize the system for large-scale production of our target metabolite - decursin. Our biomanufacturing process will be the synthesis of decursinol from umbelliferone. Later on, decursinol can be converted chemically and cheaply to decursin.

Our models can be divided into four categories:

1. Diffusion of umbelliferone (UMB), a precursor of decursin, into E. coli BL21.
2. The fermentation process and production of the enzymes for the anabolic pathway: Prenyltransferase (PT), XimD, and XimE.
3. Computational approach to site-directed mutagenesis of the XimE protein to improve its affinity to a specific substrate.
4. The enzymatic reactions and kinetics of PT, XimE, and XimD to convert UMB to DMS and DMS to decursinol (DEC).

Before we set out to design the biomanufacturing process itself, our project started with a single basic assumption: the precursor UMB can diffuse into the cell, where the enzymatic reaction takes place. To examine this assumption, we wanted to model the diffusion process (1). To this end, we turned to the work done by iGEM CLS_CLSG_UK 2020 and expanded upon their efforts. We made the needed adjustments according to our needs and used the prediction of almost-immediate diffusion as output. We then wanted to establish a way to predict production rates, by focusing on three factors: bacterial population size, enzyme expression rates, and their kinetics. This led us to the development of a fermentation model (2), whose prediction of total protein expression was used as an output. Regarding the kinetics of the enzymes, previous work has shown that XimE and XimD are not substrate-specific^[1]. To address this, and improve XimE’s catalytic efficiency in the decursin pathway, we employed an approach of using computational tools to model the substrate binding affinity of different mutations (3). The assumption of diffusion-based kinetics (1) and the protein concentration from (2) were used as input for an enzyme kinetics model (4). This model took into account the co-dependence of the enzymes in the pathway using statistical tools and was used to predict production rates. Furthermore, the insights regarding the limiting factor and byproducts resulted in a redesign of our plasmid system and were used to assess its feasibility. Afterward, the final model was fitted to the experimental data from the wet lab, to improve its predictions. All of the graphs displayed here are generated using our MATLAB code, which is open source and available here: Gitlab link.

1. Diffusion of Umbilifferon into E. coli BL21 Model

Introduction

---

In the anabolic pathway of DEC, the first metabolite is UMB which should be supplied externally, and be able to diffuse the membrane of E. coli cells to start the process. Therefore, we wanted to show the occurrence of diffusion using known equations and mathematical derivatives.

Definitions

---

Table 1: Description of Model Symbols and Values

Parameters	Units	Description	Value
$ J $	$ \frac{mol}{sec\cdot m^{2}} $	Diffusion influx	$-$
$ D $	$ \frac{m^{2}}{sec} $	Diffusion coefficient	$-$
$ C_{out} $	$ \frac{mg}{L} $	Umbelliferone concentration outside the cell	$ 100 $[1]
$ C_{in} $	$ \frac{mg}{L} $	Umbelliferone concentration inside the cell	$-$
$ z $	$ m $	Thickness of E. coli membrane	$ 4\cdot10^{-9} $[5]
$ A $	$ m^{2} $	Surface area of E. coli cell	$ 6\cdot10^{-12} $[4]
$$ \beta $$		Partition coefficient	$ 10^{1.896} $[6]
$ n $	$ mol $	Number of umbelliferone moles	$-$
$ V $	$ m^{3} $	Volume of E. coli cell	$ 10^{-18} $[4]
$ t $	$ sec $	Time
$ k $	$ \frac{m^{2}\cdot kg}{K\cdot sec^{2}} $	Boltzmann constant	$1.380649\cdot10^{-23}$[7]
$ T $	$ K $	Temperature	$ 305 $
$ R $	$ m $	Maximum radius of umbelliferone molecule	$ 5.06\cdot10^{-10} $ [6]
$ \eta $	$ Pa\cdot sec $	Viscosity of medium	$ 796\cdot10^{-6} $[33]

Assumptions:

---

1. The solutions on both sides of the membrane are homogeneous. By this, we mean that the solute is evenly mixed on either side.
2. The membrane does not change its width.
3. UMB molecule is a spherical particle with a radius of R.
4. UMB is diffusing through a liquid with a low Reynolds number suggesting a laminar flow of the medium.
5. Radius of a molecule (R) is a maximum projection radius describing the worst-case scenario of lower diffusion.

Fick's law for one-dimensional free diffusion^[3]

---

$$ (1.1)\ J=-D\frac{dC}{dx} $$ $$ (1.2) \ J=\frac{dn}{dt}\cdot\frac{1}{A}\rightarrow\frac{dn}{dt}=-DA\frac{dC}{dx} $$

By the stated assumptions 1 and 2, the above equations can be rewritten as follows:

$$ (1.3) \ \frac{dn}{dt}=-DA\frac{(C_{in}-C_{out})}{z} $$

To consider the membrane permeability it is necessary to consider the partition coefficient, which is the ratio of the solubility of the molecule in lipid and water and derived experimentally.

$$ (1.4) \ \frac{dn}{dt}=A\frac{\beta D}{z}(C_{out}-C_{in}) $$

Dividing both sides by V:

$$ (1.5) \ \frac{dn}{dt}\cdot\frac{1}{V}=A\frac{\beta D}{Vz}(C_{out}-C_{in})\rightarrow\frac{dC}{dt}=A\frac{\beta D}{Vz}(C_{out}-C_{in}) $$

Dividing both sides by $(C_{out}-C_{in})$:

$$ (1.6) \ \frac{dC}{dt}\cdot\frac{1}{(C_{out}-C_{in})}=A\frac{\beta D}{Vz} $$

Solving first order equation:

$$(1.7) \ \int\frac{dC}{(C_{out}-C_{in})}=\int A\frac{\beta D}{Vz}dt $$

Under the assumption that the only variable is $C_{in}$:

$$ (1.8) \ -ln(C_{out}-C_{in})=A\frac{\beta D}{Vz}\cdot t+a $$

Raising both sides of the equation by the power of exp:

$$ (1.9) \ (C_{out}-C_{in})=be^{-A\frac{\beta D}{Vz}\cdot t} $$
$$ (1.10) \ C_{in}=-be^{-A\frac{\beta D}{Vz}\cdot t}+C_{out}$$

When $b=e^{a}$

At the time 0 sec concentration of umbelliferone inside the cell is 0 $\frac{mol}{m^{3}}$:

$$(1.11) \ 0=-b+C_{out}\rightarrow b=C_{out} $$

Inputing equation (11) into equation (10):

$$ (1.12) \ C_{in}=C_{out}(-e^{-A\frac{\beta D}{Vz}\cdot t}+1) $$

The value for diffusion coefficient wasn't found, therefore we decided to evaluate it using the Stokes-Einstein Law for diffusion^[7]:

$$ (1.13) \ D=\frac{kT}{6\pi\eta R} $$

Inputting equation (13) into equation (12), by using assumption (5):

$$ (1.14) \ C_{in}=C_{out}(-e^{-A\frac{\beta kT}{6\pi\eta RVz}\cdot t}+1) $$

Inputting values into the final equation:

$$ (1.15) \ C_{in}=C_{out}(-e^{-12.4754\cdot t}+1) $$

According to the work by Bu, X.L et.al ^[1], we assume an initial concentration of $100 \frac{mg}{L}$ outside of the cell. Thus $\text{C}_{\text{out}} = 100 \frac{mg}{L}$.

Result and Conclusions

---

Using these values we get the following graph from our simulation:

The graph shows the concentration of umbelliferone inside the cell in mg/L vs. time in sec. In less than one second the concentration inside the cell reaches the initial concentration outside the cell. Therefore, the model has shown that diffusion occurs at a time scale of less than a second under all our assumptions (fig.2). Consequently we can use this conclusion and neglect the diffusion time in future models and developments.

2. The Fermentation Process and Enzyme Production

Introduction

---

We aimed to model the growth behavior of our bacteria in a fermentor. We used this model to predict the maximum amount of protein that will be produced by the bacteria. This amount will be significant in our enzyme kinetics model. To model the fermentation process we use the Monod growth kinetics^[9]. This theory links substrate (glucose), product (protein), and bacteria concentrations.

Definitions

---

Table 2: Description of Model Symbol and Values as Inferred from Literature

Symbol	Units	Description	Value
$ \mu $	$ \cfrac{1}{hour} $	Specific growth rate	$-$
$ X $	$ \cfrac{g_{bacteria}}{liter} $	Biomass concentration	$-$
$ t $	$ hour $	Time	$-$
$ P $	$ \cfrac{g_{protein}}{liter} $	Protein concentration	$-$
$ S $	$ \cfrac{g_{glucose}}{liter} $	Substrate (glucose) concentration	$-$
$ q_{p} $	$ \cfrac{g_{protein}}{g_{bacteria}\cdot hr} $	Specific protein formation	$-$
$ m_{s} $	$ \frac{g_{glucose}}{g_{bacteria}\cdot hr} $	Maintenance coefficient	$ 0.07 $ [9]
$ K $	$ \frac{1}{hr} $	Specific protein destruction	$ 0.01 $ [11]
$ \mu_{max} $	$ \frac{1}{hr} $	Maximum growth rate	$ 0.27 $ [12]
$ k_{s} $	$ \frac{g}{liter} $	Substrate concentration when $\mu=0.5\mu_{max}$	$ 180\cdot10^{-7} $
$ Y_{x/s} $	$ \cfrac{g_{bacteria}}{g_{glucose}} $	Yield of bacteria per substrate	$ 0.37 $ [13]
$ Y_{p/s} $	$ \cfrac{g_{protein}}{g_{glucose}} $	Yield of protein per substrate	$ 0.2 $ [14]
$ Y_{p/x} $	$ \cfrac{g_{protein}}{g_{bacteria}} $	Yield of protein per bacteria	$ 0.46 $ [12]

Assumptions

---

1. We assume that the most expensive part of the process will be the separation process between the bacteria and decursinol based on our consultation with Prof. Yuval Shoham. Under this assumption, we chose to use a batch process. This will allow a higher concentration of decursinol per liter in the manufacturing process and will reduce production costs
2. Under these conditions we assume the limiting factor will be glucose, as it will eventually run out.

Equations

---

We describe our system with a set of ordinary differential equations (ODEs)^[15]:

The following equation describes the change in biomass concentration over time:

$$ (2.1) \ \dfrac{dX}{dt}=\mu\cdot X$$

We use the Monod equation, which describes the specific growth rate:

$$ (2.2) \ \mu=\mu_{max}\cdot\frac{S}{K_{s}+S}$$

Then we define the following ODE that describes the change in protein concentration:

$$ (2.3) \ \dfrac{dP}{dt}=q_{p}\cdot X-K\cdot P$$

Calculate a constant that describes the specific protein formation:

$$ (2.4) \ q_{p}=Y_{p/x}\cdot\mu$$

The needed energy for glucose-consuming cellular activities can be described in terms of glucose (substrate) consumed per time. The energy required for cell growth:

$$ (2.5) \ \cfrac{\mu\cdot X}{Y_{x/s}}$$

The energy required for maintenance of the bacteria:

$$ (2.6) \ m_{s}\cdot X$$

The energy required for the synthesis of enzymes (product):

$$ (2.7) \ \cfrac{q_{p}\cdot X}{Y_{p/s}}$$

Adding it all together, the rate of glucose consumption is:

$$ (2.8) \ -\cfrac{dS}{dt}=\cfrac{\mu\cdot X}{Y_{x/s}}+m_{s}\cdot X+\cfrac{q_{p}\cdot X}{Y_{p/s}}$$

Results and Conclusions

---

As expected, in a batch process the glucose eventually runs out. (fig.4). Once there is no more glucose in the system, bacterial growth stops. That causes protein production to stop as well. Because there are no more new proteins, the relative effect of protein natural degradation increases and a reduction of total protein levels can be seen. The model predicts peaks of bacteria (fig.3) and protein (fig.5) concentrations at around 4 hours after induction. The predicted protein concentration is ~0.86g/L. The predicted protein concentration is on par with data from the literature for the BL21 strain, with an induction time similar to our experiment^[12] .

3. Directed Mutagenesis of the XimE Protein

In order to make the production of decursin cost-effective, throughout the modeling and biosynthetic design, we tried to focus on making the anabolic pathway more efficient. One of the ways used nowadays is directed protein mutagenesis, a robust method to design proteins with desirable functions, in our case greater production of decursin. Therefore, for this purpose computational tools were used to identify mutations that will improve our proteins and would help to shorten the time of future lab work.

The desirable target in the case of decursin production would be the remodeling of one of the key enzymes in the pathway, XimE. It is the final enzyme in our anabolic pathway, converting epoxidated 7-demethylsuberosin, (Fig.6), into decursinol. Using directed mutagenesis, we would like to achieve greater binding affinity of the substrate, which in turn, under our assumptions, leads to greater production of decursinol.
The computational approach of docking aims to predict the binding pose and affinity of a small molecule ligand and protein through scoring algorithms.

Method and Tools

---

Our workflow consisted of performing structural prediction of a protein with a known single point mutation, conducting molecular docking of predicted mutant proteins of the same type with unchangeable ligand, and eventually analyzing and comparing results to reveal the desired mutation and its further impact on wet lab work.

First, we found the 3D structure and amino-acid sequence of XimE in the Protein Data Bank (PDB) and manually performed point mutations in desired positions to build the stock of mutant sequences destined for structural predictions. While accessing AlphaFold2 from ChimeraX ^[16], we retrieved the predicted structures of mutant proteins of the same type (XimE). AlphaFold^[17] is a computational approach capable of predicting protein structures to near experimental accuracy, whereas ChimeraX^[16] runs its script automatically on Google Colab servers. During our trials with docking between target ligand and XimE we noticed great differences in scores between WT protein retrieved from PDB, aka crystallographic structure, and WT protein predicted with AlphaFold, although both have the same amino acid sequence. Therefore, we decided to use the predicted wild-type (WT) protein from AlphaFold for comparison with predicted mutant proteins for our further investigations.

The next step was protein-ligand docking. It is a method in structural molecular biology, whose goal is to predict the optimal binding mode of a ligand with a protein of known three-dimensional structure. There are various programs that can be used for protein-ligand docking, however, more comfortable the fastest and the most dependable for us appeared to be CB-Dock2 used for blind docking because the protein of interest lacked the verified position of the active site. The program first automatically detects binding sites and customizes the docking box and then employs the latest version (1.2.0) of AutoDock Vina for the docking itself.^[18][19] As a result, the software draws out the list of ranked binding sites of the ligand according to the Vina score, predicted binding affinity calculated in kcal/mol, which in turn can be translated into binding free energy ∆G.

Finally, we analyzed and compared all the scores and visualized the following interaction between mutant protein and target ligand using PyMol, a molecular visualization system on an open-source foundation.

You can find all of our files, and artifacts here: Supplements

Position of Mutations

---

XimE is a dimer made of two identical subunits, hence for the simplicity of the processes we ran the mutations on a single chain (chain a) as indicated in the literature^[21][22][23].

Before the molecular docking, we needed to choose the suitable amino acids on which the mutation would be performed. To understand the region of the cavity where the ligand binds, we ran a wild-type (WT) simulation with our ligand and received an estimated location. We integrated the amino acids around the selected area with the literature and selected three amino acids that were most examined in mutagenesis experiments. In their paper, Bei-Bei H et al.^[21] have identified Glutamic acid at position 136 (E136), Histidine at position 102 (H102), and Tyrosine at position 153 (Y153) as possible key residues for the pyran ring formation. Chuchu Jiang et al.^[22] stated that mutagenesis in E136 or H102 didn’t eliminate the enzyme activity, reinforcing the choice.

For Y153 we also based on the work of Xu-Liang Bu et al. ^[23] that performed in vitro mutagenesis in Y153.

Moreover, we performed the initial docking of the WT protein with the target ligand and received the following list of residues participating in the bond as listed in Fig.6. It can be shown that the three selected amino acids are displayed.

Results

---

The resulting scores obtained from CB-Dock2 were compared and summed in the supplements attached, to view them click here. The leading scores from each mutated amino acid are presented below in table 1.

It can be seen from the table that the initial score of the WT is -6.0 kcal/mol, and the favorable mutations’ scores are lower. This matches expectations as a more negative score indicate a higher affinity between the protein and the ligand^[24].

The best score was obtained from the mutation H102P, with a score of -9.2 kcal/mol.

Table 3: Top Vina scores for the mutagenesis of each of the selected amino acid.

Mutation	Vina score [kcal/mol]
WT	-6.0
Y153R	-7.8
E136P	-8.4
H102P	-9.2

(a)

(b)

(c)

Proline is a much more rigid amino acid than others, and due to its characteristics, it often improves the stability of the protein as a whole^[25], which may ensure a stronger and more stable bond between the XimE and the ligand.

Conclusion

---

Directed mutagenesis is a significant tool for improving the properties of proteins and the affinity between enzymes and their substrates. By using computational tools, it is possible to investigate a wide variety of mutations resulting in reducing time, manpower, and laboratory equipment. After targeting three amino acids- E136, Y153, and H102 we received different scores, with H102P being the lowest. As Proline is known for stabilizing tertiary structures of proteins it is a possible explanation as to why H102P is the optimal mutation.

Model Implications

---

The computational approach described here was desirable for us, because of the feasibility of testing the predictions and results in the wet lab. In our case, we were able to design primers that allow us to generate mutated enzymes based on the model results. By establishing the production rate of decursin as a criterion for comparison, we can assess the results of different mutations to the WT. The experiments we planned to validate our modeling work would be a good starting point for our future work.

4. The Enzymatic Reactions and Kinetics of PT, XimE, and XimD

Introduction

---

We define the system as an expression system with three enzymes: PT, XimE, and XimD. PT catalyzes the reaction from UMB to DMS and XimE together with XimD catalyzes the reaction from DMS to DEC. To optimize the system for the production of DEC, a precursor of decursin, we needed to model the enzymatic reactions and characterize their dynamics.

Estimation of Enzyme Concentration

---

Using our fermentation model we were able to estimate the total amount of protein produced in the system. According to our design, the bacterial population expresses all three enzymes. The next step was to understand the distribution between the different enzymes that catalyze the enzymatic reactions, this will allow us to calculate the concentration of each enzyme in the system. We used experimental data found by Bei-Bei et al. (2019), which held experiments on a two-enzyme system that expresses only XimE and XimD.

Now, we can model two different systems, plug-in experimental data from the two-enzyme system, and based on that draw conclusions regarding the three-enzyme system.

The three-enzyme system: a system with an expression of XimE, XimD, and PT.

The two-enzyme system: a system with an expression of XimE, and XimD.

Definitions

---

Table 4: Description of Model Symbol and Values as Inferred from Literature

Symbol	Units	Description	Value
$\text{XimD length}$	$\text{No. of amino acids}$	Amount of amino acids in a XimD protein	$474$ [27]
$\text{XimE length}$	$\text{No. of amino acids}$	Amount of amino acids in a XimE protein	$124$[27]
$\text{PT length}$	$\text{No. of amino acids}$	Amount of amino acids in a PT protein	$400$[27]
$N_{A}$	$ \frac{\text{molecules}}{\text{mol}} $	Avogadro constant	$\sim6.022\cdot10^{23}$
$\text{[Total Protein]}$	$ \frac{mg}{ml} $	Concentration of protein produced in the cell according to presented fermentation model [1]	$0.5$
$m_{\text{aa}}$	$Da$	Average weight of an amino acid	$110$ [26]
$MW_{\text{XimD}}$	$ Da $	Molar mass of a single XimD protein molecule	$51,010$[27]
$MW_{\text{XimE}}$	$ Da $	Molar mass of a single XimE protein molecule	$13,637$[27]
$m_{\text{PT}}$	$ Da $	Mass of a single PT protein molecule	$44,221$[27]
$v_{\text{cell}}$	$ \mu m^{3} $	Volume of a bacterium cell	$1$[29]
$m_{2,\text{Total XimD}}$	$ mg $	Total mass of XimD protein in the two enzyme system	Calculated in model
$m_{2,\text{Total XimE}}$	$ mg $	Total mass of XimE protein in the two enzyme system	Calculated in model
$m_{3,\text{Total XimD}} $	$ mg $	Total mass of XimD protein in the three enzyme system	Calculated in model
$m_{3,\text{Total XimE}}$	$ mg $	Total mass of XimE protein in the three enzyme system	Calculated in model
$[\text{XimD}]_{2} $	$ \mu M $	Enzyme concentration in an two enzyme system	Calculated in model
$ [\text{XimE}]_{2}$	$ \mu M $	Enzyme concentration in an two enzyme system	Calculated in model
$ [\text{XimD}]_{3}$	$ \mu M $	Enzyme concentration in an three enzyme system	Calculated in model
$[\text{XimE}]_{3} $	$ \mu M $	Enzyme concentration in an three enzyme system	Calculated in model
$ m_{\text{total protein}}$	$ mg $	Total mass of enzymes in a single bacterium	Calculated in model
$N$	$ - $	Number of molecules in volume V	Calculated in model
$V$	$ \mu M $	Total mass of enzymes in a single bacterium	Volume

Assumptions

---

1. The total protein mass in an system with $n$ enzymes is distributed among the proteins introduced in our plasmid systems:$ m_{\text{total protein}}={\displaystyle \sum_{i=1}^{n}m_{\text{total protein i}}}$
2. Our fermentation model suggests peak protein concentration of $0.86\frac{mg}{ml}$, this is supported by the literature: “For cells induced at around $OD6oo = 0.6$ and harvested 3 h later, a typical yield is about $0.5-1mg$ total protein per $ml$ of culture”^[25]. (OD refers to the cell density). We decided to take a conservative value: $[\text{Total Protein}]=0.5\frac{mg}{ml}$.
3. The ratio between the weights of both proteins is the ratio between their respective lengths.
4. Assume equality between average mass [$Da$] and molar mass $[\frac{g}{mol}]$
5. We use an estimation $\frac{1}{\mu m^{3}}=1nM$ following this approximation:
$$C=\frac{N}{AV}=\frac{1}{6\cdot10^{23}\mu m^{3}}\approx\frac{1}{10^{24}mol^{-1}\cdot1\mu m^{3}}=$$ $$=\frac{1mol}{10^{24}(10^{-6}m)^{3}}=\frac{1mol}{10^{6}m^{3}}=\frac{1mol}{10^{6}\cdot1000liter}=\frac{1mol}{10^{9}l}=10^{-9}M=1nM$$

Two Enzyme System: XimE & XimD

---

Because our expression system contains two enzymes, it holds that:

$$ (4.1)\ m_{\text{total protein}}=[\text{Total Protein}]\cdot v_{\text{cell}}=m_{\text{XimE}}+m_{\text{XimD}} $$

By our assumption, the distribution of mass is equal to the ratio between enzyme lengths:

$$ (4.2)\ \frac{m_{2,\text{Total XimD}}}{m_{2,\text{Total XimE}}}=\frac{\text{XimDlength}}{\text{XimE length}}$$

This yields:

$$ (4.3)\ m_{2,\text{Total XimD}}=\frac{\text{XimD length}}{\text{XimD length}+\text{XimE length}}\cdot m_{\text{total protein}} $$ $$ (4.4)\ m_{2,\text{Total XimE}}=\frac{\text{XimE length}}{\text{XimD length}+\text{XimE length}}\cdot m_{\text{total protein}}$$

We calculate the number of molecules from each kind to get the concentration of each enzyme in a single bacterium.
For XimD:

$$ (4.5)\ [\text{XimD}]_{2}=\frac{m_{2,\text{Total XimD}}\cdot\frac{N_{A}}{MW_{\text{XimD}}}}{v_{\text{cell}}}=\frac{4600}{1\mu m^{3}}\underbrace{=}_{\text{(Assumption 5)}}4.6\,\mu M $$

For XimE:

$$ (4.6)\ [\text{Xime}]_{2}=\frac{m_{2,\text{Total XimE}}\cdot\frac{N_{A}}{MW_{\text{XimE}}}}{v_{\text{cell}}}=\frac{4400}{1\mu m^{3}}\underbrace{=}_{\text{(Assumption 5)}}4.4\,\mu M $$

Three Enzyme System: PT & XimE & XimD

---

In this case:

$$ (4.7)\ m_{\text{total protein}}=m_{\text{XimE}}+m_{\text{XimD}}+m_{\text{PT}}=0.5mg $$

In a similar fashion, considering PT Length:XimD Length:XimE Length=3.2:3.8:1, we find:

$$ [\text{XimD}]_{3}=2.7\,\mu M $$ $$ [\text{XimE}]_{3}=2.8\,\mu M $$ $$ [\text{PT}]_{3}=2.7\,\mu M $$

Enzyme Kinetics Model

---

After estimating the concentration of the enzymes in our system we continued our modeling work to answer two questions, the first being how much of the UMB is converted to DEC. The second would be the temporal aspect: understanding how long the enzymatic reaction will take. Both of these aspects together would allow us to quantify a production rate. This understanding of the key factors in the biomanufacturing process will allow us to create a more optimized biological circuit design, and accordingly provide us with predictions for the cost and time of the production.

Michaelis–Menten Kinetics

Definitions

---

Table 5: Description of Model Symbol and Values as Inferred from Literature

Symbol	Units	Description	Value
$\text{[XimE]}$	$\mu M$	XimE enzyme concentration	$2.8 \text{[previous model]}$
$\text{[XimD]}$	$\mu M$	XimD enzyme concentration	$2.7 \text{[previous model]}$
$\text{[PT]}$	$\mu M$	Prenyltrasnfarase enzyme concentration	$2.7 \text{[previous model]}$
$v_{\text{UMB}}$	$\frac{\mu M}{sec}$	Rate of formation of UMB	$-$
$v_{\text{DMS}}$	$\frac{\mu M}{sec}$	Rate of formation of DMS	$-$
$v_{\text{DEC}}$	$\frac{\mu M}{sec}$	Rate of formation of DEC	$-$
$v_{\text{MAR}}$	$\frac{\mu M}{sec}$	Rate of formation of marmesin	$-$
$k_{m,\text{PT}\cdot\text{UMB}}$	$\mu M$	Michaelis constant	$20$[30]
$k_{1,\text{cat}}$	$\frac{1}{sec}$	Rate of formation DMS from UMB	$0.435$[30]
$[\text{UMB}]$	$\mu M$	UMB concentration	$-$
$[\text{DMS}]$	$\mu M$	DMS concentration	$-$
$[\text{DEC}]$	$\mu M$	DEC concentration	$-$
$[\text{MAR}]$	$\mu M$	Marmesin concentration	$-$

Assumptions

---

The entirety of the metabolite concentration is consumed in the reaction and converted to the products defined in each reaction.

[UMB] in the cell is constant, due to the diffusion being immediate as seen in our diffusion model

Equations

---

We defined the following reaction:

$$ (4.8)\ \text{PT}+\text{UMB}\overset{k_{1}}{\underset{k_{-1}}{\rightleftarrows}}\text{PT}\cdot\text{UMB}\overset{k_{1cat}}{\rightarrow}\text{PT}+\text{DMS} $$

based on the Michaelis-Menten equation^[30] which states that the rate of a reaction is proportional to the concentrations of the enzyme and the substrate:

$$ (4.9)\ v=\frac{d[P]}{dt}=\frac{k_{\text{cat}}\cdot[S]\cdot[E]}{k_{m}+[S]} $$

we reached the following:

$$ (4.10)\ v_{\text{UMB}}=-\frac{k_{1,\text{cat}}\cdot[\text{UMB}]\cdot[\text{PT}]}{k_{m,\text{PT}\cdot\text{UMB}}+[\text{UMB}]} $$

Gibbs Distribution Model

---

The penultimate reaction in the biosynthesis transforms DMS to DEC. It is catalyzed by two enzymes co-dependently: flavin-dependent monooxygenase (XimD) and SnoaL-like cyclase (XimE). Together, they catalyze the formation of the pyran ring (6-endo-tet cyclization) necessary for the synthesis of pyranocoumarins such as decursinol.

In vitro enzymatic assays were used ^[1] to show that XimD alone (without the presence of XimE) catalyzes the formation of furan ring (5-exo-tet cyclization) which produces marmesin (MAR) instead of DEC.

It is only when both XimE and XimD are present, that the pyran ring is formed and decursinol is produced.

Microbial host models, expressing both XimE and XimD, clarified that the two reaction paths occur simultaneously. As we wanted to optimize our system for the production of DEC, we wanted the system to prefer the reaction that produces DEC

A suggested mechanism proposes an unstable epoxide intermediate, that can spontaneously decay to marmesin in the absence of XimE.

Definitions

---

Table 6: Description of Model Symbol and Values as Inferred from Literature

Symbol	Units	Description	Value
$N$	$-$	Number of molecules	$-$
$c$	$\mu M$	Concentration of reaction complex	$-$
$c_{0}$	$\mu M$	Initial concentration	$-$
$\beta$	$\frac{1}{kJ}$	Thermal influence	$\frac{1}{k_{B}T}$
$k_{B}$	$\frac{kJ}{K}$	Boltzmann constant	$1.380649\cdot10^{-26}$
$T$	$K$	Temperature	$-$
$\Delta\varepsilon$	$kJ$	Binding energy	$-$
$\varepsilon$	$-$	A constant expressing the dependency between XimE and XimD	In the magnitude of $K_{D,3}$
$Z$	$-$	Partition function	$-$
$k_{2,cat}$	$\frac{1}{s}$	Rate of formation marmesin from DMS	$0.1$ (Estimating model coefficients)
$k_{3,cat}$	$\frac{1}{s}$	Rate of formation DEC from DMS	$0.44$ (Estimating model coefficients)
$K_{D,1}$	$\mu M$	Dissociation constant of XimD with DMS	$0.1$ (Estimating model coefficients)
$K_{D,2}$	$\mu M$	Dissociation constant of XimE,XimD system with DMS	$0.1$ (Estimating model coefficients)

Statistical Equations

---

First, we defined the following reaction:

$$ (4.11)\ \text{XimD}+\text{DMS}\overset{k_{2cat}}{\rightarrow}\text{MAR} $$ $$ (4.12)\ \text{XimD}+\text{DMS}+\text{XimE}\overset{k_{3cat}}{\rightarrow}\text{DEC}+\text{MAR} $$

We use a state variable description approach of protein-substrate binding. ^[31]

Each substrate molecule can be bound to each one of the two enzymes (XimD and XimE). This gives us four optional states:

1. No reaction: free DMS without XimD or XimE
2. Exclusively XimD: XimD is bound to DMS and catalyzes DMS to MAR
3. XimD and XimE: Both catalyze the reaction by binding to the DMS and catalyze DMS to DEC
4. Exclusively XimE: Does not have a catalytic activity alone (because it binds to an intermediate created by XimD). Therefore, we express this dependency between the enzymes with a new parameter called $\varepsilon\sim K_{D,3}>>K_{D,1}$

Thus, we define a partition function. This function takes into account every state that our system could be in:

$$ (4.13)\ Z=1+\frac{[\text{XimD}]}{K_{D,1}}+\frac{[\text{XimD}][\text{XimE}]}{K_{D,1}\cdot K_{D,3}}\varepsilon+\frac{[\text{XimE}]}{K_{D,3}} $$

While we perform the following reduction:$ \frac{\varepsilon}{K_{D,1}\cdot K_{D,3}}=\frac{1}{K_{D,2}}$

According to the Gibbs model^[3], the average number of molecules in a bounded state is:

$$ (4.14)\ N=c\frac{\frac{1}{c_{0}}e^{-\beta\Delta\varepsilon}}{Z} $$

To make our model depend on a single parameter, we define the following notation:

$$ (4.15)\ \frac{1}{\frac{1}{c_{0}}e^{-\beta\Delta\varepsilon}}=c_{0}e^{\beta\Delta\epsilon}=K_{D} $$

We define the rate of a reaction as the number of molecules in bound state times the catalytic rate constant:

$$ (4.16)\ v=\frac{d[P]}{dt}=k_{cat}\cdot N $$

And plug (4.6) and (4.7) into (4.8) to receive:

$$ (4.17)\ v=\frac{k_{cat}\cdot c}{K_{D}\cdot Z} $$

Then, the change in the concentration of the metabolites follows:

$$ (4.18)\ v_{\text{DMS}}=\frac{k_{1,\text{cat}}\cdot[\text{UMB}]\cdot[\text{PT}]}{k_{m,\text{PT}\cdot\text{UMB}}+[\text{UMB}]}-\frac{k_{2,\text{cat}}}{K_{D,2}}\frac{[\text{DMS}]\cdot[\text{XimD}]\cdot[\text{XimE}]}{Z}-\frac{k_{3,\text{cat}}}{K_{D,1}}\frac{\cdot[\text{DMS}]\cdot[\text{XimD}]}{Z} $$ $$ (4.19)\ v_{\text{DEC}}=\frac{d[\text{DEC]}}{dt}=\frac{k_{2,\text{cat}}}{K_{D,2}}\frac{[\text{DMS}]\cdot[\text{XimD}]\cdot[\text{XimE}]}{Z} $$ $$ (4.20)\ v_{\text{MAR}}=\frac{d[\text{MAR]}}{dt}=\frac{k_{3,\text{cat}}}{K_{D,1}}\frac{\cdot[\text{DMS}]\cdot[\text{XimD}]}{Z} $$

Where the change in DMS concentration is defined as:

$$ (4.21)\ v_{\text{DMS}}=-v_{\text{UMB}}-v_{\text{DEC}}-v_{\text{MAR}} $$

Estimating Model Coefficients

---

To estimate the coefficients, we use the results found by Bei-Bei et al^[30]: Cloning of XimD and XimE under two pET plasmids resulted in a DEC:MAR ratio of 1:1. Assuming the ratio of $\frac{dp}{dt}$ correlates with the final concentration, we get:

$$ (4.22)\ 1=\frac{[\text{DEC}]_{\text{final}}}{[\text{MAR}]_{\text{final}}}\implies1=\frac{v_{\text{DEC}}}{v_{\text{MAR}}}=\frac{\frac{k_{2},cat}{K_{D,2}}}{\frac{k_{3,cat}}{K_{D,1}}}\cdot[\text{XimE}]_{\text{basal}}\implies\frac{\frac{k_{2},cat}{K_{D,2}}}{\frac{k_{3,cat}}{K_{D,1}}}=\frac{1}{[\text{XimE}]_{2}} $$

Define:

$$\frac{k_{2,cat}}{K_{D,2}}=k_{\text{eff,2}}$$ $$\frac{k_{3,cat}}{K_{D,1}}=k_{\text{eff,3}}$$

As mentioned before:

$$ \frac{k_{3,cat}}{K_{D,1}}=k_{\text{eff},3}=4.4\frac{k_{2,cat}}{K_{D,2}}=4.4k_{\text{eff},2} $$

Therefore:

$$\frac{k_{3,cat}}{K_{D,1}}=k_{\text{eff,3}}=4.4\frac{k_{2,cat}}{K_{D,2}}=4.4k_{\text{eff,2}}$$

For that reason we chose arbitrary values while maintaining 2 conditions:

1. The ratio $\frac{k_{\text{eff,3}}}{k_{\text{eff,2}}}=4.4$
2. We keep the values to be in the same order of magnitude as $k_{1,\text{cat}}$

Results and Conclusions

---

Using these values we get the following graph from our simulation:

Figure 11a shows, that the entirety of our precursor, UMB is converted to the products, and results in saturation of DEC and MAR. The model also shows the feasibility of the process, the calculated production rate of DEC is 0.3 grams per liter per hour. The model predicts that the expression of the three enzymes under the same promoter will result in a higher concentration of the byproduct marmesin compared to the target product decursinol (fig.11a). Moreover, we observe that when the molar fraction of XimE increases, a higher concentration of decursinol is obtained compared to marmesin (fig.11b). This prediction prompted us to modify the design of our plasmid, we included a constitutive expression of XimE and an inducible expression of the other two enzymes. By controlling the expression levels of XimD & PT, we ensure a surplus of XimE. There is no correlation between [DEC] [MAR] and [XimD]. The ratio between the metabolites is proportional to the concentration of XimE.

Fitting the model to experimental data

Definitions

Table 6: Description of this Section Symbols and Values as Taken from qPCR Results

Symbol	Units	Description	Value
$\text{XimE}_{\text{Copy number}}$	$-$	Number of XimE transcripts	$-$
$\text{XimD}_{\text{Copy number}}$	$-$	Number of XimD transcripts	$-$
$\text{XimE}_{\text{CT}}$	$-$	Cycle number in which XimE cDNA level reached threshold of 0.3	$18.9810$
$\text{XimD}_{\text{CT}}$	$-$	Cycle number in which XimD cDNA level reached threshold of 0.3	$16.2080$

Assumptions

---

RNA level (measured by qPCR) is equivalent to the protein levels.

Equations

---

According to our design, we have a plasmid containing XimD and XimE under T7 promoters (ximE_ximD). The expression of XimD is regulated by the Lac repressor while the expression of XimE is constitutive. We measured the expression levels in our two-enzyme system in the wet lab. This was done by running a qPCR of samples containing our engineered bacteria under the induction of IPTG. For detailed results please refer to our results. To calculate the ratio between the expression levels of XimD and XimE, we used the threshold defined for the qPCR process. For each gene, we calculated a cycle threshold (CT) value in which it reached the threshold of fluorescence. At these CT values, the cDNA level is equal.

Thus:

$$ (4.23)\ \text{XimE}_{\text{Copy number}}\cdot2^{\text{XimE}_{\text{CT}}}=\text{XimD}_{\text{Copy number}}\cdot2^{\text{XimD}_{\text{CT}}}$$ $$ (4.24)\ \frac{\text{XimD}_{\text{Copy number}}}{\text{XimE}_{\text{Copy number}}}=\frac{2^{\text{XimE}_{\text{CT}}}}{2^{\text{XimD}_{\text{CT}}}}=\frac{2^{18.9810}}{2^{16.2080}}\approx6.8$$

6.8 is the molar ratio between XimD and XimE when IPTG is present. To get mass ratio we will multiply 6.8 with the ratio between the mass of a single XimD and the mass of a single XimE:

$$ (4.25)\ \text{XimD}:\text{XimD}=6.8\frac{\text{mol}}{\text{mol}}\cdot\frac{MW_{\text{XimD}}}{MW_{\text{XimE}}}=6.8\cdot\frac{51,010}{13,637}\approx25.4$$

consequently, we deduce that XimD:XimE is 25.4:1 when IPTG is present. We redo the process detailed in the two-enzyme system, and get:

$$[XimE]_{2}=0.23\mu M$$ $$[XimD]_{2}=5.5\mu M$$

We ran the model again with the new concentrations. However, now we are modeling a two enzymes system which means we remove equation 4.10 from our ODE system.

Figure 12 shows a very high concentration of the unwanted product MAR. These wet lab results greatly emphasize to us the importance of creating a construct where the concentration of XimE will be as high as possible. We discuss ways to achieve this goal in the model implications.

Model Implications

---

After studying the results of the different models, we can proceed and test the predictions in the wet lab. To test our predictions regarding how different concentrations of XimE and XimD affect the end products of the manufacturing.
We planned several plasmid designs, which can be compared according to the final concentrations of the end products.
Each design allows for different expression rates of the XimE and XimD enzymes:

1. XimE: By replacing the RBS of XimE, with an RBS designed to have a higher expression rate with the RBS calculator ^[32], we expect to have a higher final concentration of decursinol. This can also lead to the development of the model to take the RBS strength into account.
2. XimD: By replacing either the promoter or RBS of XimD with a weaker candidate, we can asses our prediction regarding the lack of correlation between $\frac{[\text{DEC}]}{[\text{MAR}]}$ and $[\text{XimD}]$.
3. Varying concentrations of inducers: by including a repression system, we could use the different inducer concentrations to conduct experiments with varying enzyme concentrations under the same plasmid. This will provide a high-throughput method of testing our model.

In the scope of iGEM, we focused on approach 3, which was the first one to be implemented in our design.