Stability analysis of genetic regulatory network with additive noises

Abstract

Background

Genetic regulatory networks (GRN) can be described by differential equations with SUM logic which has been found in many natural systems. Identification of the network components and transcriptional rates are critical to the output behavior of the system. Though transcriptional rates cannot be measured in vivo, biologists have shown that they are alterable through artificial factors in vitro.

Results

This study presents the theoretical research work on a novel nonlinear control and stability analysis of genetic regulatory networks. The proposed control scheme can drive the genetic regulatory network to desired levels by adjusting transcriptional rates. Asymptotic stability proof is conducted with Lyapunov argument for both noise-free and additive noises cases. Computer simulation results show the effectiveness of the control design and robustness of the regulation scheme with additive noises.

Conclusions

With the knowledge of interaction between transcriptional factors and gene products, the research results can be applied in the design of model-based experiments to regulate gene expression profiles.

Background

Genetic networks regulate sophisticated biological functions by interacting genes and proteins and support homeostasis in metabolism and coordinate events during the developmental program. Research on stability analysis and regulation/control of these genetic networks are particularly important. Pioneer experimental studies in construction of genetic networks to manipulate protein levels or even to construct gene circuits with repressor functions have been carried out [1–9]. These experiments have demonstrated interesting properties of GRNs of E. coli in the presence of specified repressors. With different repressors, these GRNs include single or multiple interactions between genes and proteins. In a single gene regulatory network [1], the negative feedback that is integrated in the system decreases cell-to-cell fluctuations in protein concentration measurements. Distribution of the regulated protein concentration is proportional to the degradation rate of the gene network. In a two-gene regulation network [3], bi-stability is shown by coupling two proteins with negative regulation in the synthesis of each other. Stability analysis of this bi-gene regulation network is also presented based on parameter bifurcations. The significance of this experiment is that the transition between two stable states of the GRN is much sharper with respect to the intracellular stimuli, i.e., performance can be adjusted by special input. In a tri-gene regulatory model [2], three genes are regulated with cyclic repressibility. The system exhibited self-sustained oscillations over the entire growth phase of the host E. coli cells for certain biochemical parameters. These biological experimental results have shown that genetic networks can be regulated by a scheme of local promotor control, i.e. the number, type and placement of regulatory protein binding sites. However, quantitative analysis of the regulation function has not been studied.

The objective of our current study is to develop a mathematical model of the tri-gene regulation network and extend the theoretical stability analysis to the case with measurement noises. A novel control scheme is proposed to change the state of a genetic regulatory network by adjusting transcriptional rates. The proposed control scheme can provide biologists useful design techniques in model-based experiments to predict protein levels in genetic regulatory networks by adjusting specific regulatory factors. We will use the biological scheme of adjusting regulatory factors reported in research articles [9–14]. These regulatory factors will be adjusted based on errors between the measured gene products, mRNA and protein, concentration levels and their desired values to regulate the gene expressions profiles. The proposed control is based on the concentrations of m RNAs and proteins which can be easily measured by current molecular biology techniques. Therefore, it's useful when some of the transcription rates are unmeasurable.

We will introduce the general model of GRNs to explain the work, present the control design and the stability analysis , and then give a simulation example and conclusion.

Mathematical model of genetic regulatory networks

A gene regulatory network for a eukaryote is shown in figure 1. Potential inputs of the system include a variety of developmental and environmental stimuli. System outputs are the synthesized proteins. Inputs activate a complex chain of intracellular reactions that activate a regulatory molecule, transcriptional factor, to translocate from the cytoplasm into the nucleus. In the nucleus, active transcription factors recognize a specific segment of DNA, termed a promotor. The promotor informs RNA polymerase, which binds closely to the promotor, where to start transcribing genetic information on DNA to mRNA. The mRNA molecules then leave the cell nucleus and enter the cytoplasm where proteins are synthesized in the presence of transfer RNAs. When the translated protein is capable of interacting with its own or other gene's transcription factor (denoted by the dashed arrow in figure 1, a regulatory or feedback loop is formed. Such transcriptional regulation is the typical method utilized by cells to control gene expression [9, 15]. Feedback can occur in either a positive (activator) or a negative (repressor) direction. In a gene regulatory network, natural regulatory molecules are either activated transcription factors or proteins that activate transcription factors. Gene expression is very sensitive to changes in the composition of regulatory factors, which limits attempts to control gene expression. If artificial regulatory factors [10–14] can be used to change the activity of the transcription factors and the affinity of DNA binding, we can affect the transcriptional rate in GRNs and control the output behavior of the GRNs. Thus, in order to derive novel agents that regulate gene expressions, we must first understand GRN output behavior.

