\( \tau-leaping \) model

In the reaction system which produces CDCA-FXR-RXR, the number of reactive molecules involved is relatively small and results in obvious randomness. So we tried to find something different from the rate equation for these uncertain processes.

Gillespie algorithm can accurately simulate the long-term behavior of the system, but it is inefficient because it assumes that no reaction will take place in the time interval \([t,t+\tau]\) and one reaction will take place in the infinitesimal time interval \([t+\tau,t+\tau+d\tau]\), which makes it necessary to simulate every step.

Considering the two points above, in order to simulate the system behavior more accurately and improve the simulation speed, we chose "\(\tau-leaping\) algorithm" improved from Gillespie algorithm.

Establishment of model

When a system has M reaction channels \(\left\{R_1, \ldots \ldots, R_M\right\}\) and N substances \(\left\{S_1, \ldots \ldots, S_N\right\}\) , the state vector \(X(t)=\left(X_1(t), \ldots \ldots, X_N(t)\right)^T\) can represent its dynamic state at time t, and \(X_i(t)(i=1, \ldots, N)\) represents the number of substances \(S_{i}\) at time t. \(a_{j}(x)\) represents the tendency function of the reaction channel \(R_{j}\) at time t, \(V_{j}=\left(V_{1 j}, \ldots, V_{N j}\right)^T\) represents the state change vector, where \(V_{i j}\) represents molecular changes number of substance \(S_{i}\) when \(R_{j}\) occurs.

For reaction \(R_{1}\):

\( D N A_{FXR} \stackrel{T c_{\text {FXR}}}{\longrightarrow} m R N A_{FXR} \)

Its state change vector is:

\(V_{1}=\left[\begin{array}{lllllllll} 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right]\)

and tendency function

\(a_{1}=N_{\text {plasmid }} \times T c_{FXR}\)

\(\tau\) must firstly satisfy the jumping condition and should be as large as possible to speed up. Choose tau as:

\( \tau = \underset{j \in[1, M]}{\min}\left\{\frac{\epsilon a_{0}(x)}{\left|\sum_{i=1}^{N} \xi_{i}(x) b_{j i}(x)\right|}\right\} \)

Where \(\xi(x)=\sum_{j=1}^{M} a_{j}(x) V_{j}, a_{0}=\sum_{j=1}^{M} a_{j}(x), b_{j i}(x)=\partial a_{j}(x) / \partial x_{i}\)

Then M independent Poisson distribution random numbers are used to approximate the reaction times \(K_{j}=P_{j}\left\{a_{j}(x), \tau\right\}\) of all reaction channels in the time interval.

Then the increasing amount of system state \(\lambda=\sum_{j=1}^{M} K_{j} v_{j}\) is defined, so the dynamic process of the system can always be expressed as:

\( X_{i}(t+\tau)=X_{i}+\sum_{j=1}^{M} v_{j i} P_{j}\left(a_{j}(x), \tau\right)(i=1, \ldots, N) \)

Table 1. Reactions

Reaction channel
\(R_{1}\)	\(D N A_{FXR} \stackrel{T c_{\text {FXR}}}{\longrightarrow} m R N A_{FXR}\)
\(R_{2}\)	\(D N A_{RXR} \stackrel{T c_{\text {RXR}}}{\longrightarrow} m R N A_{RXR}\)
\(R_{3}\)	\(mR N A_{FXR} \stackrel{T s_{\text {FXR}}}{\longrightarrow} FXR\)
\(R_{4}\)	\(mR N A_{RXR} \stackrel{T s_{\text {RXR}}}{\longrightarrow} RXR\)
\(R_{5}\)	\(m R N A_{FXR} \stackrel{D m_{FXR}}{\longrightarrow} \emptyset\)
\(R_{6}\)	\(m R N A_{RXR} \stackrel{D m_{RXR}}{\longrightarrow} \emptyset\)
\(R_{7}\)	\(FXR \stackrel{D s_{FXR}}{\longrightarrow} \emptyset\)
\(R_{8}\)	\(RXR \stackrel{D s_{RXR}}{\longrightarrow} \emptyset\)
\(R_{9}\)	\(CDCA + FXR \stackrel{k_{1}}{\longrightarrow} CDCA-FXR\)
\(R_{10}\)	\(CDCA-FXR + RXR \stackrel{k_{2}}{\longrightarrow} CDCA-FXR-RXR\)
\(R_{11}\)	\(ddRFP-a_{1} + ddRFP-b_{1} \stackrel{k_{3}}{\longrightarrow} ddRFP-a_{1}-b_{1}\)

