Team:NAU-CHINA/Model/SeNPs Production Model

SeNPs Production Model

1 Abstract

2 Symbolic description & Hypothesis

3 Growth of E.coli

4 SeNPs yield of a single E. coli

5 Selection of promoter combinations

6 Turning waste into treasure

7 Sensitivity and robustness analysis

1 Abstract





Not every process in gene circuit can be quantitatively observed and monitored, so the response of biological system cannot be accurately predicted. Therefore, it is important to verify the feasibility of experiments and provide guidance through mathematical modeling. In this section, we build the E.coli growth model, deterministic model and stochastic simulation model with MATLAB.

As for E.coli growth, we use Logistic model and introduce decay rate. Since SeNPs is mainly produced during the logarithmic growth phase of bacteria, we determine when bacteria finish this phase and the bacterial concentration at the corresponding time according to the grow data. The productivity of single bacterium is determined by stochastic simulation, using tau-leaping algorithm. In our design, gene ssuE and sefA don’t share the same promoter. At different intensities of the promoter combinations, the SeNPs yield under the corresponding conditions is used as an evaluation criterion to determine the optimal promoter combination. Finally, combined with the deterministic equation, the effect of "turning waste into treasure" is explored with soil extracts as raw material.

Figure 1.1 Flow chart

2 Symbolic description & Hypothesis





Symbolic description

Symbol Explanation Value Units
PART 1 Growth of E.coli
$N_{0}$ Initial OD600 value
$K$ Saturated OD600 value
$r$ Growth rate
$d$ Decay rate
PART 2 SeNPs yield of a single E. coli [In the cytoplasm]
$N_{plasmid}$ Average copy number of plasmids per cell[1] 200 Molecule
$Tc_{ssuE}$ Transcription rate of the DNA of ssuE[2] s-1
$Tc_{sefA}$ Transcription rate of the DNA of sefA[2] s-1
$Ts_{ssuE}$ Translation rate of the mRNA of ssuE[2] 0.02079 s-1
$Ts_{sefA}$ Translation rate of the mRNA of sefA[2] 0.10438 s-1
$Dm_{ssuE}$ Degradation rate of the mRNA of ssuE[3] Half life:3.96min s-1
$Dm_{sefA}$ Degradation rate of the mRNA of sefA[3] s-1
$D_{SsuE}$ Degradation rate of the SsuE[4] Half life: 15min s-1
$D_{SefA}$ Degradation rate of the SefA[4] s-1
$K_{SsuE}$ The rate at which the enzyme SsuE will reduce $SeO_{3}^{2-}$ to Se0 0.0003 s-1
$K_{wrapping}$ The rate at which the protein SefA wraps Se0 into SeNPs 0.00001 s-1
PART 3 Turning waste into wealth [In the periplasm]
$K_{SerABC}$ The rate at which enzyme SerABC converts selenate to selenite 0.027
$V_{enter}$ The rate at which ions from soil extracts enter the periplasm[5] 0.0238
$\alpha$ The rate at which selenate and selenite decrease in the periplasm 0.0000038

Hypothesis

·The copy number of plasmids is kept as a constant.
·The degradation rates of mRNA and protein are constant.
·The molecular species are uniformly distributed inside the cell.
·Some of the chemical reactions are at equilibrium or steady state.
·We don’t consider the effects of other reactions on Se0 or SeNPs.
·Ignore the phenomenon that selenate ions directly enter the cytoplasm from the periplasm.

3 Growth of E.coli





Logistic model

The Logistic model was first proposed by the Dutch biological mathematician P.F.Verhaulst in 1837. It was used to describe the law of population growth of organisms. In order to make the model more consistent with the growth trend of E.coli, decay rate d is introduced to describe the decline trend of the population of E.coli cells after decay stage, which makes the growth curve more accurate.

The comparisons between the traditional Logistic model and the new Logistic model are as follows:

Traditional:

Introducing the decay rate d:

where $N_{0}$ represents the initial OD600 value of E.coli, $K$ represents the saturation OD600 value, $r$ represents the growth rate, and $d$ represents the decay rate.

The growth data of E.coli in 12, 15 and 18mM sodium selenite solution are recorded. Among three groups of sodium selenite solutions with different concentrations, each group is further divided into three groups according to the low, middle and high intensity of promoter, with a total of nine groups of data. Logistic parameters of each group of data are obtained through fitting as follows:

Table 1 E. coli growth fitting data
Concentration(mM) Intensity K(Saturation OD600 value) N0(Initial OD600 value) r(Growth rate) d(Decay rate) R-square
12 Low 0.6973 0.045 0.5509 3.69E-11 0.9896
Middle 0.73 0.057 0.5621 4.25E-09 0.9926
High 0.7397 0.059 0.5876 3.24E-09 0.9818
15 Low 0.6911 0.045 0.4216 1.04E-10 0.9908
Middle 0.775 0.057 0.5038 3.78E-09 0.9747
High 0.7735 0.059 0.5281 2.53E-08 0.9719
18 Low 0.8327 0.045 0.3381 1.92E-09 0.9683
Middle 0.635 0.057 0.4177 5.50E-11 0.985
High 1.022 0.059 0.3256 6.66E-09 0.9541

Figure 3.1 Growth curve of E.coli

Through the nine groups above, it is observed that E.coli finish the logarithmic growth phase,during which it produces SeNPs, at about 9.5h. Therefore, the corresponding OD600 value at this time is selected as the final concentration.

According to the literature, the relationship between OD600 value and the population of E.coli is as follows:

$$y=8708.7e^{36.43x}$$

where y represents the population of bacteria and x represents OD600 value. Then according to every fitting equation, the results are as follows:

Table 2 Population of E.coli in different conditions
Concentration(mM) Intensity OD600 value of 9.5h Number(*1015)
12 Low 0.6472 0.15147
Middle 0.6909 0.74227
High 0.7089 1.4319
15 Low 0.5478 0.0040437
Middle 0.7013 1.0844
High 0.7161 1.8573
18 Low 0.4884 0.00046419
Middle 0.5328 0.0023449
High 0.5873 0.017024

4 SeNPs yield of a single E. coli





Why do we choose tau-leaping algorithm?

The reactions of producing SeNPs take place in the cytoplasm and the number of reaction molecules involved is usually relatively small, causing obvious randomness. So we need something different than the reaction rate equation in 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+dtau], 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 choose "tau-leaping algorithm" improved from Gillespie algorithm.

Different from Gillespie algorithm where the main question is "Which reaction channel will fire in the next specified time interval τ?"; the tau-leaping algorithm asks the other question, "How often does each reaction channel fire in the next specified time interval τ?"

Background knowledge of tau-leaping

When a system has M reaction channels {$R_{1}$,……,$R_{M}$}and N substances {$S_{1}$,……,$S_{N}$}, the state vector $X(t)=(X_{1}(t),……,X_{N}(t))^{T}$ can represent its dynamic state at time t, and $X_{i}(t)(i=1,...,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}=(V_{1j},…, V_{Nj})^{T}$ represents the state change vector, where $ V_{ij}$ represents molecular changes number of substance $S_{i}$ when $R_{j}$ occurs.
For reaction $R_{1}$ :

$$DNA_{ssuE}\xrightarrow[\quad\quad]{Tc_{ssuE}}{mRNA_{ssuE}}$$

Its state change vector is $V_{1}=[0\,0\,1\,0\,0\,0\,0\,0\,0]$ ,
and tendency function $a_{1}=N_{plasmid}×Tc_{ssuE}$ .

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

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

where $\xi(x)=\sum_{j=1}^M a_{j}(x)V_{j}$,$a_{0}=\sum_{j=1}^M a_{j}(x)$, $b_{ji}(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}\begin{Bmatrix}a_{j}(x),\tau\end{Bmatrix}$ 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_{ji}P_{j}(a_{j}(x),\tau)(i=1,…,N)$$

Figure 4.1 Gene circuit design

Table 3 Reactions
Reaction channel
$R_{1}$ $DNA_{ssuE}\xrightarrow[\quad\quad] {Tc_{ssuE}}mRNA_{ssuE}$
$R_{2}$ $DNA_{sefA}\xrightarrow[\quad\quad] {Tc_{sefA}}mRNA_{sefA}$
$R_{3}$ $mRNA_{ssuE}\xrightarrow[\quad\quad] {Ts_{ssuE}}ssuE$
$R_{4}$ $mRNA_{sefA}\xrightarrow[\quad\quad] {Ts_{sefA}}sefA$
$R_{5}$ $mRNA_{ssuE}\xrightarrow[\quad\quad] {Dm_{ssuE}}\emptyset$
$R_{6}$ $mRNA_{sefA}\xrightarrow[\quad\quad] {Dm_{sefA}}\emptyset$
$R_{7}$ $SsuE\xrightarrow[\quad\quad] {D_{SsuE}}\emptyset$
$R_{8}$ $SefA\xrightarrow[\quad\quad] {D_{SefA}}\emptyset$
$R_{9}$ $SeO_{3}^{2-}\xrightarrow[\quad\quad] {K_{SsuE} } Se^{0} $
$R_{10}$ $ Se^{0}+SefA\xrightarrow[\quad\quad]{K_{wrapping}}SeNPs$