Figure 1

A gene regulatory process. A gene regulatory process includes input (chemical signals), output(proteins) and feedback (regulatory) factors.

Two types of mathematical models have been developed to understand the working scheme of genetic networks: (1) Boolean networks and (2) sets of differential equations [9, 16–19]. Boolean network model describes the expression of each gene with two states: ON or OFF, and the state of a gene is determined by a Boolean function of the states of other related genes. The differential equation model describes concentrations of mRNAs and proteins as continuous functions, which can provide more accurate and detailed dynamic information. Since genetic networks are high dimensional and nonlinear in nature, it is logical to consider such genetic networks from a nonlinear dynamic viewpoint, i.e., nonlinear differential equations. From a control point of view, main purpose of our mathematical model is to predict and manipulate the dynamic behaviors by analyzing available measurements. Researches of control design for Boolean networks have been carried out [20, 21], while for continuous differential model, there are comparatively few literature references available [22].

In this section, the studied genetic network model is described by nonlinear differential equations

$$\begin{array}{l}{\dot{m}}_i(t)=-k_{mi}{m}_i(t)+b_i(p_1(t),p_2(t),\dots p_m(t))\\ {\dot{p}}_j(t)=-k_{pj}{p}_i(t)+k_{dji}{m}_i(t),\end{array}$$

(1)

where m _i (t) and p _j (t) are the concentrations of m RNAs and proteins in a genetic regulatory network, i = 1, 2, … , n, and j = 1, 2, … , m. Parameters k _mi and k _pj are degradation rates of the m RNAs and proteins, k _dji are assumed to be constants describing the link between m RNAs and proteins. If the jth protein regulates the ith gene as illustrated by the dashed arrow in figure 1, there is a regulatory link in the network. Such regulation effects of proteins to gene expressions are described by wrapped nonlinear terms b _i (p₁(t), p₂(t), … , p _m (t)), which are nonlinear functions of p₁(t), p₂(t), … , p _m (t), and can be described by SUM logic (additive effect) of each protein to the specified gene. Detailed description of the relationship between α _i s, b _i s and transcriptional factors have been previously published [9, 17, 18]. Adjustment of the regulatory factors will change parameters α _ij in the mathematical depiction written in equation (2).

