Team:SCU-China/Enzyme

MODELING

SCU-China Header


ENZYME KINETIC SIMULATION

a

Introduction

Our project this year concentrates mainly on the metabolic circuit of cordycepin (COR) and pentostatin (PTN). It involves three enzymes, Cns1, Cns2 and Cns3-HisG domain. Since all our work is related to these three enzymes, it is necessary to depict the overall picture of the enzymatic system.

Generally, our aim of constructing the enzyme kinetic simulation is divided into two parts. Firstly, we would like to theoretically give the model that testifies the necessity of using fusion expression of Cns1 and Cns2 in experiments, which can give powerful support and confidence to the later experiments of fusion expression. Additionally, we would like to discover how the diffusion rate influences the production rate. Secondly, since our goal is to co-ferment cordycepin and pentostatin furthest, we have to build a model that consists of all the elements involved in the metabolic circuit of cordycepin and pentostatin. More importantly, we need to find the correlation between the concentration of cordycepin and pentostatin, as well as give the optimal substrate ratio and optimal expression strength of the enzymes.

Fusion Expression Simulation

Our task is to deal with the circuit below. In the circuit, Cns2 converts 3'-AMP to the enol state of 2'-C-3'-dA, which then naturally converts to the keto state. Enzyme cns3 produces cordycepin using the product from Cns2.

Figure. 1 Enzyme catalytic circuit of cordycepin production, involving the enzyme Cns1 and Cns2. Enol means the enol form of the substance, and keto means the keto form.

This modeling section aims to compare the efficiency of the enzymatic reactions using a fusion-expression strategy to that using a separation-expression strategy, which in result will give theoretical support for the corresponding experiments.

It is reported in the literature that there exists interaction between Cns1 and Cns2 in vivo by yeast two-hybrid analysis (Xia Y et al., 2017). The phenomenon of interaction between two or more enzymes is not peculiar. It is common to see some multienzyme complex in nature, including pyruvate dehydrogenase, rubisco and even respiratory electron-transport chain, which is recently proven to be a megacomplex. The common consensus about the multienzyme complex is that it allows segregating certain biochemical pathways into one place in the cell, which increases the catalytic efficiency and reduces side effects. In nature, the expression of enzymes in the multienzyme complex is generally independent, after which they meet and combine into a multienzyme complex.

In the theory of structural biology, the fusion-expression strategy connects the two enzymes during expression so that the substrate of the first enzyme is sequentially converted to the product of the second enzyme, which is more efficient than the separation-expression strategy. However, it has not been mathematically proven in enzyme kinetics. Additionally, the evaluation of the linker connecting the enzymes is obscure before doing specific enzyme-activity experiments.

Here, we would like to build an enzyme kinetic model to simulate the separation-expression strategy and the fusion-expression strategy. Determination of the linker is conducted in computational structural biology model.

It is reported that keto-enol tautomerism can be simulated as a first-order reaction (Huang L et al., 2007) and can be represented as:

dCdt=keqC\frac{dC}{dt} = -k_{eq}C

in which CC is the concentration of the keto form.

  • General Assumption: Our enzyme kinetic simulation is assumed in vivo, so all the enzymatic reaction is reversable and 2'-C-3'-dA is considered to be all in keto form (keq=[Enolform][Ketoform]1k_{eq} = \frac{[Enol form]}{[Keto form]} \ll 1).

Distance between molecules in solution

The difference between these two strategies is the relative concentration of the substrate of the second enzyme, since the concentration is lower in seperate-expression strategy because the product of the first enzyme (which is at the same time the substrate of the second enzyme) need to diffuse to the second enzyme.

Considering Cns1 and Cns2 enzyme diffusing in the solvent, the initial concentration of the enzymes is uniform in the whole system. When the substrate of the first enzyme is converted into the product, it needs to diffuse to the second enzyme.

Therefore, the diffusing distance is the key issue in our simulation of separation expression.

Firstly, we should consider the distance

  • Assumption 1: Reaction environment can be treated as solution and the solution is uniform, which means that the concentration of the solvent at one position is the same with that at any other positions in the solution.
  • Assumption 2: All the enzymes can be assumed as spheres.

The average distance between molecules in a solution and the average radius of a molecule is calculated by Kozanecki et al. (Kozanecki M et al., 2016)

It is frequently useful to know the average volume of solution occupied by each molecule, or more directly, the average distance separating molecules in solution. This is a simple calculation based only on the molar concentration. In a 1 M solution, there are 6×10236×10^{23} molecules/L, = 0.6 molecules/nm3\text{nm}^3, or the volume per molecule is V=1.66 nm3nm^3/molecule at 1 M. For a concentration C, the volume per molecule is V=1.66/C. So, the distance between molecule in average is:

d=V13=1.18/C13d = V^{\frac{1}{3}} = 1.18/C^{\frac{1}{3}}

What we really want is a physically intuitive parameter for the
size of the protein. If we assume the protein has the simplest
shape, a sphere, we can calculate its radius. We will refer to this
as RminR_\text{min}, because it is the minimal radius of a sphere that could contain the given mass of protein:

Rmin=(3V/4π)13=0.066M13(for M in Dalton,Rminin nanometer)R_\text{min} = (3V/4\pi)^{\frac{1}{3}} = 0.066M^{\frac{1}{3}} (\text{for M in Dalton}, R_{min} \text{in nanometer})

Since the concentration of the enzyme in the cell is in the μmol\mu \text{mol}
magnitude, the distance between protein molecules is in the hundred-nanometer magnitude. In our project, the mass of Cns1, Cns2, Cns3 is around 80, 40, 100 kD respectively, so the radius of the enzyme is around 3 nm, which is in the unit-nanometer magnitude. Compared to the distance between enzymes, the radius is negligible.

  • Assumption 3: In the simulation of the fusion-expression, the enzyme is treated as a point, rather than a sphere.

To calculate the distance between Cns1 and Cns2 in the mixed solution of the system, we make another assumption to simplify the model.

  • Assumption 4: The two enzymes are evenly distributed in the whole solution.

Thus, if we assume the concentration of Cns1 to be Ccns1C_\text{cns1} , and concentration of Cns2 to be Ccns2C_\text{cns2} , there are totally (Ccns1+Ccns2)×6.02×1023(C_\text{cns1}+C_\text{cns2})\times 6.02\times 10^{23} molecules in 1 L system. As the calculation above, the distance between two enzymes can be represented by:

dcns1,cns2=1.18/(Ccns1+Ccns2)13d_{cns1,cns2} = 1.18/(C_{cns1}+C_{cns2})^{\frac{1}{3}}.

Compared to the separation -expression strategy, the fusion-expression strategy enables relatively no distance between Cns1 and Cns2, as they are connected using the linker.

Diffuson Simulation

During the diffusion process, the concentration of the substrate of the second enzyme is related to the product of the first enzyme by the diffusion modelge. Since all the molecule are evenly distributed, the three-dimensional system can be simplified to one-dimensional system because we only need to consider the relation between one Cns2 enzyme and the nearest Cns1 enzyme.

Here, we use Fick's law to describe the diffusion process.

  • Assumption 5: A concentration gradient is the only force driving substrate motion in the system, ignoring the influence of charge interactions or the shape of the molecule. (Eun C et al., 2014)

According to Fick's first law of diffusion, (Rice S et al., 1985)

