Introduction
In our experiments, we build biosensors responds to different levels of ambient lactate, oxygen and pH and is regulated by the transcription factors LldR, FNR and CadC, respectively that use RFP as reporter protein. Then, we added the lysis gene to these biosensors to give out proteins inside the flora. With the purpose of trying to establish models that describe the RFP- time curve of biosensor and biosensor flora number – time curve, we firstly build an ordinary differential equation model based on the dynamics of gene regulation and the input function of a gene regulated was established to describe and predict the RFP- time curve. Then, based on these above reaction mechanisms, we then Consider the impact of lysis genes to modify the traditional logistics equation in order to better reflect the RFP- time curve of biosensor flora number.
Model Assumptions
To establish a more reasonable mathematical model, the following assumptions are made:
1. The timescales for binding and transcription reactions are substantially faster than translation.
2. Full contact of all substances involved in the reaction process.
3. The enzymatic reaction process conforms to the Michaelis-Menten equation.
4. The binding reaction of the transcription factor to the environmental signal and the binding reaction of the transcription factor to the promoter will reach equilibrium in a short period of time and can be regarded as a balancing process in subsequent analysis.
5. The growth of bacteria conforms to the assumptions of the logistics equations
Model Building
A model which represents multi-component and time-varying dynamic system is widely used in various biological problems. Among all of the models, differential equation is the top of the list. To illustrate it, the regulatory network can be represented by an array of ordinary differential equations, where the interaction between a multitude of molecules (such as mRNA or protein) are quantified by the law of mass action[1]. These equations designate the level of each protein or mRNA as a function of other components in the development of the system. These models usually include time-related variables, such as the concentration of protein and mRNA, other parameters, such as protein degradation parameters, are also included.
Reaction mechanism
We used oxygen, pH, and lactate as unique indicators of colorectal cancer tumor tissue to construct bacterial biosensing systems capable of distinguishing unique organ environments. To represent the behavior of the biosensor, we used ordinary differential equations to describe the production rate of the red fluorescent protein (RFP) controlled by different biosensor modulators, including activators for the hypoxia and pH biosensors and repressors for the lactate biosensors.
For hypoxia biosensors, we used a hypoxia-sensing promoter (pPepT) that is primarily regulated by the transcriptional activator, fumarate and nitrate reduction regulatory protein (FNR) In the absence of oxygen, FNR generates a transcriptionally active homodimer. However, the cluster is degraded in the presence of oxygen, which dissociates the FNR dimer into inactive monomers. (fig 1)
Thus, by measuring red fluorescent protein (RFP) expressed under the control of the pPepT on a plasmid, we detected elevated levels of fluorescence in response to hypoxic conditions.
For lactate biosensor, We derive it from the native lldPRD operon, to detect lactic acid fermentation by host mammalian cells[3]. This lactate-sensing system was constructed on two plasmids: 1. a lactate-inducible reporter plasmid which drives expression of a gene of interest. 2. a repressor plasmid, which produces the repressor LldR that dimerizes to inhibit expression of the reporter gene unless bound to lactate. (fig2)
Hence, in response to increasing concentrations of lactate in the culture medium, we observed a corresponding increase in RFP.
For pH biosensor, we engineered the pH-sensitive promoter pCadC among other systems[4] that is regulated by a membrane-tethered activator protein (CadC), which shows increased activity in acidic medium compared with medium at a neutral pH.
Therefore, comparing with medium at a neutral pH, we observed a relative high level of fluorescence of RFP in the acidic medium.(fig3)
Based on the above reaction mechanism, we will first write the dynamics of the whole process and then build a mathematical modeling to both gain insight into and make predictions about their behavior.
The dynamics of the promoters
The dynamics of the promoters are given by the following set of reactions:
$$ \begin{equation} \begin{aligned} & FNR_2 \underset{k_{-O_2}}{\stackrel{k_{O_2}}{\rightleftharpoons}} 2FNR\\ & pPepT+FNR_2 \stackrel{k_{p1}}{\longrightarrow} pPepT^+\\ & LldR_2 \underset{k_{-L}}{\stackrel{k_{L}}{\rightleftharpoons}} 2LldR\\ & pLldR+LldR_2 \stackrel{k_{l1}}{\longrightarrow} pLldR^-\\ & pCadC+CadC \stackrel{k_{c1}}{\longrightarrow} pCadC^+\\ \end{aligned} \end{equation} $$Where: $ pPepT^+ $ represent the states of pPepT bound with activator.
$ pLldR^- $ represent the states of pLldR bound with repressor.
$ pCadC^+ $ represent the states of pCadC bound with activator.
$ FNR $ represents FNR monomer.
$ FNR_2 $ represents FNR dimer, which is the activator of $ pPepT $.
$ LldR $ represents LldR monomer.
$ LldR_2 $ represents LldR dimer, which is the repressor of $ pLldR $.
$ CadC $ represents CadC monomer.
$ CadC_2 $ represents CadC dimer, which is the activator of $ pCadC $.
That is the dynamics of the promoters.
The dynamics of the transcription, translation
Then we move to the set of reactions govern the transcription, translation of each gene/protein:
$$ \begin{equation} \begin{aligned} pPepT^+ & \stackrel{k_{p2}}{\longrightarrow} pPepT^+ + m_{O_2}\\ pPepT & \stackrel{k_{p2'}}{\longrightarrow} pPepT + m_{O_2}\\ m_{O_2} + & \stackrel{k_{p3}}{\longrightarrow} m_{O_2}+RFP\\ pLldR & \stackrel{k_{l2}}{\longrightarrow} pLldR + m_{L}\\ m_{L} & \stackrel{k_{l3}}{\longrightarrow} m_{L}+RFP\\ pCadC^+ & \stackrel{k_{c2}}{\longrightarrow} pCadC^+ + m_{pH}\\ pCadC & \stackrel{k_{c2'}}{\longrightarrow} pCadC + m_{pH}\\ m_{pH} & \stackrel{k_{c3}}{\longrightarrow} m_{pH}+RFP\\ \end{aligned} \end{equation} $$where $m_{O_2}$ represents the mRNA which is transcribed of gene segment with pPepT as its promoter;
$m_{L}$ represents the mRNA which is transcribed of gene segment with pPepT as its promoter;
$m_{pH}$ represents the mRNA which is transcribed of gene segment with pPepT as its promoter;
$ RFP $ represents the translation product: fumarate and nitrate reduction regulatory protein(RFP);
The dynamics of the decay
Finally, the mRNA, protein, and the promoter which is bound with activator/repressor decay exponentially, and all tagged forms of decay enzymatically according to the following set of reactions:
$$ \begin{equation} \begin{aligned} RFP &\stackrel{\gamma}{\longrightarrow} \varnothing\\ m_{O_2} &\stackrel{\gamma_{p4}}{\longrightarrow} \varnothing\\ m_{L} &\stackrel{\gamma_{l4}}{\longrightarrow} \varnothing\\ m_{pH} &\stackrel{\gamma_{c4}}{\longrightarrow} \varnothing\\ pPepT^+ &\stackrel{\gamma_{p5}}{\longrightarrow} pPepT\\ pLldR^- &\stackrel{\gamma_{l5}}{\longrightarrow} pLldR\\ pCadC^+ &\stackrel{\gamma_{c5}}{\longrightarrow} pCadC\\ \end{aligned} \end{equation} $$ODE model I: Expression of RFP gene
According to the above reaction mechanism, this reactions rates can be expressed using the Law of Mass Action. We use [x] to denote the concentration of substance x. Therefore, constitutive expression of RFP is mathematically described by an ordinary differential equations (ODE) that captures the species dynamics.
Below is the reaction dynamics of the hypoxia biosensor
$$ \begin{equation} \left\{ \begin{aligned} &\frac{d\left[FNR_2\right]}{d t}=k_{-O_2}[FNR]^2-k_{p_1}[pPepT]\left[FNR_2\right] \\ &\frac{d\left[FNR\right]}{d t}=k_{O_2}\left[FNR_2\right] \\ &\frac{d\left[pPepT\right]}{d t}=-k_{p_1}[pPepT]\left[FNR_2\right]+\gamma_{p_5}\left[pPepT^{+}\right] \\ &\frac{d\left[pPepT^{+}\right]}{d t}=k_{p_1}[pPepT]\left[FNR_2\right]-\gamma_{p_5}\left[pPepT^{+}\right] \\ &\frac{d\left[m_{O_2}\right]}{d t}=k_{p_2}\left[pPepT^{+}\right]+k_{p'_2}\left[pPepT^{+}\right]-\gamma_{p_4}\left[m_{O_2}\right] \\ &\frac{d[R F P]}{d t}=k_{p 3}\left[m_{O_2}\right]-\gamma[R F P] \end{aligned} \right. \end{equation} $$Next is the reaction dynamics of the lactic acid biosensor
$$ \begin{equation} \left\{ \begin{aligned} &\frac{d\left[LldR_2\right]}{d t}=k_{-L}[LldR]^2-k_{l_1}[pLldR]\left[LldR_2\right] \\ &\frac{d\left[LldR\right]}{d t}=k_{L}\left[LldR_2\right] \\ &\frac{d\left[pLldR\right]}{d t}=-k_{l_1}[pLldR]\left[LldR_2\right]+\gamma_{l_5}\left[pLldR^{-}\right] \\ &\frac{d\left[pLldR^{-}\right]}{d t}=k_{l_1}[pLldR]\left[LldR_2\right]-\gamma_{l_5}\left[pLldR^{-}\right] \\ &\frac{d\left[m_L\right]}{d t}=k_{l_2}\left[pLldR^{-}\right]-\gamma_{l_4}\left[m_L\right] \\ &\frac{d[RFP]}{dt}=k_{L3}\left[m_L\right]-\gamma[RFP] \end{aligned} \right. \end{equation} $$Last but not least, this is the reaction dynamics of the lactic pH biosensor
$$ \begin{equation} \left\{ \begin{aligned} &\frac{d\left[pCadC\right]}{d t}=-k_{c_1}[pCadC]\left[CadC\right]+\gamma_{c_5}\left[pCadC^{+}\right] \\ &\frac{d\left[pCadC^{+}\right]}{d t}=k_{c_1}[pCadC]\left[CadC\right]-\gamma_{c_5}\left[pCadC^{+}\right] \\ &\frac{d\left[m_{pH}\right]}{d t}=k_{c_2}\left[pCadC^{+}\right]+k_{c'_2}\left[pCadC^{+}\right]-\gamma_{c_4}\left[m_{pH}\right] \\ &\frac{d[R F P]}{d t}=k_{c 3}\left[m_{pH}\right]-\gamma[R F P] \end{aligned} \right. \end{equation} $$Simplification of Model
Because the change of RFP concentration over time is our most concerned index. However, the above equations of ordinary differential equations require a large amount of computation and a large number of solving parameters, and the time cost of complete simulation is too high. Therefore, we can simplify the above equations and only focused on the relationship between RFP concentration and time to reduce the computational cost [5].
Considering a gene regulated by a single regulator, this transcriptional interaction is described as “transcription factor X regulates gene Y”. In the absence of input signals, X is inactive and no Y is generated. Upon the presence of an external signal, the transcription factor rapidly changes to the active form and binds to the promoter of gene Y. Gene Y starts to be transcribed and mRNA is translated, leading to the accumulation of proteins. The concentration change of Y is described by the dynamic equation:
$$ \begin{equation} \frac{dy}{dt}=\beta-\gamma y. \end{equation} $$Here, $ \beta $ is the rate at which the gene produces protein Y, and $ \alpha $ is the sum of the degradation rate and dilution rate of the protein. y is the concentration of protein Y. The difference between production and degradation/dilution results in a change in the concentration of Y.
In some of the more accurate models, $ \beta $ would vary with the concentration of transcription factor change. This can be described by the strength of the effect of a transcription factor on the transcription rate of its target gene.
Let us consider first the production rate of protein Y controlled by a single transcription factor X. The number of molecules of protein Y produced per unit time $ f(x^*) $ is a function of the concentration of X in its active form $ X^*$. Then, the input functions of many real genes can be expressed in the form of Hill functions [6]:
$$ \begin{equation} f(x^*) = \begin{cases} & \frac{\beta}{1+(\frac{K}{x^*})^{n}},\quad (X^* \text{ is an activator)}\\ & \frac{\beta}{1+(\frac{x^*}{K})^{n}},\quad (X^* \text{ is an repressor)}\\ \end{cases} \end{equation} $$Where, $ \beta $ is the maximum expression level of the promoter;
$ n $ controls the steepness of the function;
$ x^* $ is the concentration of X in its active form $ X^*$;
$ K $ defines the concentration of active transcription factor, and it is related to factors such as the chemical affinity between a transcription factor and its site on the promoter.
Combining with above equations, we can get the dynamic equation that describes the concentration change of Y, namely:
$$ \begin{equation} \frac{dy}{dt} = \begin{cases} & \frac{\beta}{1+(\frac{K}{ x^*})^{n}}-\gamma y,\quad (X^* \text{ is an activator)}\\ & \frac{\beta}{1+(\frac{ x^*}{K})^{n}}-\gamma y,\quad (X^* \text{ is an repressor)}\\ \end{cases} \end{equation} $$Since transcription factor X requires an environmental signal to transition to the active state, in order to describe this process, we denote $x_T$ as the total transcription factor X concentration, $x$ as the transcription factor concentration X without binding environmental signal, and $x^*$ as the transcription factor concentration after binding environmental signal (namely active mode, $X^*$). S is the concentration of environmental signal in the environment.
From the law of conservation of mass, we have:
$$\begin{equation} x_T = x^* + x \end{equation}$$Let $k_{on}$ denote the binding rate of the environmental signal and the transcription factor, and $k_{on}$ denote the separation rate of the environmental signal and the transcription factor. We have:
$$\begin{equation} \frac{dx^*}{dt}= k_{on} xS - k_{off}x^*\end{equation}$$In general, $k_{on}>>k_{off}$, hence, the reaction will quickly achieve chemical equilibrium relation, at this point, $\frac{dx^*}{dt}=0$.[7] We can get:
$$\begin{equation} x^* = \frac{x_TS}{S+k_s}\end{equation}$$where, $k_s=\frac{k_{off}}{k_{on}}$ is the equilibrium constant of the signal-transcription factor binding reaction. From $ x^*=x_T-x $, we have:
$$\begin{equation} x = \frac{x_Tk_s}{S+k_s}\end{equation}$$Based on the above analysis, we developed a simplified system of ordinary differential equations that describe the dynamic expression of RFP(y) in response to varying levels of environmental lactate concentration (L; equation (1)), oxygen concentration (O2; equation (2)) and pH (H^+ concentration) (H+; equation (3)), regulated by transcriptional factors LldR, FNR and CadC, respectively.
$$ \begin{equation} \begin{aligned} \left\{ \begin{array}{lc} \frac{dy}{dt}=\frac{\alpha_y}{1+\frac{K_1}{1+K_{L}L}}-\gamma_{y}y,\\ \frac{dy}{dt}=\frac{\alpha_y}{1+K_2(1+K_{O_{2}}O_{2})}-\gamma_{y}y,\\ \frac{dy}{dt}=\frac{\alpha_y}{(1+K_3)H^{+}+K_{H}}-\gamma_{y}y,\\ y(0) = 0 \end{array} \right. \end{aligned} \end{equation} $$where: $\alpha_y$ is maximum transcription rate/maximum production rate;
$\gamma_y$ is total degradation/dilution rate of RFP;
$ K_1 $ is Binding affinity of transcription factor LIdR to promoter plldP;
$ K_L $ is Binding affinity of transcription factor LIdR to lactic acid;
$ K_2 $ is Binding affinity of transcription factor FNR to promoter pTeT;
$ K_{O_2} $ is affinity of transcription factor FNR for oxygen binding;
$ K_3 $ is Binding affinity of transcription factor CadC to promoter pCadC;
$ K_H $ is affinity of the transcription factor CadC for hydrogen ion binding.
Then, We use this equation to describe the dynamic expression of RFP(y) in different biosensors.
ODE model II: Expression of lysis gene
In this model, a new lysis gene was added to the gene fragment of the biosensor. Therefore, we hope to establish a model to describe and predict the changes of biosensor population with time after the addition of lysis genes.
We denote $r$ as the growth rate of CR flora, $N$ as the number of CR flora. If we assume that $r$ doesn't change over time, then the simplest model is that:
$$\begin{equation}\frac{dN}{dt} = rN\end{equation}$$In practice, for the limit of nutrients and living space in the bacterial culture environment, there will be a maximum upper limit number of flora $N_{max}$, and at the same time, as the number of flora increases, the growth rate of the flora should decrease, so $r$ is related to $N$. Under the assumption that the natural growth rate is $r_0$ we can obtain:
$$\begin{equation}r(N) = aN+r_0\end{equation}$$When $ N = N_{max} $, namely, the number of flora reaches its maximum, at which point the growth rate of the flora should be 0, so it can be obtained:
$$\begin{equation}r(N_{max})=aN_{max}+r_0=0\end{equation}$$Therefore, we have:
$$\begin{equation}a=-\frac{r_0}{N_{max}}\end{equation}$$Thus, a logistics equation that describe the growth of the flora can be obtained[8]:
$$\begin{equation}\frac{dN}{dt}=r_0\left(1-\frac{N}{N_{max}}\right)N\end{equation}$$The above model is a description and prediction of the number of bacterial populations in natural environments over time. However, this equation doesn’t apply to the biosensor population in this experiment . Therefore, we try to modify this equation.
The time-dependent relationship of the translation product produced by the biosensor is shown in the formula 10. Therefore, we can obtain the following model for the lactic acid biosensor with lysis gene:
$$ \begin{equation} \begin{cases} \displaystyle \frac{dN}{dt} &= r_0\left(1-\frac{N}{N_{max}}\right)N-\frac{\alpha_y}{1+\frac{K_1}{1+K_{L}L}}N\\ N(0) &= N_0\\ \end{cases} \end{equation} $$Similarly, we can obtain mathematical models describing the pH biosensor and hypoxia biosensor with lysis genes:
$$ \begin{equation} \begin{cases} \displaystyle \frac{dN}{dt} &= r_0\left(1-\frac{N}{N_{max}}\right)N-\frac{\alpha_y}{(1+K_3)H^{+}+K_{H}}N\\ N(0) &= N_0\\ \end{cases} \end{equation} $$ $$ \begin{equation} \begin{cases} \displaystyle \frac{dN}{dt} &= r_0\left(1-\frac{N}{N_{max}}\right)N-\frac{\alpha_y}{1+K_2(1+K_{O_{2}}O_{2})}N\\ N(0) &= N_0\\ \end{cases} \end{equation} $$Then, we use this equation to describe the dynamic of number in different biosensors.
Model Simulation
Determined parameter values
Model I: Expression of RFP gene
For the above biosensors, we first constructed lactic acid biosensors with pLldR as their promoters, and detected curves of RFP concentration over time under the experimental conditions of lactic acid concentration of 0Mm, 0.5mm, 1mM, and 1.5mm, respectively. Three groups of parallel experiments were performed for each concentration gradient.
For lactate biosensor, the model is as follows:
$$ \begin{equation} \frac{dy}{dt} = \frac{\alpha_y}{1+\left(\frac{K_1}{1+K_LL}\right)}-\gamma_y y, \quad y(0) = 0 \end{equation} $$If we denote $y_j^i $ as the experimental data of the $ i $ th parallel group at $ j $th sampling time, and $ N $is the number of data samples in each parallel group, then the data sampling time can be denoted as $\{t_j\}_{j=1}^N $. Since our experiment was conducted in triplicate control group, $i \in\{1,2,3\}$,.
Our goal is to select the above parameters reasonably so as to minimize the error between the curve of the formula 14 and the experimental data. Geometrically, our model curve needs to be the closest to the experimental data curve, mathematically, we need to find the parameter $\alpha_y^*, K_1^*, K_L^*, \gamma_y^* $that minimizes the error between the model and the experimental data and satisfies the following conditions:
$$ \begin{equation} \min \sum_{j=1}^{N} \left(\tilde{y_j}-\bar{y_j}\right)^2, \quad \bar{y_j} = \frac{\sum_{i=1}^3y_j^i}{3}\end{equation} $$ $$\begin{equation} \left\{ \begin{aligned} \displaystyle \frac{dy}{dt} &= \displaystyle \frac{\alpha_y^*}{1+\left(\frac{K_1^*}{1+K_L^*L}\right)}-\gamma_y^* y \\ \tilde{y_j} &= y(t_j) \\ \end{aligned} \right. \end{equation} $$The objective function is in fact called least square method (LSM). The essence of this problem is an optimization problem. So, for the target function: $ \sum_{j=1}^{N} \left(\tilde{y_j}-\bar{y_j}\right)^2 $, we apply the particle swarm optimization algorithm(PSO)[9],To find the global optimal parameter $\alpha_y^*, K_1^*, K_L^*, \gamma_y^* $, the following is the iterative solution procedure of the particle swarm optimization algorithm.
Its basic algorithm is as follows: Let $S$ be the number of particles in the swarm, each having a position $x_i \in R^n$ in the search-space and a velocity $v_i \in R^n$. Let $p_i$ be the best known position of particle $ i $ and let $ g $ be the best known position of the entire swarm. A basic PSO algorithm is then:
The values $b_{lo}$ and $b_{up}$ represent the lower and upper boundaries of the search-space respectively. The $w$ parameter is the inertia weight. The parameters $\phi_p$ and $\phi_g$ are often called cognitive coefficient and social coefficient. The termination criterion can be the number of iterations performed, or a solution where the adequate objective function value is found. The parameters $w$, $\phi_p$, and $\phi_g$ are selected by the practitioner and control the behavior and efficacy of the PSO method.
From the above animation, we can see that after nearly 70 iterations, our algorithm finally finds the parameters $\alpha_y^*, K_1^*, K_L^*, \gamma_y^* $ that minimize the error between the model and the experimental data, which is the curve of the model that is closest to the experimental data curve.
After that, we selected the average RFP concentration of three parallel groups as the data point, plotted the test data curve that generated the change of RFP concentration over time, and plotted the model curve and the test data curve at the same time, as shown in the following figure:
Then, we constructed pH biosensors with pCadC as their promoters, and detected curves of RFP concentration over time under the experimental conditions of $H^+$ concentration in 10^{-5.3}, 10^{-5.8}, 10^{-6.3}, 10^{-7.3} respectively. Three groups of parallel experiments were performed for each concentration gradient.
Similarly, we use the above method and PSO to find best parameters and plot time curves of the models.
Finally, we constructed hypoxia biosensors with pPepT as their promoter, and detected curves of RFP concentration over time under the experimental conditions of $ O_2 $ volume fraction of 0%, 20%, respectively. Three groups of parallel experiments were performed for each volume fraction gradient.
Similarly, we use the above method and PSO to find best parameters and plot time curves of the models.
Model II: Expression of lysis gene
For the above biosensors, we first constructed lactic acid biosensors with pLldR as their promoters, and connected the red fluorescent protein gene and the φ174E lysis gene, then detected curves of $OD_{600}$ over time under the experimental conditions of lactic acid concentration in 0Mm, 0.5mm, 1mM, and 1.5mm, respectively.
Three groups of parallel experiments were performed for each concentration gradient.
For lactate biosensor, the model is as follows:
$$ \begin{equation} \begin{cases} \displaystyle \frac{dN}{dt} &= r_0\left(1-\frac{N}{N_{max}}\right)N-\frac{\alpha_y}{1+\frac{K_1}{1+K_{L}L}}N\\ N(0) &= N_0\\ \end{cases} \end{equation} $$Similarly, we still use PSO and LSM to global optimization to get the optimal parameters, then, the obtained optimal model and the experimental data points after the average of the three experiments were plotted in the same coordinate system, and the curves were shown in the figure below.
In the figure above, the circles with different colors represent the data points obtained by experimental groups with different lactic acid concentrations, the dashed lines represent the curves of experimental groups with different lactic acid concentrations, and the solid lines with different colors represent the curves of the obtained model.
Then for pH biosensors with pCadC as their promoters, then we detected curves of $OD_{600}$ over time under the experimental conditions of $H^+$ concentration in 10^{-5.3}, 10^{-5.8}, 10^{-6.3}, 10^{-7.3} respectively. Three groups of parallel experiments were performed for each concentration gradient.
Finally, we constructed hypoxia biosensors with pPepT as their promoter, and detected curves of $OD_{600}$ over time under the experimental conditions of $ O_2 $ volume fraction of 0%, 20%, respectively. Three groups of parallel experiments were performed for each volume fraction gradient.
The estimation of parameters can be found in the Appendix A. Parameters, models are built by MATLAB as their tool. Codes of these models can be found in Appendix B. Codes.
Sensitivity analysis
In this part, we analyze the sensitivity of the model I. Since the deviation between the model obtained by the lactate biosensor with the lactate concentration of 0.5mM and the actual data is the smallest, we use the model corresponding to the lactate biosensor with the lactate concentration of L=0.5mM for analysis.
According to Appendix A. Parameters, the optimal model Parameters obtained by us are:
| Parameter | $\alpha_y^*$ | $K_1^*$ | $K_L^*$ | $\alpha_y^*$ |
| Optimized Value | 9585.988858886 | 2919.519950676 | 0.1000000000000 | 0.1000000000000 |
Therefore, we first try to see how the model changes when $\alpha_y$varies in the range of $[90\%\alpha_y^*,110\%\alpha_y^*]$.
Firstly, we devide $[ 90 \ % \ alpha_y ^ *, 110 \ % \ alpha_y ^ *] $ in 20 parts uniformly, to give out 20 values of $\ alpha_y ^i $, where $i\ in\ {1 ,\cdots, 20 }$. Then we plot time curves of these differential equations.
$$\begin{equation} \frac{dy^i}{dt}=\frac{\alpha_y^i}{1+\frac{K_1^*}{1+K_{L}^*L}}-\gamma_{y}^*y^i, \quad i=1,\dots,20\\ \end{equation} $$After that, we plotted the resulting function curve in a unified coordinate system, with time as the X-axis, RFP concentration as the Y-axis, and $\alpha_y^ i $value as the z-axis, as shown below:
Below is a plot of the relationship between the experimental data points and the changing function curve. The different colors of the data points indicate the different number of groups in the experimental group when they were collected. The black solid line represents the best parameter curve obtained by PSO, and the different colored dotted lines represent the curve after the deviation of the parameters
After that, we select 200 time points $t_j$ equally spaced, where $j\in\{1,\dots,200\}$, to calculate the maximum deviation ratio $\eta_j$ of the curve corresponding to $t_j$ from the optimal curve at different times for different values of $\alpha_y^ i $, which can be calculated as follows:
$$\begin{equation} \eta_j = \max_{i\in\{i,\dots,20\}}\left(\frac{y^i(t_j)-y(t_j)}{y(t_j)}\right)^2 \times 100\% \end{equation} $$Finally, we plot $(t_j,\eta_j)$ with smooth curves to obtain the curve of the maximum deviation ratio curve over time.
As can be seen in the figure above, the maximum deviation ratio $\eta_j$ of our curve increases monotonically with time. The maximum deviation ratio of $\eta_j$ remains within 1% even when the parameters change within 10%, which shows that our model has high stability.
We also carry out a similar sensitivity analysis for $K_1^*$, $K_{L}^*$, $\gamma_{y}^*$ at the same time, the variation range of the curve still remains in the range of $[90\%,110\%]$, therefore, we can draw the differential equation curve as shown in the following figure.
Finally, we use $\ alpha_y ^ * $, $K_1 ^ * $, $K_ {L} ^ * $, $\ gamma_ { y ^ *} $ respectively to draw the corresponding curve of optimal curve of maximum deviation and plot these under the same axis, the result is as follows:
From the above figure, we can see that the maximum deviation ratio of our curve $\eta_j$ remains within 1.5% even when the parameters are varied by 10%, indicating that our model is stable.
Similarly, sensitivity analysis was also carried out for other concentrations of lactic acid biosensors and two other types of biosensors, which again proved that our model had high stability.
Reference
[1] Guldberg, C.M.; Waage, P. (1879). "Ueber die chemische Affinität" [On chemical affinity]. Journal für praktische Chemie. 2nd series (in German). 19: 69–114. doi:10.1002/prac.18790190111. Reprinted, with comments by Abegg in Ostwalds Klassiker der Exacten Wissenschaften, no. 104, Wilhelm Engleman, Leipzig, 1899, pp 126-171
[2] Chien, Tiffany, et al. "Enhancing the tropism of bacteria via genetically programmed biosensors." Nature biomedical engineering 6.1 (2022): 94-104.
[3] Weghoff, M. C., Bertsch, J. & Muller, V. A novel mode of lactate metabolism in strictly anaerobic bacteria. Environ. Microbiol. 17, 670–677 (2015).
[4] Stirling, F. et al. Synthetic cassettes for pH-mediated sensing, counting, and containment. Cell Rep. 30, 3139–3148 (2020).
[5] Crack, J. C. et al. Signal perception by FNR: the role of the iron-sulfur cluster. Biochem. Soc. Trans. 36, 1144–1148 (2008).
[6] Alon, U. An Introduction to Systems Biology: Design Principles of Biological Circuits 3–19 (Chapman & Hall/CRC, 2007).
[7] Michaelis–Menten Equation. In: Encyclopedia of Genetics, Genomics, Proteomics and Informatics. Springer, Dordrecht. https://doi.org/10.1007/978-1-4020-6754-9_10317 (2008).
[8]Gang, X. U. , H. Wen , and W. U. Kun . "Primal chaos data system and feedback control research of Logistic population increase model." Journal of Natural Science of Heilongjiang University (2003).
[9] Mohammad Reza Bonyadi, Zbigniew Michalewicz; Particle Swarm Optimization for Single Objective Continuous Space Problems: A Review. Evol Comput 2017; 25 (1): 1–54. doi: https://doi.org/10.1162/EVCO_r_00180