$$b_{ij}(p_j(t))=\{\begin{array}{l}{\alpha }_{ij}\frac{{(p_j(t)/\beta )}^H}{1+{(p_j(t)/\beta )}^H}\text{ case 1 }\\ {\alpha }_{ij}-{\alpha }_{ij}\frac{{(p_j(t)/\beta )}^H}{1+{(p_j(t)/\beta )}^H}\text{ case 2}\end{array}$$

(2)

where case 1 represents that transcription factor (protein) j is an activator of gene i, and in case 2, transcription factor j is a repressor of gene i. Parameters H in equation (2) is the Hill coefficient, β is a positive constant, and α _ij is the dimensionless transcription rate of transcription factor j to gene i. Equation (1) can be rewritten as the following form by applying (2)

$$\begin{array}{l}\dot{m}(t)=Am(t)+Bg(p(t))+l_0\\ \dot{p}(t)=Cp(t)+Dm(t)\end{array}$$

(3)

where $m(t)=[m_1(t),\dots m_i(t)\dots m_n(t)]\in {\Re }^n,p(t)=[p_1(t),\dots p_j(t)\dots p_m(t)]\in {\Re }^m,A=diag\{-k_{m_i}\}\in {\Re }^{n\times n},C=diag\{-k_{p_j}\}\in {\Re }^{m\times m},D=diag\left\{k_{dji}\right\}\in {\Re }^{m\times n}.$ Vector $l_0\in {\Re }^n$ is determined by

$$l_{0i}={\displaystyle \sum _{{\Omega }_{ij}}{\alpha }_{ij}},$$

(4)

within Ω _ij , which is a set of repressors for gene i. Matrix $B\in {\Re }^{n\times m}$ describes the interconnection of transcription factors and genes in the network with element α _ij ,

$$B_{ij}=\{\begin{array}{ll}0\hfill & \text{no regulation for protein }j\text{ to gene }i,\hfill \\ {\alpha }_{ij}\hfill & \text{case 1,}\hfill \\ -{\alpha }_{ij}\hfill & \text{case 2,}\hfill \end{array}$$

and $g\left(p\left(t\right)\right)={\left[g\left(p_1\left(t\right)\right)\dots g\left(p_j\left(t\right)\right)\dots g\left(p_n\left(t\right)\right)\right]}^T\in {\Re }^m,g\left(p_j\left(t\right)\right)=\frac{{\left(p_j\left(t\right)/\beta \right)}^H}{1+{\left(p_j\left(t\right)/\beta \right)}^H}.$ Some special properties of the nonlinear function g(p(t)) should be pointed out:

• g(p _j ) > 0 always holds with p _j > 0 and the equal sign holds only when p _j =0;

• g(p _j ) is a monotonic increasing function, i.e. $\frac{{\partial }_g(p_j)}{\partial p_j}\ge 0;$

• g(p _j ) satisfies a sector condition, which can be obtained from the mean-value theorem, $g\left(e_j+p_j\right)-g\left(p_j\right)=\frac{\partial g}{\partial p_j}|p_j=\xi \cdot e_j,\xi \in \left[p_j,p_j+e_j\right],$ and thus, $0\le e_j[g\left(e_j+p_j\right)-g\left(p_j\right)]\le {\gamma }_j{e}_j^2,$ where ${\gamma }_j=\max \left\{\frac{H{\left(p_j/\beta \right)}^{H-1}}{\beta {\left(1+{\left(p_j/\beta \right)}^H\right)}^2}|p_j\in \left[p_j,p_j+e\right]\right\}.$ The vector case of the above equation can be expressed as

  $$0\le e^T\left[g\left(e+p\right)-g\left(p\right)\right]\le \gamma e^{T}e,$$

  (5)

where γ = max {γ₁, γ _j …, γ _m }.

Combine the state m(t), p(t) together, the state space model for control design and stability analysis is written as follows:

$$\begin{array}{l}\dot{x}=f\left(t,x\right)=A_{c}x+B_{c}u+l,u=-g\left(y\right)\\ y=C_{c}x,\end{array}$$

(6)

where $x^T=\left[m^T{p}^T\right],l^T=\left[l_0^{T}0\right],x,l\in {\Re }^{n+m},C_c=[0_{m\times n}{I}_{m\times m}]\in {\Re }^{m\times (n+m)}$

$$A_c=\left[\begin{array}{cc}A& 0\\ D& C\end{array}\right]\in {\Re }^{\left(n+m\right)\times \left(n+m\right)},\text{and}{B}_c=\left[\begin{array}{c}-B\\ 0\end{array}\right]\in {\Re }^{\left(n+m\right)\times n}.$$

Methods: adaptive control of genetic networks in a noise-free case

Nonlinear adaptive control has been applied to many systems to improve performance. The control objective here is to drive the current state of genetic networks to desired values m^* and p^*. In order to make the control design biologically meaningful, the following assumptions are introduced.

Assumption 1: Transcription rates of mRNAs in the studied genetic network are adjustable.

This assumption guarantees the possibility of adjusting the transcription rates and drive the current states to the desired values.

Assumption 2: State x^*T = [m^*T p^*T] is the stable steady state generated by the same genetic network in equation (3) with desired transcription rates.

This assumption guarantees that m^* and p^* are achievable and has biological meaning for a real genetic network. The stability also guarantees that once the state is driven to x^*, it will stay there.

The control design proposed here includes two parts: i) stability analysis of x^* generated by the genetic network with the desired degradation, transcription rate, the given sector parameter γ; ii) control design to drive the current states to the desired x^*. Assume x^* is generated by

$$\begin{array}{l}\dot{x}=A_{c}x+B_c^{*}u+l,u=-g\left(y\right),\\ y=C_{c}x.\end{array}$$