J(z,t)=Dc(z,t)zJ(z,t) = -D \frac{\partial c(z,t)}{\partial z}

  • JJ is the diffusion flux, of which the dimension is amount of substance per unit area per unit time. JJ measures the amount of substance that will flow through a unit area during a unit time interval.
  • DD is the diffusion coefficient or diffusivity. Its dimension is area per unit time.
  • cc is the concentration, of which the dimension is amount of substance per unit volume.
  • zz is position, the dimension of which is length.

In the equation, the minus sign indicates that the flow occurs from the concentrated to the dilute region of the solution. (Macheras P et al., 2006)

If we treat the situation in a discretized way, the Fick's first law of diffusion can be interpreted in:

J(z,t)=DΔc(z,t)Δz=Dc(z+1,t)c(z,t)(z+1)zJ(z,t) = -D \frac{\Delta c(z,t)}{\Delta z} =-D\frac{c(z+1,t)-c(z,t)}{(z+1)-z}

Therefore, the flux of a position is related to the concentration gradient of this position.

When we consider the change of the concentration as the function of the time, we can get:

c(z,t)t=xJ=x(Dxc(z,t))=D2c(z,t)x2\frac{\partial c(z,t)}{\partial t} = - \frac{\partial}{\partial x} J = \frac{\partial}{\partial x}(D\frac{\partial}{\partial x} c(z,t)) = D\frac{\partial ^2 c(z,t)}{\partial x^2}

This formula is called Fick's second law of diffusion. Given the initial condition that c(0,0)=cl(t)c(0,0) = c_l(t) and c(δ,0)=cr(t)c(\delta,0) = c_r(t), the differential ordinary equation can be solved as the form below. (M.H. Jacobs., 1935)

c(z,t)=cl(t)+[cr(t)cl(t)]zδ+2πi=1cr(t)cosiπcl(t)isin(iπzδ)exp(Dδ2i2π2t)c(z,t) = c_l(t)+[c_r(t)-c_l(t)]\frac{z}{\delta}+\frac{2}{\pi}\sum_{i=1}^{\infty}{\frac{c_r(t)cosi\pi -c_l(t)}{i}}sin(i\pi \frac{z}{\delta})exp(-\frac{D}{\delta^2}i^2\pi^2 t)

In this equation, cr(t)c_r(t) is the substrate concentration of the second enzyme around the second enzyme, while cl(t)c_l(t) is the substrate concentration of the second enzyme in the context of the first enzyme. These two concentrations are independent to the diffusion process, but dependent on the production and consumption of substance.

If we fix zz to δ\delta, we can get:

c(z=δ,t)=cr(t)+2πi=1cr(t)cosiπcl(t)isin(iπ)exp(Dδ2i2π2t)c(z=\delta,t) = c_r(t)+\frac{2}{\pi}\sum_{i=1}^{\infty}{\frac{c_r(t)cosi\pi -c_l(t)}{i}}sin(i\pi )exp(-\frac{D}{\delta^2}i^2\pi^2 t)

Since sin(iπ)=0sin(i\pi) = 0,
c(z=δ,t)=cr(t)c(z=\delta,t) = c_r(t)

In this case, this equation is meaningless at the position of the second enzyme, but if we consider the position near the second enzyme, which makes the sin(iπ)0sin(i\pi) \neq 0 but very small. Based on this manipulation, the equation above is meaningful. Since the equation is too complicated to fit in our case, we'd like to simplify it to this form.

c(z=δ,t)=cr(t)kdiffcr(t)cl(t)exp(δ)c(z=\delta,t) = c_r(t) - k_\text{diff}\frac{c_r(t)-c_l(t)}{exp(\delta)}

Since δ\delta is the distance between the source and the point of interest, if we do not fix δ\delta, but set it to a parameter xx, we can get:

c(x,t)=cr(t)kdiffcr(t)cl(t)exp(x)c(x,t) = c_r(t) - k_\text{diff}\frac{c_r(t)-c_l(t)}{exp(x)}

Because cr(t)c_r(t) and cl(t)c_l(t)
is independent to the diffusion process, they can be represented to cr(t)c_r(t) and cl(t)c_l(t) , which is consistent to the equations in enzyme kinetics. c(x,t)c(x,t) is the updated concentration of the substance near the second enzyme which is affected by diffusion. So the equation looks like:

dcrdt=kdiffclcrexp(x)\frac{dc_r}{dt} = k_\text{diff}\frac{c_l-c_r}{exp(x)}

This equation describes that the derivation of the concentration in the low concentration area is proportional to the difference between the concentrations and nonlinearly inversely proportional to the distance xx. The coefficient kdiffk_\text{diff} represents the diffusion rate of substances in the solution. Since we have no access to the value of this parameter, we set it to an optional value and investigate the influence on the model by this value.

Enzyme Kinetics Simulation of Fusion Expression

Fusion Expression Simulation with limited resource

Fusion-expression strategy enables the enzyme Cns1 and Cns2 to be expression closely in distance. The product of the first enzyme instantaneously functions as the substrate of the second enzyme and is converted to the end product, cordycepin.

Based on all the assumptions and prerequisite derivation. we can create the following diagram based on typical enzyme kinetics:

Figure. 2 Elucidation of the primary reactions in enzymatic kinetics. S1 represents the substrate of Cns2, 3’-AMP. P1 is the product of Cns2, 2’-C’3’-dA, and it is the substrate of Cns1. P2 represents the product of Cns1, cordycepin. Since we consider the reactions in vivo, all the reactions are reversible.

Based on mass action law, we can now list the equations:

d[S1]dt=k1[Ecns2][S1]+k1[Ecns2S1]\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1]

d[Ecns2]dt=k1[Ecns2][S1]+k1[Ecns2S1]+k2[Ecns2S1]k2[P1][Ecns2]\frac{d[E_\text{cns2}]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1] + k_2[E_\text{cns2}S_1] - k_{-2}[P_{1}][E_\text{cns2}]

d[Ecns2S1]dt=k1[Ecns2][S1]k1[Ecns2S1]k2[Ecns2S1]+k2[P1][Ecns2]\frac{d[E_\text{cns2}S_1]}{dt} = k_1[E_\text{cns2}][S_1] - k_{-1}[E_\text{cns2}S_1] - k_2[E_\text{cns2}S_1] + k_{-2}[P_{1}][E_\text{cns2}]

d[P1]dt=k2[Ecns2S1]k2[P1][Ecns2]k3[Ecns1][P1]+k3[Ecns1S2]\frac{d[P_1]}{dt} = k_2[E_\text{cns2}S1]-k_{-2}[P_1][E_\text{cns2}]-k_3[E_\text{cns1}][P_1]+k_{-3}[E_\text{cns1}S_2]

d[Ecns1]dt=k3[Ecns1][P1]+k3[Ecns1S2]+k4[Ecns1S2]k4[P2][Ecns1]\frac{d[E_\text{cns1}]}{dt} = -k_{3}[E_\text{cns1}][P_{1}] + k_{-3}[E_{cns1}S_2] + k_4[E_{cns1}S_2] - k_{-4}[P_2][E_{cns1}]

