Stability analysis of genetic regulatory network with additive noises

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][2][3][4][5][6][7][8][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 pro-teins 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][10][11][12][13][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 mRNAs 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][11][12][13][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.
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][17][18][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 where m i (t) and p j (t) are the concentrations of mRNAs 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 mRNAs and proteins, k dji are assumed to be constants describing the link between mRNAs and proteins. If the jth protein regulates the ith gene as illustrated by the dashed arrow in figure 1, there is a regulatory link A gene regulatory process.

Intracellular molecular
in the network. Such regulation effects of proteins to gene expressions are described by wrapped nonlinear terms b i (p 1 (t), p 2 (t), … , p m (t)), which are nonlinear functions of p 1 (t), p 2 (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).
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) where Vector is determined by within Ω ij , which is a set of repressors for gene i. Matrix describes the interconnection of transcription factors and genes in the network with element α ij , and 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.
• g(p j ) satisfies a sector condition, which can be obtained Combine the state m(t), p(t) together, the state space model for control design and stability analysis is written as follows: where

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.
γ , . and 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 where is the desired transcription rates. Stability of the genetic network means: (1) State x * is the equilibrium of equation (7), i.e.
(2) Starting at any initial states x 0 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)  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: 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][24][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 where i = 1, 2, … , n, j = 1, 2, … , m, e k is the kth 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.  with the adaption control chosen as From the above Lyapunov argument, with a positive definite V and negative definite 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 where are piecewise continuous in t, the systems can be viewed as the nominal system 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. 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 satisfies the following condition: With the additive noise disturbance (t, Φ) considered, the time derivative of V 2 is rewritten as , :

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].
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.

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 col-Control parameters with measurement noises. 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.