(7)

where $B_c^{*}$ is the desired transcription rates. Stability of the genetic network means: (1) State x^* is the equilibrium of equation (7), i.e. $A_c{x}^{*}-B_c^{*}g\left(C_c{x}^{*}\right)+l=0;$ (2) Starting at any initial states x₀ close to x^*, the trajectory will converge to x^* as time goes to infinity. If we define e^* = x − x^*, the error e → 0 as t → ∞.

Lemma 1: The system (7) is globally exponentially stable if the transfer function matrix G(s) = C _c (sI − A _C )⁻¹$B_c^{*}$is strictly positive real, i.e. with the controllable and observable pair$\{A_c,B_c^{*},C_c\}$, there exist positive matrices P and Q, P = P^T, such that

$$P{A}_c+A_c^{T}P=-Q,P{B}_c^{*}=C_c^T.$$

(8)

Proof: With constant x^*, the error dynamics can be expressed as

$$\begin{array}{c}\dot{e}=A_{c}x+B_c^{*}u+l=A_c\left(e+x^{*}\right)-B_c^{*}g\left[C_c\left(e+x^{*}\right)\right]+l\\ =A_{c}e-B_c^{*}\left[g\left(C_{c}e+C_c{x}^{*}\right)-g\left(C_c{x}^{*}\right)\right].\end{array}$$

(9)

Choosing the Lyapunov function $V_1(e)=\frac{1}{2}{e}^{T}Pe,$ its time derivative ${\dot{V}}_1(e)$ is obtained by the following equations.

$$\begin{array}{c}{\dot{V}}_1=\frac{1}{2}\left({\dot{e}}^{T}Pe+e^{T}P\dot{e}\right)\\ =\frac{1}{2}\left\{e^T{A}_c^T-{\left[g\left(C_{c}e+C_c{x}^{*}\right)-g\left(C_c{x}^{*}\right)\right]}^T{B}_c^{*T}\right\}Pe\\ +\frac{1}{2}{e}^{T}p\left\{A_{c}e-B_c^{*}\left[g\left(C_{c}e+C_c{x}^{*}\right)-g\left(C_c{x}^{*}\right)\right]\right\}\\ =\frac{1}{2}{e}^T\left[P{A}_c+A_c^{T}P\right]e-e^{T}P{B}_c^{*}\left\{g\left[C_{c}e+C_c{x}^{*}\right]-g\left(C_c{x}^{*}\right)\right\}\\ =-\frac{1}{2}{e}^{T}Qe-{\left(C_{c}e\right)}^T\left\{g\left[C_{c}e+C_c{x}^{*}\right]-g\left(C_c{x}^{*}\right)\right\}\le -\frac{1}{2}{e}^{T}Qe\end{array}$$

(10)

which is negative definite since (C _c e)^T {g[C _c e + C _c x^*)]—g(C _c x^*)} > 0 by applying equation (5). Thus, the error dynamics will converge to zero exponentially by satisfying the strict positive real condition given as equation (8) on the system parameters.

Remark 1: Lemma 1 gives a sufficient condition on parameter settings for a stable genetic network. The stability is determined by system parameters and the sector condition of the nonlinear feedback function.

Lemma 1 provides an easy way to check the stability of a genetic network, since it is easy to get the transfer function for linear time invariant system and check whether the transfer function is strictly positive real or not.

Remark 2: While applying the small gain theorem, with the consideration of the sector parameter γ, we can get the sufficient and necessary condition: ${\sup }_{\omega }\left|G\left(j\omega \right)\right|<1/\gamma ,\omega \in \Re .$ Since the necessary condition is not involved in the following control design, we ignore the proof here. Related information of it can be found in [23–25].

Based on the above stability analysis, the following theorem gives the control design that drives the current state to the desired state by adjusting the transcription rates of the system.

Theorem 1: Assume transcription rates of gene i in system (7) can be adjusted by the control law

$${\dot{\alpha }}_{ij}=-k_{ij}{\displaystyle \sum _{i,j,k=1}^{i=n,j=m,k=n+m}{e}_k{P}_{ki}{g}_j,}$$

(11)