According to the central rule,
$R_{1}$,$R_{2}$,$R_{3}$,$R_{4}$ represents the transcription and translation processes involved in enzyme SsuE and protein SefA synthesis;
$R_{5}$,$R_{6}$,$R_{7}$,$R_{8}$ represents the degradation process of four species($mRNA_{ssuE}$, $mRNA_{sefA}$, $SsuE$, $SefA$);
$R_{9}$ represents that $ SeO_{3}^{2-}$ converted to $Se^{0}$ by the enzyme SsuE;
$R_{10}$ represents the formation of SeNPs after $Se^{0}$ wrapped by SefA;
We don’t consider the degradation of $Se^{0}$ and $SeNPs$ because they are stable.

Determinations

State vectors: According to the reaction equations above, the substances involved in our model can be represented as state vectors:

X(t)=(DNAssuE, DNAsefA, mRNAssuE, mRNAsefA, SsuE, SefA, $SeO_{3}^{2-}$, Se 0, SeNPs)

State change vectors:

V1=[0 0 1 0 0 0 0 0 0]
V2=[0 0 0 1 0 0 0 0 0]
V3=[0 0 -1 0 1 0 0 0 0]
V4=[0 0 0 -1 0 1 0 0 0]
V5=[0 0 -1 0 0 0 0 0 0]
V6=[0 0 0 -1 0 0 0 0 0]
V7=[0 0 0 0 -1 0 0 0 0]
V8=[0 0 0 0 0 -1 0 0 0]
V9=[0 0 0 0 0 0 -1 1 0]
V10=[0 0 0 0 0 -1 0 -1 1]