d[Ecns1S2]dt=k3[Ecns1][P1]k3[Ecns1S2]k4[Ecns1S2]+k4[P2][Ecns1]\frac{d[E_\text{cns1}S_2]}{dt} = k_3[E_\text{cns1}][P_{1}] - k_{-3}[E_\text{cns1}S_2] - k_4[E_\text{cns1}S_2] + k_{-4}[P_2][E_\text{cns1}]

d[P2]dt=k4[Ecns1S2]k4[P2][Ecns1]\frac{d[P_2]}{dt} = k_4[E_\text{cns1}S_2] - k_{-4}[P_2][E_\text{cns1}]

Based on the derivation of Michaelis-Menten equation, from the equation of d[Ecns2]dt\frac{d[E_\text{cns2}]}{dt} and d[Ecns2S1]dt\frac{d[E_\text{cns2}S_1]}{dt}, we can find that d[Ecns2]dt+d[Ecns2S1]dt=0\frac{d[E_\text{cns2}]}{dt}+\frac{d[E_\text{cns2}S_1]}{dt} = 0, which means that the value of [Ecns2]+[Ecns2S1][E_\text{cns2}]+[E_\text{cns2}S_1] is a constant. Here, we defined this constant as Eini2E_\text{ini2}, which represents the constant of the initial concentration of cns2 enzyme. Likewise, considering to cns1 enzyme, value of [Ecns1]+[Ecns1S2][E_\text{cns1}]+[E_\text{cns1}S_2] is a constant, which is defined as Eini1E_\text{ini1}.

Therefore, the equations convert to:

d[S1]dt=k1[Ecns2][S1]+k1(Eini2[Ecns2])\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}(E_\text{ini2}-[E_\text{cns2}])

d[Ecns2]dt=k1[Ecns2][S1]+(k1+k2)(Eini2[Ecns2])k2[P1][Ecns2]\frac{d[E_\text{cns2}]}{dt} = -k_1[E_\text{cns2}][S_1] + (k_{-1}+k_2)(E_\text{ini2}-[E_\text{cns2}]) - k_{-2}[P_{1}][E_\text{cns2}]

d[P1]dt=k2(Eini2[Ecns2])k2[P1][Ecns2]k3[Ecns1][P1]+k3(Eini1[Ecns1])\frac{d[P_1]}{dt} = k_2(E_\text{ini2}-[E_\text{cns2}])-k_{-2}[P_1][E_\text{cns2}]-k_3[E_\text{cns1}][P_1]+k_{-3}(E_\text{ini1}-[E_\text{cns1}])

d[Ecns1]dt=k3[Ecns1][P1]+(k3+k4)(Eini1[Ecns1])k4[P2][Ecns1]\frac{d[E_\text{cns1}]}{dt} = -k_{3}[E_\text{cns1}][P_{1}] + (k_{-3}+k_4)(E_\text{ini1}-[E_\text{cns1}]) - k_{-4}[P_2][E_\text{cns1}]

d[P2]dt=k4(Eini1[Ecns1])k4[P2][Ecns1]\frac{d[P_2]}{dt} = k_4(E_\text{ini1}-[E_\text{cns1}]) - k_{-4}[P_2][E_\text{cns1}]

function [f] = dXdT(t, x)

k1F = 1;
k1R = 0.1;
k2F = 2;
k2R = 0.1;
k3F = 1;
k3R = 0.1;
k4F = 1;
k4R = 0.1;
E01 = 1;
E02 = 1;

s1=x(1);
e1=x(2);
p1=x(3);
e2=x(4);
p2=x(5);

%ds1dt = 0;
ds1dt = -k1F*e1*s1 + k1R*(E01-e1);
de1dt = -k1F*e1*s1 + (k1R+k2F)*(E01-e1) - k2R*p1*e1;
dp1dt = k2F*(E01-e1) - k2R*p1*e1 - k3F*e2*p1 + k3R*(E02-e2);
de2dt = -k3F*e2*p1 + (k3R+k4F)*(E02-e2) - k4R*p2*e2;
dp2dt = k4F*(E02-e2) - k4R*p2*e2;
f = [ds1dt; de1dt; dp1dt; de2dt; dp2dt];

end

tspan = 0:0.001:10;
initial = [1,1,0,1,0];

[t,x] = ode45( @dXdT, tspan, initial);

Solving the ordinary differential equation system will give us the concentration deviation of the substrate, intermediate product and end product, which is 3'-AMP, 2'-C'-3'-dA and cordycepin respectively. In this simulation, the resource is limited, which means that the substrates will deplete without any supplement.