where i = 1, 2, … , n, j = 1, 2, … , m, e _k is the k th element in the error vector between the current state and the desired state, i.e., e = x—x^*, g _j = g(p _j (t)), P _ij is the element in the positive definite matrix P defined in equation (8). The system (6) will converge to the desired state x^* asymptotically as time goes to infinity.

Proof: From equations (6, 7), the error dynamics is

$$\begin{array}{c}\dot{e}=A_{c}x+B_{c}u+l=A_{c}e-B_{c}g\left(y\right)+B_c^{*}g\left(y*\right)\\ =A_{c}e-\left(B_c-B_c^{*}\right)g\left(y\right)-B_c^{*}\left[g\left(y\right)-g\left(y*\right)\right]\end{array}$$

(12)

Choose Lyapunov candidate $V\left(e,{\tilde{\alpha }}_{ij}\right)=\frac{1}{2}{e}^{T}Pe+{\displaystyle {\sum }_{i,j=1}^{i=n,j=m}\frac{1}{2k_{ij}}{\tilde{\alpha }}_{ij}^2},$ where ${\tilde{\alpha }}_{ij}\text{s}$ are defined as ${\tilde{\alpha }}_{ij}=B_c-B_c^{*}={\alpha }_{ij}-{\alpha }_{ij}^{*}.$ Recall that the bottom block with m × n dimension in both matrices are zeroes, ${\tilde{\alpha }}_{ij}$ is only related to the transcription rates in the system. Thus, the time derivative of the Lyapunov function can be derived as follows.

$$\begin{array}{c}\dot{V}=\frac{1}{2}\left({\dot{e}}^{T}Pe+e^{T}P\dot{e}\right)+{\displaystyle \sum _{i,j=1}^{i=n,j=m}\frac{1}{k_{ij}}{\tilde{\alpha }}_{ij}{\dot{\tilde{\alpha }}}_{ij}}\\ =\frac{1}{2}\left\{e^T{A}_c^T-g{\left(y\right)}^T{\left(B_c-B_c^{*}\right)}^T-{\left[g\left(y\right)-g\left(y*\right)\right]}^T{B}_c^{*T}\right\}Pe\\ +\frac{1}{2}{e}^{T}P\left\{A_{c}e-\left(B_c-B_c^{*}\right)g\left(y\right)-B_c^{*}\left[g\left(y\right)-g\left(y*\right)\right]\right\}+{\displaystyle \sum _{i,j=1}^{i=n,j=m}\frac{1}{k_{ij}}{\tilde{\alpha }}_{ij}\dot{\tilde{\alpha }}}\\ =\frac{1}{2}{e}^T\left[P{A}_c+A_c^{T}P\right]e-e^{T}P{B}_c^{*}\left\{g\left[C_{c}e+C_c{x}^{*}\right]-g\left(C_c{x}^{*}\right)\right\}\\ -e^{T}P\left(B_c-B_c^{*}\right)g\left(y\right)+{\displaystyle \sum _{i,j=1}^{i=n,j=m}\frac{1}{2k_{ij}}{\tilde{\alpha }}_{ij}{\dot{\tilde{\alpha }}}_{ij}}\\ \le -\frac{1}{2}{e}^{T}Qe-{\displaystyle \sum _{i,j,k=1}^{i=n,j=m,k=n+m}{e}_k{P}_{ki}{\tilde{\alpha }}_{ij}{g}_j+{\displaystyle \sum _{i,j=1}^{i=n,j=m}\frac{1}{k_{ij}}{\tilde{\alpha }}_{ij}{\dot{\tilde{\alpha }}}_{ij}}}\\ \le -\frac{1}{2}{e}^{T}Qe\end{array}$$

(13)

with the adaption control chosen as ${\dot{\tilde{\alpha }}}_{ij}=-{\dot{\alpha }}_{ij}=k_{ij}{\displaystyle {\sum }_{i,j,k=1}^{i=n,j=m,k=n+m}{e}_k{P}_{ki}{g}_j}.$

From the above Lyapunov argument, with a positive definite V and negative definite $\dot{V}$ all the errors decrease to zero. This concludes the proof that the tracking error of the genetic network from current state to the desired state is globally asymptotically stable with the adaptive control given in (11), i.e., the states are driven to desired levels. As the tracking errors converge to zero, the adjustable transcription rates α _ij will converge to constants.