\(R_{1},\, R_{2},\, R_{3},\, R_{4}\, \)represents the transcription and translation processes involved in the synthesis of protein FXR and protein RXR.

\(R_{5},\, R_{6},\, R_{7},\, R_{8}\, \)represents the degradation process of four species(\(mRNA_{FXR}, mRNA_{RXR}, FXR, RXR\)).

\(R_{9}\, \)represents the process of CDCA combining with FXR to become CDCA-FXR.

\(R_{10}\, \)represents the process of CDCA-FXR combining with RXR to become CDCA-FXR-RXR.

\(R_{11}\, \)represents the bottom noise reaction of the combination of two ddRFP.

Solution of model

State vector: According to the reaction equation above, the substances involved in our model can be expressed as state vectors:

\( \mathrm{X}(\mathrm{t})=\left(\mathrm{DNA}_{\text {FXR }}, \mathrm{DNA}_{\text {RXR }}, \mathrm{mRNA}_{FXR}, \mathrm{mRNA}_{\text {RXR}}\right. ,\\ FXR, RXR, \)

\(CDCA, CDCA-FXR, CDCA-FXR-RXR,ddRFP-a_{1}-b_{1 }) \)

State change vectors:

Table 2. Tendency functions

The reaction channel	Tendency functions	Explanation
\(R_{1} \& R_{2}\)	\(F_{1}=N_{plasmid}×Tc_{FXR} \) \( F_{2}=N_{plasmid}×Tc_{RXR}\)	The pET vector existed as a plasmid of low copy number in host E. coli, which is roughly distributed in the interval [15-20], so we assumed that the average plasmid copy number is a constant 20.
\(R_{3} \& R_{4}\)	\(F_{3}=[mRNA_{FXR}×Ts_{FXR}]\) \(F_{4}=[mRNA_{RXR}×Ts_{RXR}]\)	In E. coli, the rate of transcription is basically 10-100 nt/s, the elongation rate is nearly 5-8 fold higher than that of E. coli ternary complexes, so the number of mRNA molecules transcribed into DNA per second could be calculated according to the length of mRNA bases of FXR-ddRFPA1 and RXR-ddRFPB1.
\(R_{5} \& R_{6}\)	\(F_{5}=[mRNA_{FXR}]×Dm_{FXR}\) \(F_{6}=[mRNA_{RXR}]×Dm_{RXR}\)	The degradation rates of mRNA and protein are related to the half-life and linearly related to their quantity. Therefore, the corresponding degradation rate can be obtained by the half-life time (taking molecules/s as the unit).
\(R_{7} \& R_{8}\)	\(F_{7}=[FXR]×Ds_{FXR}\) \(F_{8}=[RXR]×Ds_{RXR}\)	The degradation rates of mRNA and protein are related to the half-life and linearly related to their quantity. Therefore, the corresponding degradation rate can be obtained by the half-life time (taking molecules/s as the unit).
\(R_{9}\)	\(F_{9}=[ CDCA]×[FXR]×k_{1}\)	According to the concentration and volume provided by the Wet Lab, the initial number of CDCA corresponding to different concentrations of solution could be calculated according to the number of molecules = concentration*volume*NA.
\(R_{10}\)	\(F_{10}=[ CDCA-FXR]×[RXR]×k_{2}\)
\(R_{11}\)	\(F_{11}=[FXR]×[RXR]×k_{3}\)

Solving steps

1)Initialize system conditions,t=0.

2)Calculate tendency function \(a_{j}(j=1,2,…,M)\), calculate \(a_{0}(x)=\sum a_{j}(x)\).

3)We take \(0 <\epsilon < 1\).

4)\(\tau = \underset{j \in[1, M]}{\min}{\begin{Bmatrix}\epsilon a_{0}(x)\over|\sum_{i=1}^N \xi_{i}(x)b_{ji}(x)|\end{Bmatrix}} \)