Figure. 3 Concentration of 3'-AMP, 2'-C'-3'-dA, and cordycepin in fusion expression model with limited resource (1M initial 3'-AMP).

In the figure, the final concentration of the end product, cordycepin, under this set of parameters, is around 0.89. The value of the cordycepin concentration is meaningless, since all the parameters are set optionally. However, it provides a model which can be used in comparison of the separation-expression strategy with the same set of parameters.

Fusion Expression Simulation with abundant resource

The substrate of cns2 is 3'-AMP, while the substrate of cns3-HisG domain is adenosine, which both are the basic substances and liable to produce in a cell. If we provide abundant resource, the concentration of the substrates are considered to be constant.

Therefore, we need to alter the equations of the substrate concentration. We fix the concentration of S1S_1 (3'AMP) to be constant, and let other equations unchanged.

d[S1]dt=0\frac{d[S_1]}{dt} = 0

The result is below:

Figure. 4 Concentration of 3'-AMP, 2'-C'-3'-dA, and cordycepin in fusion expression model with abundant resource (1M 3'-AMP of all times).

From the figure above, we can see that the concentration of the intermediate product 2'-C-3'-dA rises slowly while the end product cordycepin rises much more quickly.

Enzyme Kinetics Simulation of Separation Expression with Diffusion

Considering the separation-expression strategy, we add the diffusion process into the diagram and we can create the following diagram:

Figure. 5 Elucidation of the primary reactions in enzymatic kinetics adding diffusion of P1 from Cns2 side to Cns1 side. S1 represents the substrate of Cns2, 3'-AMP. P1l is the product of Cns2, 2'-C'-3'-dA, before diffusion while P1r is the concentration after diffusion. P2 represents the product of Cns1, cordycepin. Since we consider the reactions in vivo, all the reactions are reversible.

Separation Expression Simulation with limited resource

Adding the term of the diffusion process, we can now list the ordinary differential equations:

d[S1]dt=k1[Ecns2][S1]+k1[Ecns2S1]\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1]

d[Ecns2]dt=k1[Ecns2][S1]+k1[Ecns2S1]+k2[Ecns2S1]k2[P1l][Ecns2]\frac{d[E_\text{cns2}]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1] + k_2[E_\text{cns2}S_1] - k_{-2}[P_{1l}][E_\text{cns2}]

d[Ecns2S1]dt=k1[Ecns2][S1]+k1[Ecns2S1]k2[Ecns2S1]+k2[P1l][Ecns2]\frac{d[E_\text{cns2}S_1]}{dt} = k_1[E_\text{cns2}][S_1] + - k_{-1}[E_\text{cns2}S_1] - k_2[E_\text{cns2}S_1] + k_{-2}[P_{1l}][E_\text{cns2}]

d[P1l]dt=k2[Ecns2S1]k2[P1l][Ecns2]Kdiff[P1l][P1r]exp(dcns1,cns2)\frac{d[P_{1l}]}{dt} = k_2[E_\text{cns2}S_1] - k_{-2}[P_{1l}][E_\text{cns2}] - K_\text{diff}\frac{[P_{1l}]-[P_{1r}]}{exp(d_\text{cns1,cns2})}

d[P1r]dt=Kdiff[P1l][P1r]exp(dcns1,cns2)k3[Ecns1][P1r]+k3[Ecns1S2]\frac{d[P_{1r}]}{dt} = K_\text{diff}\frac{[P_{1l}]-[P_{1r}]}{exp(d_\text{cns1,cns2})} - k_3[E_\text{cns1}][P_{1r}] + k_{-3}[E_\text{cns1}S_2]

d[Ecns1]dt=k3[Ecns1][P1r]+k3[Ecns1S2]+k4[Ecns1S2]k4[P2][Ecns1]\frac{d[E_\text{cns1}]}{dt} = -k_{3}[E_\text{cns1}][P_{1r}] + k_{-3}[E_\text{cns1}S_2] + k_4[E_\text{cns1}S_2] - k_{-4}[P_2][E_\text{cns1}]

d[Ecns1S2]dt=k3[Ecns1][P1r]k3[Ecns1S2]k4[Ecns1S2]+k4[P2][Ecns1]\frac{d[E_\text{cns1}S_2]}{dt} = k_3[E_\text{cns1}][P_{1r}] - k_{-3}[E_\text{cns1}S_2] - k_4[E_\text{cns1}S_2] + k_{-4}[P_2][E_\text{cns1}]

d[P2]dt=k4[Ecns1S2]k4[p2][Ecns1]\frac{d[P_2]}{dt} = k_4[E_\text{cns1}S_2] - k_{-4}[p_2][E_\text{cns1}]

From the equation of d[Ecns2]dt\frac{d[E_\text{cns2}]}{dt} and d[Ecns2S1]dt\frac{d[E_\text{cns2}S_1]}{dt}, we can find that d[Ecns2]dt+d[Ecns2S1]dt=0\frac{d[E_\text{cns2}]}{dt}+\frac{d[E_\text{cns2}S_1]}{dt} = 0[Ecns2]dt+d[Ecns2S1]dt=0{[E_\text{cns2}]}{dt}+\frac{d[E_\text{cns2}S_1]}{dt} = 0, which means that the value of [Ecns2]+[Ecns2S1][E_\text{cns2}]+[E_\text{cns2}S_1] is a constant. Here, we defind this constant as Eini2E_{ini2}, which represents the constant of the initial concentration of cns2 enzyme. Likewise, considering to cns1 enzyme, value of [Ecns1]+[Ecns1S2][E_\text{cns1}]+[E_\text{cns1}S_2] is a constant, which is defined as Eini1E_\text{ini1}.

Therefore, the equations convert to:

d[S1]dt=k1[Ecns2][S1]+k1(Eini2[Ecns2])\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}(E_\text{ini2}-[E_\text{cns2}])

d[Ecns2]dt=k1[Ecns2][S1]+(k1+k2)(Eini2[Ecns2])k2[P1l][Ecns2]\frac{d[E_\text{cns2}]}{dt} = -k_1[E_\text{cns2}][S_1] + (k_{-1}+k_2)(E_\text{ini2}-[E_\text{cns2}]) - k_{-2}[P_{1l}][E_\text{cns2}]

d[P1l]dt=k2(Eini2[Ecns2])k2[P1l][Ecns2]kdiff[P1l][P1r]exp(dcns1,cns2)\frac{d[P_{1l}]}{dt} = k_2(E_\text{ini2}-[E_\text{cns2}]) - k_{-2}[P_{1l}][E_\text{cns2}] - k_\text{diff}\frac{[P_{1l}]-[P_{1r}]}{exp(d_\text{cns1,cns2})}

d[P1r]dt=kdiff[P1l][P1r]exp(dcns1,cns2)k3[Ecns1][P1r]+k3(Eini1[Ecns1])\frac{d[P_{1r}]}{dt} = k_\text{diff}\frac{[P_{1l}]-[P_{1r}]}{exp(d_\text{cns1,cns2})} - k_3[E_\text{cns1}][P_{1r}] + k_{-3}(E_\text{ini1}-[E_\text{cns1}])

d[Ecns1]dt=k3[Ecns1][P1r]+(k3+k4)(Eini1[Ecns1])k4[P2][Ecns1]\frac{d[E_\text{cns1}]}{dt} = -k_{3}[E_\text{cns1}][P_{1r}] + (k_{-3}+k_4)(E_\text{ini1}-[E_\text{cns1}]) - k_{-4}[P_2][E_\text{cns1}]

d[P2]dt=k4(Eini1[Ecns1])k4[P2][Ecns1]\frac{d[P_2]}{dt} = k_4(E_{ini1}-[E_{cns1}]) - k_{-4}[P_2][E_{cns1}]


function [f] = dXdT(t, x)

k1F = 1;
k1R = 0.1;
k2F = 2;
k2R = 0.1;
k3F = 1;
k3R = 0.1;
k4F = 1;
k4R = 0.1;
E1 = 1;
E2 = 1;

s1=x(1);
e1=x(2);
pl=x(3);
pr=x(4);
e2=x(5);
p2=x(6);

kdiff = 10;
d = 1.118/((e1+e2)**(1/3));

ds1dt = -k1F*e1*s1 + k1R*(E1-e1);
de1dt = -k1F*e1*s1 + (k1R+k2F)*(E1-e1) - k2R*pl*e1;
dpldt = k2F*(E1-e1) - k2R*pl*e1 - (kdiff)*(pl-pr)/exp(d);
dprdt = (kdiff)*(pl-pr)/exp(d) - k3F*e2*pr + k3R*(E2-e2);
de2dt = -k3F*e2*pr + (k3R+k4F)*(E2-e2) - k4R*p2*e2;
dp2dt = k4F*(E2-e2) - k4R*p2*e2;
f = [ds1dt; de1dt; dpldt; dprdt; de2dt; dp2dt];

end

tspan = 0:0.001:40;
initial = [1,1,0,0,1,0];

[t,x] = ode45( @dXdT, tspan, initial);

The concentration of the substrate, intermediate products and end product is plotted below:

Figure. 6 Concentration of 3'-AMP, 2'-C'-3'-dA(L), 2'-C'-3'-dA(R), and cordycepin in separate expression model with limited resource (1M initial 3'-AMP). The figure is plotted when kdiff=1.0, k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1

The figure above shows that the concentration of 2'-C'3'-dA is different after diffusion from that before diffusion. The peak of the orange line and grey line is also inconsistent. The peak of the grey line is slightly rightward to the peak of the orange line.

Furthermore, we would like to investigate the influence of the diffusion rate kdiffk_\text{diff} on the concentration of the intermediate products and the end product.


Figure. 7 Concentration of cordycepin in separate expression model with limited resource (1M 3’-AMP of all times). COR_kiff=0.1 means that this line represents the cordycepin concentration when kdiff is 0.1. The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1

In the figure above, 2'-C'3'-dA concentration before diffusion shows a higher level when the diffusion rate is low, and it drops more quickly when the diffusion rate is high, resulting in a low concentration level. Reverse to the 2'-C'3'-dA concentration before diffusion, 2'-C'3'-dA concentration after diffusion presents increasing concentration level when the diffusion rate rises. These figures illustrate how diffusion influence the enzyme kinetics.