Boundedness with additive noises

When signals are sensed, signal distortion, transmission delay and noise are unavoidable. In this section an additive measurement noise will be considered, and the distorted measurements of mRNA and protein concentration levels will be used in the adaption law.

Considering the aforementioned system in equation (6) with additive noises

$$\dot{x}=f\left(t,x\right)+d\left(t,x\right),$$

(14)

where $f:{\Re }^{+}\times {\Re }^n\to {\Re }^n,d:{\Re }^{+}\times {\Re }^n\to {\Re }^n$ are piecewise continuous in t, the systems can be viewed as the nominal system

$$\dot{x}=f\left(t,x\right),$$

(15)

with perturbation term d(t, x). The variable u is a function of x and is omitted here for simplicity.

Assume the nominal system has an equilibrium point at the origin, if d(t, x) = 0 as x = 0, then the origin is still an equilibrium point with the disturbance d(t, x). However, in most cases, the disturbance does not satisfy this assumption, i.e. d(t, 0) ≠ 0. In this case, the origin is not an equilibrium point and no conclusion can be drawn about stability of the origin. The following Theorem shows that the best result we can expect is the uniform boundedness of the disturbed system when the origin is exponentially stable.

Theorem 2: Consider the system and adaptation laws given by equations (6), and (11), the system is uniformly bounded with disturbance $\Vert d\left(t,x\right)\Vert \le \delta$ where δ is a positive constant. The upper bound b _u and lower bound b _l of the system is determined by

$$b_u=\frac{{\alpha }_4\delta }{{\alpha }_3\mu }\sqrt{\frac{{\alpha }_2}{{\alpha }_1}},\text{and}{b}_l=\frac{{\alpha }_4\delta }{{\alpha }_3\mu },$$

where 0 < μ < 1, respectively. α _i s are positive constants defined as follows.

proof From the proof of Theorem 1, the error dynamics of the whole system is exponentially stable at the origin. Thus, there exist a Lyapunov function $V(\Phi ),\Phi ={[e,{\tilde{\alpha }}_{ij}]}^T$ satisfies the following condition:

$$\begin{array}{l}{\alpha }_1{\Vert \Phi \Vert }^2\le V\left(\Phi \right)\le {\alpha }_2{\Vert \Phi \Vert }^2,\\ \frac{\partial V}{\partial t}+\frac{\partial V}{\partial \Phi }f\left(t,\Phi \right)\le -{\alpha }_3{\Vert \Phi \Vert }^2,\\ \frac{\partial V}{\partial \Phi }\le {\alpha }_4\Vert \Phi \Vert .\end{array}$$

(16)

With the additive noise disturbance (t, Φ) considered, the time derivative of V₂ is rewritten as

$$\begin{array}{c}\dot{V}\le -{\alpha }_3{\Vert \Phi \Vert }^2+{\alpha }_4\Vert \Phi \Vert \cdot \Vert (t,\Phi )\Vert \\ \le -{\alpha }_3{\Vert \Phi \Vert }^2+{\alpha }_4\delta \Vert \Phi \Vert \\ \le -{\alpha }_3(1-\mu ){\Vert \Phi \Vert }^2-{\alpha }_3\mu {\Vert \Phi \Vert }^2+{\alpha }_4\delta \Vert \Phi \Vert \\ \le -{\alpha }_3(1-\mu ){\Vert \Phi \Vert }^2,\end{array}$$

(17)

for any $\Vert \Phi \Vert \ge \frac{\delta {\alpha }_4}{\mu {\alpha }_3}.$ This gives the lower bound as $bl=\frac{\delta {\alpha }_4}{\mu {\alpha }_3}.$ From BIBO Theorem from the reference [26], the upper bound can be determined as $bu=\frac{{\alpha }_4\delta }{{\alpha }_3\mu }\sqrt{\frac{{\alpha }_2}{{\alpha }_1}}.$

From theorem 2, it can be shown that with bounded noises, both the state estimation errors and adaptation control are bounded. With disturbance bound δ → 0, differences between the adaptation control and the desired constant transcriptional rates shrink to zero also.

Simulation results and discussion

To show the effectiveness of the proposed control, an example is simulated with Matlab. We consider the dynamics of the tri-gene network described in [2].