Table 4 Tendency functions
The reaction channel Tendency functions Explanation
$R_{1}$&$R_{2}$ $F_{1}=N_{plasmid}×Tc_{ssuE}$ $F_{2}=N_{plasmid}×Tc_{sefA}$ $N_{plasmid}$:The plasmid copy number is roughly distributed in the interval [20,200], so we assume that the average plasmid copy number is a constant 200. $Tc_{ssuE}$+$Tc_{sefA}$:In E.coli, the rate of transcription is basically 10-100 nt/s[9], so the number of mRNA molecules that can be transcribed per second of DNA can be calculated according to the length of mRNA bases of SsuE and sefA.
$R_{3}$&$R_{4}$ $F_{3}=[mRNA_{ssuE}×Ts_{ssuE}]$ $F_{4}=[mRNA_{sefA}×Ts_{sefA}]$ $Ts_{ssuE}$+$Ts_{sefA}$:The intensity of RBS determines the number of amino acids that can be translated by the mRNA per second, so the number of protein molecules that can be translated by the mRNA per second can be calculated according to the number of amino acids of the proteins SsuE and SefA.
$R_{5}$&$R_{6}$ $F_{5}=[mRNA_{ssuE}]×Dm_{SsuE}$ $F_{6}=[mRNA_{sefA}]×Dm_{sefA}$ $Dm_{ssuE}$+$Dm_{sefA}$+$D_{SsuE}$+$D_{SefA}$: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}=[ SsuE]×D_{SsuE}$
$F_{8}=[SefA]×D_{SefA}$
$R_{9}$ $F_{9}$=$[ SeO_{3}^{2-}]$×$[SsuE]$×$K_{SsuE}$ $[SeO_{3}^{2-}]$:According to the concentration and volume provided by the wet lab, the initial number of $SeO_{3}^{2-}$ corresponding to different concentrations solution can be calculated according to the number of molecules = concentration*volume*NA.
$R_{10}$ $F_{10}=[ Se^{0}]×[SefA]×K_{wrapping}$

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=\mathop{min}_{j\in[1,M]}{\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 stochastic model constructed above which is based on tau-leaping, we will solve the following two problems respectively:

1.Selection of promoter combinations

2.SeNPs production under conditions of selenium-containing soil extracts

5 Selection of promoter combinations





Different concentrations of sodium selenite solution and different combinations of promoter both affect the population of E.coli and the ability of single cell to produce SeNPs. Promoter intensity will affect DNA transcription rate and RBS intensity will affect mRNA translation rate. In our degisn, ssuE and sefA don’t share the same promoter but use the same RBS. The corresponding rates of high(J23119), middle(J23101) and low(J23105) promoters are as follows:

Table 5 Transcription rates corresponding to different promoter combinations
ssuE sefA
J23105 J23101 J23119
J23105 0.00194\0.00972 0.00194\0.02835 0.00194\0.13889
J23101 0.00566\0.00972 0.00566\0.02835 0.00566\0.13889
J23119 0.02772\0.00972 0.02772\0.02835 0.02772\0.13889

Note: the transcription rates corresponding to the three promoters are 5.6nt/s(J23105), 16.33nt/s(J23101), 80nt/s(J23119), and the translation rate corresponding to the RBS is 20aa/s;

The mRNA bases numbers of ssuE and sefA are 2886nt and 576nt respectively, and the time required to produce one SsuE and SefA is 48.1s and 9.58s respectively.

Table 6 The number of molecules of selenite
Concentration Molecules
12mM(50mL) 3.612×1020
15mM(50mL) 4.515×1020
18mM(50mL) 5.418×1020

In this part, there are 3 concentrations of sodium selenite solution, 3 kinds of promoters, which are combined into a total of 27 groups.

Because the combinations are too many to discuss, we only study 9 groups with sodium selenite solution of 12mM in detail (the results of 15mM and 18mM are shown in the table 8 and table 9).

Next, according to the tau-leaping algorithm, we solve the stochastic model under 9 combinations and select the optimal combination of promoters according to the final SeNPs output. Under different combinations of 3 promoters (9 groups in total), we obtain the changes of mRNA of ssuE and sefA, SsuE, SefA and SeNPs over time through MATLAB.

Now we take one for example. In 12mM condition and low intensity promoters(J23105) of both ssuE and sefA, the results are as follows:

Figure 5.1 The expression of mRNA of 12mM

As can be seen from the Figure 5.1, the number of mRNA molecules of both ssuE and sefA increase and then gradually stabilize under the combined function of transcription, translation and degradation, basically stable at 15~25.

Figure 5.2 The expression of protein of 12mM

As can be seen from the Figure 5.2, the molecules of SsuE are around 220, and that of SefA are roughly stable at 5-10.

SsuE has been increasing, because as an enzyme, SsuE does not consume when participating in the reaction, and the mRNA translation rate of ssuE is much higher than its own degradation rate. However, SefA shows a trend of firstly increasing but then decreasing, because Se0 is not generated at the initial stage of the reaction and is not wrapped by SefA. With the increasing of Se0, SefA begins to decrease because it is involved in the formation of SeNPs as wrapping protein of Se0.

Figure 5.3 The expression of elemental selenium and SeNPs of 12mM

As can be seen from Figure 5.3, as the reactions goes on, Se0 is almost not produced during the first 2000s. Then it increases to around 300000 and then gradually decreases. SeNPs is not produced during the first 2000s because there is no Se0,then it is produced after Se0 begins to be generated.

As you can see from all the results above:

A.We have a lot of jumping points which indicates that our algorithm has great efficiency.

B.The trend of all substances corresponds with reality.

Then we remove all the jumping points except the effective points and put nine combinations of promoters together for analysis and comparison. The results are as follows:

Figure 5.4 The expression of mRNA of 9 promoter combinations

As is shown in the Figure 5.4, some lines seem long while some seem short. The reason is that in some conditions we have more jumping points while in other conditions we have less jumping points. Long lines indicate more jumping points and smaller time interval. We can see that 9 groups have the same trend of mRNA but finally reach three levels.

Figure 5.5 The expression of protein of 9 promoter combinations

As is shown in the Figure 5.5, the molecules of SsuE reach three levels because of different promoter intensities, and the molecules of SefA firstly increase but then decrease. Moreover, the promoter intensity of ssuE will affect the molecules of SefA.

Figure 5.6 The expression of SeNPs of 9 promoter combinations

As is shown in the Figure 5.6, SeNPs can reach more than 20000 molecules in the conditions of 12mM in a single cell.

After communicating with the wet lab, we know that there is an upper production limit of SsuE in a single E.coli because excessive SsuE will increase the burden of E.coli and make it hard to survive. The upper limit is around 500, so we only discuss the combinations with a theoretical yield lower than this limit. The feasible combinations are: lowssuE +lowsefA, lowssuE+middlesefA, lowssuE+highsefA. (The same is true for 15mM and 18mM) The results are presented to the wet lab:engineering

Then combined with the population of bacteria, we obtain the total yield of SeNPs (total yield = population of E.coli *product amount produced by single E.coli) as follows:

Table 7 Toal yield of SeNPs in 12 mM(1015)
ssuE sefA
J23105 J23101 J23119
J23105 175.06 1508.07 13231.34


Table 8 Toal yield of SeNPs in 15 mM(1015)
ssuE sefA
J23105 J23101 J23119
J23105 4.84 1913.13 15790.73


Table 9 Total yield of SeNPs in 18 mM(1015)
ssuE sefA
J23105 J23101 J23119
J23105 0.62 4.90 152.41

It can be seen from Table 7, Table 8 and Table 9.

Total yield in 18 mM is much less than 12mM or 15mM. Probably because selenite overloads bacteria. In conclusion, for the 9 combinations,there are different SeNPs yields of single E.coli and different concentrations of E.coli.

The feasible combination that can achieve the maximum SeNPs yield is J23105\J23119 [low(ssuE)\high(sefA)]. However, only J23105\J23105 [low(ssuE)\low(sefA)] combination is the best choice because it can get the best quality of SeNPs(The smaller the particle size, the better the quality). The details for this can be seen in the results.

In our model, in [12mM low(ssuE)\low(sefA)] conditions, the final theoretical yield of SeNPs is 0.035g, which is calculated by the relative molecular mass of SeNPs. And the actual production from wet lab is 0.0268g. We have achieved mutual verification between wet lab and model.

Therefore, in the following part “Turning waste into treasure” -- discussing soil extracts, we choose this combination to study the actual situation.

6 Turning waste into treasure





Transformation of selenate($SeO_{4}^{2-}$) in periplasm

As for the selenium contaminated areas, the soil extracts contains selenate ($SeO_{4}^{2-}$) and selenite($SeO_{3}^{2-}$) ions. If the soil extracts is used as raw material, selenate and selenite firstly enter the periplasm from the soil extracts, and then selenate convert to selenite under the function of enzyme SerABC, and then selenite enters the cytoplasm and participates in the reactions to produce SeNPs.

Figure 6.1 The reactions in periplasm

As for reactions in the periplasm, considering that they have a certain state continuity and certainty, the Ordinary Differential Equation Model is appropriate to study the change of the concentration of various substances over time, and then describe the biochemical behavior of this part.

$${d[SeO_{4\,soil}^{2-}]\over dt} = -V_{enter}[SeO_{4\,soil}^{2-}]$$ $${d[Se_{3\,soil}^{2-}]\over dt}=-V_{enter}[SeO_{3\,soil}^{2-}]$$ $${d[SeO_{4\,periplasm}^{2-}]\over dt}=V_{enter}[SeO_{4\,soil}^{2-}]-K_{SerABC}[SeO_{4\,periplasm}^{2-}]-\alpha[SeO_{4\,periplasm}^{2-}]$$ $${d[SeO_{3\,periplasm}^{2-}]\over dt}=V_{enter}[SeO_{3\,soil}^{2-}]+K_{SerABC}[SeO_{4\,periplasm}^{2-}]-\alpha[SeO_{3\,periplasm}^{2-}]$$

where Venter is the transport rate of two ions from the soil extracts into the periplasm, and KSerABC is the conversion rate of selenate to selenite under the function of enzyme KSerABC, which is related to the concentration and activity of the enzyme, but we can sum it up as a parameter here. This part refers to the Law of Mass Action where the reaction rate is proportional to the concentration of reactants.

In order to make the simulation more suitable for periplasm, we introduce a general parameter $\alpha$ to describe the reduction of selenate and selenite caused by various reasons, such as other biochemical reactions and changes in properties and morphology in periplasm.

Since enzyme SerABC exists in E.coli itself, it is assumed that its concentration is sufficient before and after participating in related reactions. Therefore, we ignore its transcription, translation and degradation processes, and focus on the study of selenate and selenite concentration changes.

In the Ordinary Differential Equation Model of periplasm part, selenate and selenite ion concentrations of soil extracts are selected as 0.073ug/mL and 1.08ug/mL after reading literautre[10], and the changes of various substances over time in this part are obtained by solving the model as follows:

Figure 6.2 Ion concentration variation in periplasm

As shown in the Figure 6.2, the transformation of selenate is basically completed at 660s, and the total concentration of selenite is 14.568mmol/L. In order to make the results more intuitive, 50mL soil extracts is selected as the object of further study, and the molecular of selenite entering cytoplasm is $4.385×10^{20}$.

Formation of selenium-nanoparticles in cytoplasm

The result $4.385×10^{20}$ above is taken as the initial molecular number of selenite in the cytoplasm. It is estimated that orders of magnitude of E.coli is about $10^{14}$~$10^{16}$. Therefore, a single cell can use $4.385×10^{20}$ selenite. Choosing[50 ml J23105\J23105][low(ssuE)\low(sefA)] as object of the study, the final SeNPs yield is obtained as shown in the Figure 6.3.

Figure 6.3 The yield of SeNPs in single E.coli

It can be seen from the Figure 6.3 that in a certain time range and under more complicated natural conditions, SeNPs successfully complete production, and the output increase steadily and finally reach 1179.09, indicating that our experiment and established model are suitable for more practical and complex situations and can be further extended.

7 Sensitivity and robustness analysis





Stability analysis

Firstly, according to the stochastic model based on tau-leaping algorithm, we make stability analysis to TcssuE and TcsefA. Take(12mM) for example, we give them definite floating value and solve the model, then normalize the results of all the substances and compare with change of parameters. We get:

Figure 7.1 Stability analysis

Table 10 Ratio of change of stability analysis
Change of mRNA/Change of parameter -0.002 -0.001 0 0.001 0.002
-0.002 2.8441 7.7614 0 6.0175 5.3627
-0.001 2.7760 18.4493 0 3.7946 31.1056
0 0 0 1 0 0
0.001 2.0732 31.2029 0 5.4562 14.3253
0.02 6.0325 16.1927 0 1.0040 1.2125

From the analysis above, our model has strong stability.

Robustness analysis

We also make robustness analysis to prove the strong robustness of our model. When solving the model, we add a random perturbation within the scope of 0.05 for TcssuE and TcsefA. The same as above, we get:

Figure 7.2 Robustness analysis

Table 11 Ratio of change of robustness analysis
Parameter floating value -0.01 -0.005 0 0.005 0.01
-0.01 0.0131 -0.0018 -0.0015 -0.0005 -0.0025
-0.005 0.0505 -0.0011 -0.0024 -0.0024 -0.0001
0 -0.01941 -0.0013 0 0.0010 -0.0057
0.005 0.0261 0.0018 -0.0059 -0.0041 -0.0030
0.01 0.0253 -0.0031 -0.0018 -0.0058 -0.0072

From the analysis above, our model has strong robustness.




[1]Zhou T. Methods for determination of plasmid copy number[J]. Biotechnology Newsletter, 1999(01): 51-57.

[2]http://parts.igem.org/Part:BBa_J23105.

[3]Bernstein J A, et al. Global analysis of Escherichia coli RNA degradosome function using DNA microarrays[J]. Proceedings of the National Academy of Sciences, 2004, 101(9): 2758-2763.

[4]Li C. Synthetic biology[M]. Beijing: Chemical Industry Press, 2019: 182.

[5]Terence A. Brown,Alex Shrift. Selective assimilation of selenite by Escherichia coli[J]. NRC Research Press Ottawa, Canada,1982,28(3):307-309.

[6]Cao YY. Model fitting and cluster analysis of bacterial growth curve[D]. Shanghai: East China Normal University, 2020.

[7]Lin J, He J Q, Shang C L, et al. Study on the correlation between the OD value of 5 common pathogenic bacteria and the number of viable bacteria[J]. Journal of Tianjin Agricultural College, 2017, 24(04): 42-45.

[8]Gillespie D T. Approximate accelerated stochastic simulation of chemically reacting systems[J]. The Journal of chemical physics, 2001, 115(4): 1716-1733.

[9]R.Milo and R.Phillips. Cell Biology by the Numbers[M]. First edition,2015. ISBN9780815345374.

[10]Jiang G B, et al. Extraction and ion chromatographic analysis of $SeO_{3}^{2-}$ and $SeO_{4}^{2-}$ in soil[J]. Analytical laboratory, 1994(04): 59-61.