The concentration deviation of the end product cordycepin in relation to the increase of kdiffk_\text{diff} suggests that the production efficiency raises when the diffusion rate increases, while the final concentration of the end product is the same value. The red line (cordycepin concentration of fusion-expression strategy) goes much more quickly than all other lines even that of the highest diffusion rate, 10. It is reasonable because fusion-expression strategy without diffusion is considered that kdiffk_{diff} \to \infty.

Separation Expression Simulation with abundant resource

Just as the manipulation of the fusion expression model, we set the concentration of the substrate S1S_1 to a constant. The result is below:

Figure. 8 Concentration of 3'-AMP, 2'-C'-3'-dA(L), 2'-C'-3'-dA(R), and cordycepin in separate expression model with abundant resource (1M 3'-AMP of all times). The figure is plotted when kdiff=1.0, k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1

Different from the case when the resource is limited, all the products (including the intermediate products and end product) are increasing without any dropping tendency. The difference between 2'-C-3'-dA before and after diffusion becomes constantly bigger. After a while, the production of cordycepin is nearly a straight line.

Figure. 9 Concentration of cordycepin in separate expression model with abundant resource (1M 3'-AMP of all times). The figure is plotted k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1

Considering the influence of kdiffk_\text{diff} on the production of cordycepin and comparing the production of cordycepin under fusion-expression strategy and separation-expression strategy, we can derive the conclusion that the higher kdiffk_\text{diff} is the higher production efficiency is.

Summary

The result of the separation expression model, both in limited resource and abundant resource, testifies that fusion-expression strategy, compared with separation-expression strategy, has a more efficient production of cordycepin.

Based on the model and the result, we prove the priority of fusion expression in theory by enzyme kinetic. This good result gives us strong confidence and support on the experiments of fusion expression of Cns1 and Cns2.

Systematic Enzyme Kinetics

In the fusion expression simulation, we discussed the catalytic process of Cns1 and Cns2, leaving alone the function of Cns3 and its product pentostatin, which is an inhibitor to the adenosine deaminase (ADA). When we add the role of cns3-HisG domain and ADA, the systematic enzyme kinetics is represented by the following diagram.

Figure. 10 Elucidation of the primary reactions in enzymatic kinetics. S1 represents the substrate of Cns2, 3'-AMP. P1 is the product of Cns2, 2'-C'-3'-dA, and it is the substrate of Cns1. P2 represents the product of Cns1, cordycepin. I represents inhibitor, pentostatin. Since we consider the reactions in vivo, all the reactions are reversible.

In the diagram, product of Cns1, P2P_2 (cordycepin), is deaminized to 3'-dI (P3P_3 ) by adenosine deaminase (ADA), which is the EADAE_{ADA} in the diagram. However, inhibitor (pentostatin), which is the II in the diagram, can competitively inhibit the catalytic activity of ADA (Wu P et al., 2017. Li G et al., 2015). Pentostatin and ADA will from a complex, which is not functional and can not produce any substances.

Systematic Simulation with Limited Resource

As the cases in the above sections, we can derive the kinetic equations from the diagram.

d[S1]dt=k1[Ecns2][S1]+k1[Ecns2S1]\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1]

d[Ecns2]dt=k1[Ecns2][S1]+k1[Ecns2S1]+k2[Ecns2S1]k2[P1][Ecns2]\frac{d[E_\text{cns2}]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1] + k_2[E_\text{cns2}S_1] - k_{-2}[P_{1}][E_\text{cns2}]

d[Ecns2S1]dt=k1[Ecns2][S1]k1[Ecns2S1]k2[Ecns2S1]+k2[P1][Ecns2]\frac{d[E_\text{cns2}S_1]}{dt} = k_1[E_\text{cns2}][S_1] - k_{-1}[E_\text{cns2}S_1] - k_2[E_\text{cns2}S_1] + k_{-2}[P_{1}][E_\text{cns2}]

d[P1]dt=k2[Ecns2S1]k2[P1][Ecns2]k3[Ecns1][P1]+k3[Ecns1S2]\frac{d[P_1]}{dt} = k_2[E_\text{cns2}S1]-k_{-2}[P_1][E_\text{cns2}]-k_3[E_\text{cns1}][P_1]+k_{-3}[E_\text{cns1}S_2]

d[Ecns1]dt=k3[Ecns1][P1]+k3[Ecns1S2]+k4[Ecns1S2]k4[P2][Ecns1]\frac{d[E_\text{cns1}]}{dt} = -k_{3}[E_\text{cns1}][P_{1}] + k_{-3}[E_\text{cns1}S_2] + k_4[E_\text{cns1}S_2] - k_{-4}[P_2][E_\text{cns1}]

d[Ecns1S2]dt=k3[Ecns1][P1]k3[Ecns1S2]k4[Ecns1S2]+k4[P2][Ecns1]\frac{d[E_\text{cns1}S_2]}{dt} = k_3[E_\text{cns1}][P_{1}] - k_{-3}[E_\text{cns1}S_2] - k_4[E_\text{cns1}S_2] + k_{-4}[P_2][E_\text{cns1}]

d[P2]dt=k4[Ecns1S2]k4[P2][Ecns1]k5[EADA][P2]+k5[EADAS3]\frac{d[P_2]}{dt} = k_4[E_\text{cns1}S_2] - k_{-4}[P_2][E_\text{cns1}] - k_5[E_\text{ADA}][P_2] + k_{-5}[E_\text{ADA}S_3]

d[EADA]dt=k5[EADA][P2]+k5[EADAS3]+k6[EADAS3]k6[P3][EADA]k9[I][EADA]+k9[EADAI]\frac{d[E_\text{ADA}]}{dt} = -k_5[E_\text{ADA}][P_2] + k_{-5}[E_\text{ADA}S_3] + k_6[E_\text{ADA}S_3] - k_{-6}[P_3][E_\text{ADA}] - k_9[I][E_\text{ADA}] + k_{-9}[E_\text{ADA}I]

d[EADAS3]dt=K5[EADA][P2]k5[EADAS3]k6[EADAS3]+k6[P3][EADA]\frac{d[E_\text{ADA}S_3]}{dt} = K_5[E_\text{ADA}][P_2] - k_{-5}[E_\text{ADA}S_3] - k_6[E_\text{ADA}S_3] + k_{-6}[P_3][E_\text{ADA}]

d[P3]dt=k6[EADAS3]k6[P3][EADA]\frac{d[P_3]}{dt} = k_6[E_\text{ADA}S_3] - k_{-6}[P_3][E_\text{ADA}]

d[S4]dt=k7[EHisG][S4]+k7[EHisGS4]\frac{d[S_4]}{dt} = -k_7[E_\text{HisG}][S_4] + k_{-7}[E_\text{HisG}S_4]

d[EHisG]dt=k7[EHisG][S4]+k7[EHisGS4]+k8[EHisGS4]k8[EHisG][S4]\frac{d[E_\text{HisG}]}{dt} = -k_7[E_\text{HisG}][S_4] + k_{-7}[E_\text{HisG}S_4] + k_8[E_\text{HisG}S_4] - k_{-8}[E_\text{HisG}][S_4]