$$\begin{array}{c}{\dot{m}}_i(t)=-k_{mi}{m}_i+\frac{{\alpha }_i}{1+p_j^n}\\ {\dot{p}}_i=-k_{pi}{p}_i-k_{di}{m}_i\end{array}$$

(18)

where i sequence is lacl, tetR, cl and j sequence is cl, lacl, tetR, k _mi = 1, k _pi = 1, n = 2 and α is the control parameter. With parameters chosen as α _i = 2.5, the genetic network has a globally exponential equilibrium at $m_i^{*}=p_i^{*}=1.1150,$ which is the steady state. Evolution trajectories of m RNA and protein concentration levels are shown in figure 2.

Figure 2

Uncontrolled State of mRNA and protein levels. (A) Uncontrolled state of mRNAs of lacl, tetR, cl in the tri-gene regulatory in a noise-free case. The initial states of computer simulation start from initial status as [3, 2, 1]. (B) Uncontrolled proteins cl, lacl, tetR in the tri-gene regulatory in a noise-free case. All states converge to 1.115 with in 50 second simulation time.

The control objective is to alter the current state x = [3, 2, 1, 1, 3, 2]^T to the desired states x^* by adjusting transcription rates α_i. Assume there is no regulation at the initial state, the control law is chosen to be ${\dot{\alpha }}_i=-10e_i{g}_j.$ Evolution trajectories of m _i and p _i with control design are shown in figure 3. The control variables are shown in figure 4. It can be seen that the steady states converge to x^* and the transcription rates converge to the desired value 2.5 in the simulation. The simulation also shows that the convergence rate of the state and the parameter with the proposed control is faster than the original genetic network.

Figure 3

Controlled state of mRNA and protein levels in a noise-free case. (A) Controlled state of mRNAs of lacl, tetR, cl in the tri-gene regulatory in a noise-free case. The initial states of computer simulation start from initial status as [3, 2, 1]. (B) Controlled proteins cl, lacl, tetR in the tri-gene regulatory in a noise-free case. All states in (A) and (B) converge to 1.115 less than 20 second simulation time.

Figure 4

Control parameters in a noise-free case. Online estimated control variables α _i . Initial values of α _i are all chosen as zeros without losing generality. Control parameters converge to 2.5 with in 20 second simulation time.

In order to show the robust of the control scheme, an additive white noise is injected to the measurements of the mRNA and protein concentration levels. Magnitude of the noises are bounded by 0.5. All parameter settings are the same as the ideal case illustrated above. Simulation results of the regulated mRNA and protein concentration levels are shown in figure 5 and the control of the transcription rates are shown in figure 6.

Figure 5

Controlled state of mRNA and protein levels with measurement noises. (A) Controlled state of mRNAs of lacl, tetR, cl in the tri-gene regulatory with measurement noises. The initial states of computer simulation start from initial status as [3, 2, 1]. (B) Controlled proteins cl, lacl, tetR in the tri-gene regulatory with measurement noises. The initial states of computer simulation start from initial status as [1, 3, 2]^T. All states in (A) and (B) converge to 1.115 with limited error.

Figure 6

Control parameters with measurement noises. Online estimated control variables α _i with measurement noises. Initial values of α _i are all chosen as zeros without losing generality. Control parameters approach to 2.5 with limited error.

Conclusions and future research

In the current study, a general regulatory network is presented and a new control scheme is proposed. This scheme controls the current states of the genetic network to desired values. We have obtained global convergence of the tracking error by online adjustment of the transcription rates. We have applied Lyapunov argument to the convergence analysis of the states and boundedness of the transcription rates. A tri-gene regulatory is simulated with the proposed algorithm. Effectiveness of the proposed controller design is verified by Matlab simulation for noise-free measurement and bounded noises.

In this research, interaction between the proteins and mRNAs are described by Hill functions with SUM logic. More detailed studies need to be carried out to determine the structure of a real GRN and extend the control scheme to more general cases. In addition, the estimation algorithms are developed with continuous measurers of the state. In real biological systems, such measures will be collected at specific time points. Estimation error caused by sampling effects should also be considered. Together, our results demonstrated that GRN could be regulated by artificially changing the transcription rates, to approach the desired gene expression levels.