5)Random numbers \(K_{j}\) are generated according to poisson distribution \(P_{j}(a_{j}(x),\tau)\) in each step, and the increment of system state is defined as \(\lambda=\sum_{j=1}^M K_{j}v_{j}\).

6)Update time \(t=t+\tau\), system state \(X(t+\tau)=X+\lambda\).

7)Iteration: When t is less than the specified time range, turn to the second step. Otherwise stop.

8)Plot and analyze the results.

Through the stochastic model we constructed above based on \(\tau-leaping\), and solved the following two problems respectively:

1. The content of each substance involved in the simulation reaction

2. Selection of optimal initial concentration of CDCA

①The content of each substance involved in the simulation reaction

Taking 50μmol/L as an example, the content of different substances with time simulated by \(\tau-leaping\) is shown as follow:

②Selection of optimal initial concentration of CDCA

The different concentrations of CDCA will affect the content of the product CDCA-FXR-RXR. In order to maximize the product content, we first simulated the content of CDCA-FXR-RXR under the initial CDCA concentration of 0μmol/L, 10μmol/L, 25μmol/L, 30μmol/L, 50μmol/L, 70μmol/L and 100μmol/L. The optimal CDCA concentration in 0-100 μmol/L was obtained by smoothing spline fitting of the above five groups of data.

The five groups of data obtained by fitting:

Table 3. Molecular number of CDCA-FXR-RXR at different concentrations of CDCA

Initial concentration of CDCA	0	10	25	30	50	70	100
Molecules	0	1374.018182	1469.763636	1453.672727	1445.945455	1442.545455	1446.945455

Smoothing spline was used to fit the data, and the curve after fitting：

Through calculation, the best initial concentration of CDCA is 25.1μmol/L. And the number of molecules is 1595.36509215649.

Stability and robustness analysis

Stability analysis

Firstly, according to the stochastic model based on \(\tau-leaping\) algorithm, we made stability analysis to \(Tc_{FXR}\) and \(Tc_{RXR}\). Taking 25μM for example, we gave them definite floating value and solved the model, then normalized the results of all the substances and compared them with change of parameters. Eventually we got:

Table 4. Ratio of change of stability analysis

Change of mRNA/Change of parameter	-0.002	-0.001	0	0.001	0.002
-0.002	-6.280870045	-2.549528644	0	-6.904973411	-10.60975722
-0.001	-0.18590313	-6.63939751	0	15.77520848	-16.22668752
0	0	0	1	0	0
0.001	-13.14600707	-6.32070643	0	-38.6678511	8.524986403
0.002	1.845752508	1.646570583	0	8.578101583	-9.149089769

From the analysis above, our model has strong stability.

Robustness analysis

We also made a robustness analysis to prove the strong robustness of our model. When solving the model, we added a random perturbation within the scope of 0.1 for \(Tc_{FXR}\) and \(Tc_{RXR}\). Similar to the stability analysis above, we got:

Table 5. Ratio of change of robustness analysis

Parameter floating value	-0.01	-0.005	0	0.005	0.01
-0.01	-0.057909622	-0.023506654	-0.00012243	-0.063663855	-0.097821962
-0.005	-0.000857013	-0.030607623	-0.124022086	0.072723711	-0.074805029
0	0.006121525	0.145202561	0	0.095128491	0.036606717
0.005	-0.060603093	-0.029138457	-0.016160825	-0.178258794	0.039300187
0.01	0.017017838	0.015181381	0.064520868	0.079090097	-0.084354608

From the analysis above, our model has strong robustness.

Ordinary Differential Equation

In conclusion, we used the \(\tau-leaping\) algorithm to simulate the operation process of the system. While in real production, small disturbances caused by environmental factors could exist anytime and anywhere. So we used Ordinary Differential Equations to simulate the Cell-free system and gave small perturbations to each parameter to explore the extent of impact to the reaction system affected by the external environment.

We used Ordinary Differential Equations to simulate the Cell-free system and gave small perturbations to each parameter to explore the extent of impact to the reaction system affected by external environment.

For the process of transcription and translation, we adopt the Michaelis-Menten equation and Law of mass action. Michaelis-Menten equation is a velocity equation that expresses the relationship between the initial velocity of an enzymatic reaction and the substrate concentration. The law of mass action tells us that the rate of a chemical reaction is proportional to the effective mass of the reactant. Therefore, the production and degradation equations of \(mRNA_{FXR}\) were simulated as follow:

\( \frac{d[mRNA_{FXR} ]}{d t}=\frac{V_{m}[DNA_{FXR}]}{K_{m}+[D N A_{FXR}]}-Ts[mRNA_{FXR}] \)

Sortase was used in the translation process of \(mRNA_{FXR}\). Under the assumption that the half-life of the translated protein per unit time is constant regardless of the concentration of the enzyme or its own concentration, using the first-order kinetic model, the translation equation of \( mRNA_{FXR} \) was as follow:

\( \frac{d[F X R]}{d t}=T s[ m R N A_{FXR}]-k_{1} E_{0}[F X R] \)

RXR is in the same way:

\( \begin{aligned} \frac{d[m R N A_{RXR}]}{d t} &=\frac{V_{m}[D N A_{RXR}]}{K_{m}+[D N A_{RXR}]}-T s[m R N A_{R X R }] \\ \frac{d[R X R]}{d t} &=T s[m R N A_{R X R }]-k_{2} E_{0}[R X R] \end{aligned} \)

During the experiment, FXR was used to combine with CDCA to form a complex. This complex bound to RXR to form a dimer. The process was as follow:

\( \begin{array}{c} CDCA + FXR \stackrel{k_{1}}{\longrightarrow} CDCA-FXR \\ CDCA-FXR + RXR \stackrel{k_{2}}{\longrightarrow} CDCA-FXR-RXR \end{array} \)

The rate of change for CDCA, FXR, and RXR are as follow:

\( \begin{array}{c} V_{C D C A}=-\frac{d_{C_{0}}}{d_{t}}=\frac{V_{\max C D C A} C_{0}}{K_{m}}=\frac{K_{c a t} E_{0} C_{0}}{K_{m}+C_{0}} \\ V_{F X R}=-\frac{d_{C_{0}}}{d_{t}}=\frac{V_{\max F X R} C_{0}}{K_{m}}=\frac{K_{c a t} E_{0} C_{0}}{K_{m}+C_{0}} \\ V_{R X R}=-\frac{d_{C_{0}}}{d_{t}}=\frac{V_{\max R X R} C_{0}}{K_{m}}=\frac{K_{c a t} E_{0} C_{0}}{K_{m}+C_{0}} \end{array} \)

\( V_{C D C A}=\frac{d_{C_{1}}}{d_{t}}=\frac{V_{\max C D C A} C_{0}}{K_{m 1}}-\frac{V_{\max R X R} C_{0}}{K_{m}}=\frac{K_{c a t} E_{0} C_{0}}{K_{m 1}+C_{0}}-\frac{K_{c a t} E_{0} C_{0}}{K_{m}+C_{0}} \)

The rate of change in FXR-CDCA-RXR is as follow:

\( V_{\text {FXR-CDCA-RXR }}=\frac{d_{C_{0}}}{d_{t}}=\frac{V_{\max \text { CDCA }-\mathrm{FXR}-\mathrm{RXR}} C_{0}}{K_{m}}=\frac{K_{\mathrm{CDCA}-\mathrm{FXR}-\mathrm{RXR}} E_{0} C_{0}}{K_{m}+C_{0}} \)

Stability Analysis

In order to prove that the reaction could proceed stably under the general environment, we added a disturbance to each reaction channel (The concentration of CDCA: 25μmol/L). The disturbance represented the changes of temperature, humidity and other environmental factors in the interior. Seven values of the deviation from the standard environment -15%, -10%, -5%, 0, 5%, 10% and 15% were taken as the disturbance frequency. Through simulation, the change of molecular number after adding disturbance was:

We took \(2.88×10^{4} s\) as the time node, the disturbance deviation as the abscissa, and the number of CDCA-FXR-RXR molecules at the time node as the ordinate. Make the following diagram:

In conclusion, we found that under different perturbations, the number of CDCA-FXR-RXR molecules changed in the same trend, and the number of CDCA-FXR-RXR molecules tended to be the same at the same time node. It indicated that the reaction could be carried out stably under different indoor environments. Therefore, our reaction had strong stability.