d[EHisGS4]dt=k7[EHisG][S4]k7[EHisGS4]k8[EHisGS4]+k8[EHisG][I]\frac{d[E_\text{HisG}S_4]}{dt} = k_7[E_\text{HisG}][S_4] - k_{-7}[E_\text{HisG}S_4] - k_8[E_\text{HisG}S_4] + k_{-8}[E_\text{HisG}][I]

d[I]dt=k8[EHisGS4]k8[EHisG][S4]k9[I][EADA]+k9[EADAI]\frac{d[I]}{dt} = k_8[E_\text{HisG}S_4] - k_{-8}[E_\text{HisG}][S_4] - k_9[I][E_\text{ADA}] + k_{-9}[E_\text{ADA}I]

d[EADAI]dt=k9[I][EADA]k9[EADAI]\frac{d[E_\text{ADA}I]}{dt} = k_9[I][E_\text{ADA}] - k_{-9}[E_\text{ADA}I]

Likewise, we can simplify this 14-equation system to 10-equation system because the initial concentrations of the enzymes are fixed. In the enzymes, EADAE_\text{ADA} is a little different, since the case that [EADA]+[EADAS3]+[EI]=Eini3[E_\text{ADA}] + [E_\text{ADA}S_3] + [EI] = E_\text{ini3}, so [EADAS3][E_\text{ADA}S_3] can be represented by Eini3[EADA][EI]E_\text{ini3} - [E_\text{ADA}] - [EI]

.

Thus, the equations can be simplified to:

d[S1]dt=k1[Ecns2][S1]+k1(Eini2[Ecns2])\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}(E_\text{ini2}-[E_\text{cns2}])

d[Ecns2]dt=k1[Ecns2][S1]+(k1+k2)(Eini2[Ecns2])k2[P1][Ecns2]\frac{d[E_\text{cns2}]}{dt} = -k_1[E_\text{cns2}][S_1] + (k_{-1}+k_2)(E_\text{ini2}-[E_\text{cns2}]) - k_{-2}[P_{1}][E_\text{cns2}]

d[P1]dt=k2(Eini2[Ecns2])k2[P1][Ecns2]k3[Ecns1][P1]+k3(Eini1[Ecns1])\frac{d[P_1]}{dt} = k_2(E_\text{ini2}-[E_\text{cns2}])-k_{-2}[P_1][E_\text{cns2}]-k_3[E_\text{cns1}][P_1]+k_{-3}(E_\text{ini1}-[E_\text{cns1}])

d[Ecns1]dt=k3[Ecns1][P1]+(k3+k4)(Eini1[Ecns1])k4[P2][Ecns1]\frac{d[E_\text{cns1}]}{dt} = -k_{3}[E_\text{cns1}][P_{1}] + (k_{-3}+k_4)(E_\text{ini1}-[E_\text{cns1}]) - k_{-4}[P_2][E_\text{cns1}]

d[P2]dt=k4(Eini1[Ecns1])k4[P2][Ecns1]k5[EADA][P2]+k5(Eini3[EADA][EI])\frac{d[P_2]}{dt} = k_4(E_\text{ini1}-[E_\text{cns1}]) - k_{-4}[P_2][E_\text{cns1}] - k_5[E_\text{ADA}][P_2] + k_{-5}(E_\text{ini3} - [E_\text{ADA}] - [EI])

d[EADA]dt=k5[EADA][P2]+(k5+k6)(Eini3[EADA][EI])k6[P3][EADA]k9[I][EADA]+k9[EADAI]\frac{d[E_\text{ADA}]}{dt} = -k_5[E_\text{ADA}][P_2] + (k_{-5}+k_6)(E_\text{ini3} - [E_\text{ADA}] - [EI]) - k_{-6}[P_3][E_\text{ADA}] - k_9[I][E_\text{ADA}] + k_{-9}[E_\text{ADA}I]

d[P3]dt=k6(Eini3[EADA][EI])k6[P3][EADA]\frac{d[P_3]}{dt} = k_6(E_\text{ini3} - [E_\text{ADA}] - [EI]) - k_{-6}[P_3][E_\text{ADA}]

d[S4]dt=k7[EHisG][S4]+k7(EiniH[EHisG])\frac{d[S_4]}{dt} = -k_7[E_\text{HisG}][S_4] + k_{-7}(E_\text{iniH}-[E_\text{HisG}])

d[EHisG]dt=k7[EHisG][S4]+(k7+k8)(EiniH[EHisG])k8[EHisG][S4]\frac{d[E_\text{HisG}]}{dt} = -k_7[E_\text{HisG}][S_4] + (k_{-7}+k_8)(E_\text{iniH}-[E_\text{HisG}]) - k_{-8}[E_\text{HisG}][S_4]

d[I]dt=k8(EiniH[EHisG])k8[EHisG][S4]k9[I][EADA]+k9[EADAI]\frac{d[I]}{dt} = k_8(E_\text{iniH}-[E_\text{HisG}]) - k_{-8}[E_\text{HisG}][S_4] - k_9[I][E_\text{ADA}] + k_{-9}[E_\text{ADA}I]

d[EADAI]dt=k9[I][EADA]k9[EADAI]\frac{d[E_\text{ADA}I]}{dt} = k_9[I][E_\text{ADA}] - k_{-9}[E_\text{ADA}I]

function [f] = dXdT(t, x)

k1F = 1;
k1R = 0.1;
k2F = 2;
k2R = 0.1;
k3F = 1;
k3R = 0.1;
k4F = 1;
k4R = 0.1;
k5F = 1;
k5R = 0.5;
k6F = 1;
k6R = 0.5;
k7F = 1;
k7R = 0.1;
k8F = 1;
k8R = 0.1;
k9F = 10;
k9R = 0.1;
E01 = 1;
E02 = 1;
E03 = 1;
E04 = 1;

s1=x(1);
e1=x(2);
p1=x(3);
e2=x(4);
p2=x(5);
e3=x(6);
p3=x(7);
s4=x(8);
e4=x(9);
i=x(10);
ei=x(11);

ds1dt = -k1F*e1*s1 + k1R*(E01-e1);
%dsdt = 0;
de1dt = -k1F*e1*s1 + (k1R+k2F)*(E01-e1) - k2R*p1*e1;
dp1dt = k2F*(E01-e1) - k2R*p1*e1 - k3F*e2*p1 + k3R*(E02-e2);
de2dt = -k3F*e2*p1 + (k3R+k4F)*(E02-e2) - k4R*p2*e2;
dp2dt = k4F*(E02-e2) - k4R*p2*e2 - k5F*e3*p2 + k5R*(E03-e3-ei);
de3dt = -k5F*e3*p2 + (k5R+k6F)*(E03-e3-ei) - k6R*p3*e3 - k9F*i*e3 + k9R*ei;
dp3dt = k6F*(E03-e3-ei) - k6R*p3*e3;
ds4dt = -k7F*e4*s4 + k7R*(E04-e4);
%ds4dt = 0;
de4dt = -k7F*e4*s4 + (k7R+k8F)*(E04-e4) - k8R*e4*s4;
didt  = k8F*(E04-e4) - k8R*e4*s4 - k9F*i*e3 + k9R*ei;
deidt = k9F*i*e3 - k9R*ei;
f = [ds1dt; de1dt; dp1dt; de2dt; dp2dt; de3dt; dp3dt; ds4dt; de4dt; didt; deidt];

end

tspan = 0:0.001:100;
initial = [1,1,0,1,0,1,0,2,1,0,0];

