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 substitute 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 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 to segregate 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 substitute of the first enzyme is sequentially converted to the product of the second enzyme, which is more efficient than 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:

$\frac{dC}{dt} = -k_{eq}C$

in which $C$ 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 ($k_{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 substitute 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 substitute 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 substitute 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×10^{23}$ molecules/L, = 0.6 molecules/$\text{nm}^3$, or the volume per molecule is V=1.66 $nm^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 = 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 $R_\text{min}$, because it is the minimal radius of a sphere that could contain the given mass of protein:

$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 $\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 $C_\text{cns1}$ , and concentration of Cns2 to be $C_\text{cns2}$ , there are totally $(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:

$d_{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 substitute 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) = -D \frac{\partial c(z,t)}{\partial z}$

• $J$ is the diffusion flux, of which the dimension is amount of substance per unit area per unit time. $J$ measures the amount of substance that will flow through a unit area during a unit time interval.
• $D$ is the diffusion coefficient or diffusivity. Its dimension is area per unit time.
• $c$ is the concentration, of which the dimension is amount of substance per unit volume.
• $z$ 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 \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:

$\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) = c_l(t)$ and $c(\delta,0) = c_r(t)$, the differential ordinary equation can be solved as the form below. (M.H. Jacobs., 1935)

$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, $c_r(t)$ is the substitute concentration of the second enzyme around the second enzyme, while $c_l(t)$ is the substitute 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 $z$ to $\delta$, we can get:

Since $sin(i\pi) = 0$,
$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\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=\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 $x$, we can get:

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

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

$\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 $x$. The coefficient $k_\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 substitute 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 substitute of Cns2, 3’-AMP. P1 is the product of Cns2, 2’-C’3’-dA, and it is the substitute 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:

$\frac{d[S_1]}{dt} = -k_1[E_\text{cns2}][S_1] + k_{-1}[E_\text{cns2}S_1]$

$\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}]$

$\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}]$

$\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]$

$\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}]$

$\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}]$

$\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 $\frac{d[E_\text{cns2}]}{dt}$ and $\frac{d[E_\text{cns2}S_1]}{dt}$, we can find that $\frac{d[E_\text{cns2}]}{dt}+\frac{d[E_\text{cns2}S_1]}{dt} = 0$, which means that the value of $[E_\text{cns2}]+[E_\text{cns2}S_1]$ is a constant. Here, we defined this constant as $E_\text{ini2}$, which represents the constant of the initial concentration of cns2 enzyme. Likewise, considering to cns1 enzyme, value of $[E_\text{cns1}]+[E_\text{cns1}S_2]$ is a constant, which is defined as $E_\text{ini1}$.

Therefore, the equations convert to:

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

$\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}]$

$\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}])$

$\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}]$

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

MATLAB Code:

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 substitute, 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 substitutes 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.

1.3.2 Fusion Expression Simulation with abundant resource

The substitute of cns2 is 3'-AMP, while the substitute 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 substitutes are considered to be constant.

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

$\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.

1.4 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 substitute 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.

1.4.1 Separation Expression Simulation with limited resource

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

$\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}]$

$\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}]$

$\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})}$

$\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]$

$\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}]$

$\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}]$

$\frac{d[P_2]}{dt} = k_4[E_\text{cns1}S_2] - k_{-4}[p_2][E_\text{cns1}]$

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

Therefore, the equations convert to:

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

$\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}]$

$\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})}$

$\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}])$

$\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}]$

$\frac{d[P_2]}{dt} = k_4(E_{ini1}-[E_{cns1}]) - k_{-4}[P_2][E_{cns1}]$