[t,x] = ode45( @dXdT, tspan, initial);

The result shows the concentration deviation of the substances involved.

Figure. 11 Concentration of 3’-AMP, cordycepin (COR), adenosine and Pentostatin (PTN) in systematic model with limited resource (1M 3’-AMP, 2M 3’-AMP and 3M 3’-AMP respectively). The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1, k5 = 1, k-5 = 0.5, k6F = 1, k-6 = 0.5, k7 = 1, k-7 = 0.1, k8 = 1, k-8 = 0.1, k9 = 10, k-9 = 0.1;

Comparison of the initial concentration of 3'-AMP (1, 2, 3 respectively) shows that when the resource is limited, equal concentration of cordycepin pentostatin will cause a fair part of cordycepin converting to nonfunctional product, 3'-dI. However, substrate ratio above 2:1 (3'-AMP: adenosine) will dramatically repress this conversion and enable more production of cordycepin.

After the overall description of the model, we would like to evaluate the parameter of k9k_9, which is the reaction constant of the combination reaction of inhibitor and ADA, which is a kind of indicator for measuring the inhibition strength of pentostatin. From the figure, we can see that as the increase of k9k_9, the max COR production as peak goes up and the decrease is inhibited.

Figure. 12 Concentration of cordycepin and pentostatin in systematic model with limited resource. Parameter k9 varies from 1 to 20. The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1, k5 = 1, k-5 = 0.5, k6F = 1, k-6 = 0.5, k7 = 1, k-7 = 0.1, k8 = 1, k-8 = 0.1, k-9 = 0.1;

The influence of k9k_9 on PTN production is shown above, as the lines increase slower in initial stage when k9k_9 is higher. However, the production level of PTN is nearly stationary.

The result gives us a conclusion that the best harvest time of cordycepin is slightly after the peak and the higher k9k_9 is, the slower conversion of cordycepin to nonfunctional product.

Systematic simulation with abundant resource

The substrates in the system are S1S_1 (3'-AMP) and S4S_4 (adenosine), which would be constant with abundant resource.

d[S1]dt=0\frac{d[S_1]}{dt} = 0

d[S4]dt=0\frac{d[S_4]}{dt} = 0

Leaving other equations unchanged and we could get the result below.

Figure. 13 Concentration of 3’-AMP, COR, 3’-dI, adenosine, and PTN in systematic model with abundant resource (1M 3’-AMP and 1M adenosine of all times). The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1, k5 = 1, k-5 = 0.5, k6F = 1, k-6 = 0.5, k7 = 1, k-7 = 0.1, k8 = 1, k-8 = 0.1, k9 = 10, k-9 = 0.1;

The result suggests that when the substrates are abundant, the concentration of pentostatin will raise more quickly as time goes because the concentration of the inhibition complex (EADAIE_\text{ADA}I) amasses and it needs less inhibitor (pentostatin) to form the complex.

Additionally, we need to discuss the optional sustitute ratio ([3’-AMP][adenosine]\frac{[\text{3'-AMP}]}{[\text{adenosine}]}) based on current parameters which enable the max production of cordycepin.

Figure. 14 Concentration of COR in systematic model with abundant resource. Initial concentration of 3’-AMP is 1M and initial concentration of adenosine varies from 0.1M to 3M (concentrations of 3’-AMP and adenosine are constant). The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1, k5 = 1, k-5 = 0.5, k6F = 1, k-6 = 0.5, k7 = 1, k-7 = 0.1, k8 = 1, k-8 = 0.1, k9 = 10, k-9 = 0.1;

When we fixed the 3'-AMP concentration to 1M and the adenosine concentration to a set of values, 0.1, 0.2, 0.5, 1.0, 2.0 and 3.0, we plotted the figure above. In the figure, as the increase of initial adenosine concentration, the producing rate of cordycepin increases. The mechanism behind the result is that more pentostatin is produced with higher adenosine concentration and the conversion of cordycepin to 3'-dI is inhibited, resulting in more cordycepin remained.

The result above is reasonable, but meaningless, because we all know that more production of pentostatin will create more cordycepin. But what if we set a bound for total resource? In nature, the total energy of a cell at a time is limited although the energy can get supplemented by absorbing external resources. It is true that if a cell spends more energy in the production of cordycepin, it will spend less energy in pentostatin. Therefore, in the model, we set the sum of the concentration of 3'-AMP and adenosine to a fixed value and the two concentrations are stationary respectively. If we change the ratio between 3'-AMP and adenosine, the two concentrations will vary but the whole concentration is constant. Based on the manipulation above, we solve the equation system and get the following result:

Figure. 15 Concentration of COR in systematic model with abundant resource. Sum of 3’-AMP and adenosine concentration is constant, 2M. Initial concentration of 3’-AMP varies from 0.5M to 1M (adenosine otherwise from 1.5 M to 0). The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1, k5 = 1, k-5 = 0.5, k6F = 1, k-6 = 0.5, k7 = 1, k-7 = 0.1, k8 = 1, k-8 = 0.1, k9 = 10, k-9 = 0.1;

From the figure above, we can find that when the ratio between 3'-AMP and adenosine is below a threshold, a high ratio will result in more cordycepin production efficiency. However, the case is inconsistent when the ratio exceeds the threshold. As the lines of the colors of light blue, green and blown, the higher ratio will generate a lower production rate. In the extreme case, when adenosine is absent, the production of cordycepin will decrease to the level when the ratio is 0.5/1.5.

For better visualization, if we use the initial substrates (3'-AMP and adenosine) ratio to represent the energy distribution on cordycepin and pentostatin, we can identify the optimal utilization of resources for cordycepin production. We measured the concentration of cordycepin 20 hours after the begining point with various cordycepin synthesis ratios (persentage of resources or energy used in the production of cordycepin).


function [f] = dXdT(t, x)

k1F = 1;
k1R = 0.1;
k2F = 2;
k2R = 0.1;
k3F = 1;
k3R = 0.1;
k4F = 1;
k4R = 0.1;
k5F = 1;
k5R = 0.5;
k6F = 1;
k6R = 0.5;
k7F = 1;
k7R = 0.1;
k8F = 1;
k8R = 0.1;
k9F = 10;
k9R = 0.1;
E01 = 1;
E02 = 1;
E03 = 1;
E04 = 1;

s1=x(1);
e1=x(2);
p1=x(3);
e2=x(4);
p2=x(5);
e3=x(6);
p3=x(7);
s4=x(8);
e4=x(9);
i=x(10);
ei=x(11);

%ds1dt = -k1F*e1*s1 + k1R*(E01-e1);
ds1dt = 0;
de1dt = -k1F*e1*s1 + (k1R+k2F)*(E01-e1) - k2R*p1*e1;
dp1dt = k2F*(E01-e1) - k2R*p1*e1 - k3F*e2*p1 + k3R*(E02-e2);
de2dt = -k3F*e2*p1 + (k3R+k4F)*(E02-e2) - k4R*p2*e2;
dp2dt = k4F*(E02-e2) - k4R*p2*e2 - k5F*e3*p2 + k5R*(E03-e3-ei);
de3dt = -k5F*e3*p2 + (k5R+k6F)*(E03-e3-ei) - k6R*p3*e3 - k9F*i*e3 + k9R*ei;
dp3dt = k6F*(E03-e3-ei) - k6R*p3*e3;
%ds4dt = -k7F*e4*s4 + k7R*(E04-e4);
ds4dt = 0;
de4dt = -k7F*e4*s4 + (k7R+k8F)*(E04-e4) - k8R*e4*s4;
didt  = k8F*(E04-e4) - k8R*e4*s4 - k9F*i*e3 + k9R*ei;
deidt = k9F*i*e3 - k9R*ei;
f = [ds1dt; de1dt; dp1dt; de2dt; dp2dt; de3dt; dp3dt; ds4dt; de4dt; didt; deidt];

end

function f = FProduct(a, b) %a is 3'-AMP, b is adenosine

tspan = 0:0.001:20;
initial = [a,1,0,1,0,1,0,b,1,0,0];

[t,x] = ode45( @dXdT, tspan, initial);
x(end,5)

end

for a = 0:0.01:1
b = 1-a;
[a,FProduct(a,b)]

end

Figure. 16 Cordycepin (COR) concentration in systematic model with abundant resource. Sum of 3’-AMP and adenosine concentration is constant, 1M. Initial concentration of 3’-AMP varies from 0 to 1M (adenosine otherwise from 1 M to 0). Cordycepin production level percentage is evaluated by the initial 3’-AMP percentage. The figure is plotted when k1 = 1, k-1 = 0.1, k2 = 2, k-2 = 0.1, k3 = 1, k-3 = 0.1, k4 = 1, k-4 = 0.1, k5 = 1, k-5 = 0.5, k6F = 1, k-6 = 0.5, k7 = 1, k-7 = 0.1, k8 = 1, k-8 = 0.1, k9 = 10, k-9 = 0.1;

We performed enzyme kinetic simulation to investigate the relationship between the synthesis of cordycepin and pentostatin and the potential methods to maximize cordycepin production. The result above gives us a clear illustration in theory that:

The result above gives us a clear illustration in theory that:

  • Pentostatin does help promote the production of cordycepin.
  • Appropriate ratio between the catalytic rate of cordycepin and pentostatin will result in the most productive level of cordycepin. The ratio can be adjusted by change of enzyme expression level if we have no access to changing the substrate level.

  • Since the goal of our project is to increase the efficiency and yield of cordycepin production, this result is helpful to the manufacturing. Although the parameters are set optional owing to lack of information, it will give a brief guideline of the optional ratio of cordycepin and pentostatin and will benefit the factories a lot when the enzymatic parameter is provided.

    Appendix

    Symbol Description
    [ ] The concentration symbol
    S1S_1 substitue adenosine
    Ecns1E_{cns1} enzyme Cns1
    Ecns2E_{cns2} enzyme Cns2
    Ecns2S1E_{cns2}S_1 enzyme-substrate complex of cns2 and S1S_1
    Ecns1S2E_{cns1}S_2 enzyme-substrate complex of Cns1 and S2S_2
    P1EnolP_{1Enol} product of enzyme Cns1, which is in its enol form
    P1KetolP_{1Keto-l} keto form of the product of Cns2, at location of enzyme Cns1
    P1KetorP_{1Keto-r} keto form of the product of Cns2, at location of enzyme Cns2, after diffusion
    P2P_2 product of enzyme Cns1
    EADAE_{ADA} enzyme deoxyadenosine (ADA)
    EADAS3E_{ADA}S_3 enzyme-substrate complex of ADA and S3(P2)S_3(P_2)
    P3P_3 product of enzyme ADA by deamination of P2P_2 (cordycepin)
    S4S_4 substrate adenosine
    EHisGE_{HisG} HisG domain of enzyme cns3, catalyzing convertion of adenosine to pentostatin
    EHisGS4E_{HisG}S_4 enzyme-substrate complex of HisG and S4S_4
    II product of enzyme HisG, pentostatin
    Parameter Description Value Reference
    k1k_1 forward rate constant of Cns2-S1 complex formation 1
    k1k_{-1} Reverse rate constant of Cns2-S1 complex formation 0.1
    k2k_2 forward rate constant of P1EnolP_{1Enol}/P1P_1 product production 2
    k2k_{-2} reverse rate constant of P1EnolP_{1Enol}/P1P_1 product production 0.1
    k3k_3 forward rate constant of Cns1-S2 complex formation 1
    k3k_{-3} Reverse rate constant of Cns1-S2 complex formation 0.1
    k4k_4 forward rate constant of P2P_2 product production 1
    k4k_{-4} reverse rate constant of P2P_2 product production 0.1
    k5k_5 forward rate constant of EADAS3E_{ADA}S_3 complex formation 1
    k5k_{-5} Reverse rate constant of EADAS3E_{ADA}S_3 complex formation 0.5
    k6k_6 forward rate constant of P3P_{3} product production 1
    k6k_{-6} reverse rate constant of P3P_{3} product production 0.5
    k7k_7 forward rate constant of HisG-S4 complex formation 3
    k7k_{-7} Reverse rate constant of HisG-S4 complex formation 0.1
    k8k_8 forward rate constant of II product production 3
    k8k_{-8} reverse rate constant of II product production 0.1
    k9k_9 forward rate constant of EADAIE_{ADA}I product production 10
    k9k_{-9} reverse rate constant of EADAIE_{ADA}I product production 0.1

    Reference

    Xia, Y. , Luo, F. , Shang, Y. , Chen, P. , Lu, Y. , & Wang, C. . (2017). Fungal cordycepin biosynthesis is coupled with the production of the safeguard molecule pentostatin. Cell Chemical Biology, S2451945617303276.

    Kozanecki, M. , Halagan, K. , Saramak, J. , & Matyjaszewski, K. . (2016). Diffusive properties of solvent molecules in neighborhood of polymer chain as seen by monte-carlo simulations. Soft Matter, 12(25), 5519.

    Huang, L. , Huang, Y. , Chi, Y. , Lin, J. , Yu, L. , & Xu, L. , et al. (2007). Study on the kinetics of keto-enol tautomerism of p-hydroxyphenylpyruvic acid using capillary electrophoresis. Journal of Chromatography A, 1175(2), 283-288.

    Eun, C. , Kekenes-Huskey, P. M. , Metzger, V. T. , & Mccammon, J. A. . (2014). A model study of sequential enzyme reactions and electrostatic channeling. The Journal of Chemical Physics, 140(10), 105101.

    Rice, S. A. . (1985). Diffusion-limited reactions. Comprehensive Chemical Kinetics.

    Macheras, P. , & Iliadis, A. . (2006). Diffusion and Kinetics. Modeling in Biopharmaceutics, Pharmacokinetics, and Pharmacodynamics. Springer New York.

    M.H. Jacobs, Diffusion Processes, in: M.H. Jacobs (Ed.), Diffusion Processes, Springer Berlin Heidelberg, Berlin, Heidelberg, 1935, pp. 1-145.

    Wu, P., Wan, D., Xu, G., Wang, G., Ma, H., Wang, T., ... & Li, Y. Q. (2017). An unusual protector-protégé strategy for the biosynthesis of purine nucleoside antibiotics. Cell chemical biology, 24(2), 171-181.

    Li, G., Nakagome, I., Hirono, S., Itoh, T., & Fujiwara, R. (2015). Inhibition of adenosine deaminase (ADA)‐mediated metabolism of cordycepin by natural substances. Pharmacology research & perspectives, 3(2), e00121.