A Unifying Mathematical Framework for Genetic Robustness,Environmental Robustness,Network Robustness and their Trade-off on Phenotype Robustness in Biological Networks Part I: Gene Regulatory Networks in Systems and Evolutionary Biology

Abstract

Robust stabilization and environmental disturbance attenuation are ubiquitous systematic properties observed in biological systems at different levels. The underlying principles for robust stabilization and environmental disturbance attenuation are universal to both complex biological systems and sophisticated engineering systems. In many biological networks, network robustness should be enough to confer intrinsic robustness in order to tolerate intrinsic parameter fluctuations, genetic robustness for buffering genetic variations, and environmental robustness for resisting environmental disturbances. With this, the phenotypic stability of biological network can be maintained, thus guaranteeing phenotype robustness. This paper presents a survey on biological systems and then develops a unifying mathematical framework for investigating the principles of both robust stabilization and environmental disturbance attenuation in systems and evolutionary biology. Further, from the unifying mathematical framework, it was discovered that the phenotype robustness criterion for biological networks at different levels relies upon intrinsic robustness + genetic robustness + environmental robustness ≤ network robustness. When this is true, the phenotype robustness can be maintained in spite of intrinsic parameter fluctuations, genetic variations, and environmental disturbances. Therefore, the trade-offs between intrinsic robustness, genetic robustness, environmental robustness, and network robustness in systems and evolutionary biology can also be investigated through their corresponding phenotype robustness criterion from the systematic point of view.

Introduction

Inherently, real physical and biological systems suffer from intrinsic perturbations and extrinsic disturbances. In the last two decades, robust stabilization and noise-filtering theories have been developed by control engineers to achieve the robust stability for tolerating intrinsic perturbations as well as to obtain noise-filtering ability against extrinsic disturbances for improving the reliability and performance of control systems.^1–5 Their applications are also extended from control systems^1–4 to signal processing systems, ⁵ communication systems, ⁶ and biological systems.^7,8 Since both engineering systems and biological systems need robust stabilization and noise-filtering abilities to tolerate intrinsic perturbations and resist extrinsic disturbances (or noises) so as to maintain their desired function or performance, there should exist some common schemes for robust stabilization and disturbance attenuation for these two kinds of systems.

At the molecular level, a gene regulatory network or protein interaction network is inherent with intrinsic parameter fluctuations due to random molecular fluctuations (or gene expression noises) and environmental disturbances. Since the study and design of gene regulatory networks and protein interaction networks have become important topics in systems biology, synthetic biology, and evolutionary biology,^9–13 the robust stabilization and noise-filtering properties of biochemical or genetic regulatory networks have attracted much attention of engineers and molecular biologists. Robust stabilization and noise-filtering ability of gene networks under intrinsic parameter fluctuations and environmental disturbances have recently been discussed from the nonlinear stochastic stabilization and noise-filtering point of view. Robustness and evolvability are found to be internal properties of biological systems. ¹² They determine a biological system's persistence and potentiality for future evolutionary changes. The biological system is evolvable if mutations in it are able to produce heritable phenotypic variations. In general, the more robust a system is, the less phenotypic variation a given number of mutations can generate, and hence the less evolvable the system is.^9,10

Further, the interplay between robustness and sensitivity in the gene regulatory network is also discussed from the nonlinear stochastic system point of view. It was found that if the sum of genetic robustness and environmental robustness is less than network robustness, ie, network robustness can confer both genetic robustness to tolerate genetic variations and environmental robustness to resist environmental disturbances, and then the phenotype of the gene regulatory network is robust under genetic variations and environmental disturbances. The trade-off between genetic robustness and environmental robustness in evolution is discussed from the viewpoint of stochastic stability robustness and noise sensitivity of the nonlinear stochastic gene network. This may be relevant to the statistical trade-off between bias and variance, the so-called bias/variance dilemma. ¹¹ Further the trade-off could be considered as an antagonistic pleiotropic action of a gene regulatory network from the systems biology perspective.

In the real evolutionary process of a population of biological networks, the random transmission and mutation of genes provide biological diversities for natural selection. In order to preserve functional phenotypes between generations, gene networks need to evolve robustly under discontinuous genetic mutations and environmental disturbances. In this study, a population of evolutionary gene networks is represented by nonlinear stochastic dynamic system with random genetic variations accumulated in evolution. Therefore, the robust stabilization of the natural favorite phenotype exerts a selection force on a population of gene networks to maintain network function. However, gene networks in population are generally adjusted by random genetic variations to generate phenotypes for new challenges in the network's evolution, ie, the evolvability. Hence, there should be some interplay between network evolvability and robust stabilization in evolutionary gene networks. The interplay between network evolvability and network robustness of biological networks has been discussed from the nonlinear stochastic point of view. ¹⁴

It was found that if network robustness can provide genetic robustness for buffering genetic variations while environmental robustness allows for resistance of environmental disturbances, then the phenotype of a biological network is robust in evolution. The trade-off between genetic robustness and environmental robustness in evolution is discussed according to the robust stochastic stabilization and noise-filtering analysis of a nonlinear stochastic gene network. The balance between network evolvability and network robustness, ie, between resisting and allowing changes in their own internal states (phenotypes) of stochastic gene networks, ¹⁵ sheds light from the systematic perspective on the mechanisms that govern the exploitation and toleration of the messiness of genetic variations and environmental disturbances in the evolutionary process. ¹⁴ In the evolutionary scenario, network robustness is an intrinsic property of evolvability and might, in the long run, improve the evolution of biological networks at all levels. This is because the accumulation of neutral mutations can result in neutral networks that provides evolutionary paths to new adaptations for the network population by random drift. ¹⁶ In this situation, network robustness and network evolvability might again be positively correlated. ¹⁷ In this study, the balance between network evolvability and network robustness in the evolutionary process will be explicitly revealed by the phenotype robustness criterion of evolutionary biological networks.

In this study, a unifying mathematical model is proposed for biological systems at different levels. According to this unifying mathematical model, the framework of phenotype robustness criteria is constructed for systems and evolutionary biology. We found that if network robustness can grant the intrinsic robustness for tolerating intrinsic parameter fluctuations, genetic robustness for buffering genetic variations, and environmental robustness for resisting environmental stimuli, then the phenotypes of these biological networks will be maintained under intrinsic parameter fluctuations, genetic variations and environmental stimuli. Moreover, the trade-off between intrinsic robustness, genetic robustness, environmental robustness, and network robustness can be revealed by the phenotype robustness criteria for biological networks at different levels. In the following sections, the aforementioned types of robustness, as well as their trade-off on phenotype robustness in systems and evolutionary biology, are discussed in sequence.

Trade-off between Intrinsic Robustness, Environmental Robustness and Network Robustness in Systems Biology

Linear gene regulatory network

Initially, for the convenience of illustration, we will consider only the following linear biochemical dynamics of a n-gene regulatory network

$$\frac{dx\left(t\right)}{dt}=Nx\left(t\right)$$

(1) where the concentration vector (network state) x(t) and stoichiometric interaction matrix N of the gene regulatory network are respectively given by

$$x\left(t\right)=\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\\ ⋮\\ ⋮\\ x_n\left(t\right)\end{array}\right],\hspace{1em}N=\left[\begin{array}{ccccc}{N}_{11}& & \cdots & & N_{1n}\\ & \ddots & & & \\ ⋮& & N_{ij}& & ⋮\\ & & & \ddots & \\ N_{n1}& & \cdots & & N_{nn}\end{array}\right]$$

in which x_i(t) denotes the concentration of the i-th gene, and N_ij denotes the regulation from gene j to gene i.

Suppose the linear gene regulatory network suffers intrinsic parameter fluctuations mainly due to random molecular fluctuations (or gene expression noises) so that stoichiometric interaction matrix N is perturbed by L random fluctuation sources as $N+\sum _{i=1}^L\Delta N_i{n}_i\left(t\right)$ , where ΔN_i denotes the effect on the gene network of the i-th random fluctuation source n_i(t), which is represented by white Gaussian noise with zero mean and unit variance to denote the stochastic part of fluctuation. That is, the stochastic part of fluctuation is absorbed by n_i(t) with dw_i(t) = n_i(t)dt, where w_i(t) is a standard Wiener process^3,4 (or Brownian motion) and the change of the gene regulatory network by the random fluctuation source n_i(t) is denoted by ΔN_i. Further, the gene regulatory network also suffers from the environmental disturbance v(t). In such situation, under intrinsic parameter fluctuations and environmental disturbances in vivo, the gene regulatory network in (1) should be modified as follows ¹⁸

$$\frac{dx\left(t\right)}{dt}=\left(N+\sum _{i=1}^L\Delta N_i{n}_i\left(t\right)\right)x\left(t\right)+v\left(t\right)$$

(2)

In the conventional notation of engineering and system science, the stochastic dynamic equation (2) for the gene regulatory network in vivo could be represented by the following Ito stochastic system^3,4

$$dx\left(t\right)=Nx\left(t\right)dt+v\left(t\right)dt+\sum _{i=1}^L\Delta N_{i}x\left(t\right)d{w}_i\left(t\right)$$

(3) where dw_i(t) = n_i(t)dt, and w_i(t) denotes the corresponding Wiener process or Brownian motion of the i-th random parametric fluctuation in the gene regulatory network. The stochastic process x(t) of the gene regulatory network in (3) can be considered as a population of gene regulatory networks over random parameter fluctuations.

Before further analysis of Ito stochastic system of gene regulatory network in (3), some definitions are given below.

Definition 1: The stochastic system in (3) with v(t) ≡ 0 is called stochastical stability, if there exists a Lyapunov (power) function V(x(t)) > 0 such that E dV (x(t))/dt ≤ 0 where E(●) denotes the expectation operation on ● over all random parameter fluctuations, ie, the average Lyapunov function of gene regulatory networks over all random parameter fluctuations (or the average power of the population of gene regulatory networks) does not increase. In this situation, the random parameter fluctuation $\sum _{i=1}^L\Delta N_i{n}_i\left(t\right)$ can be tolerated by the gene regulatory.

Definition 2: Intrinsic robustness: The ability to tolerate intrinsic parameter fluctuations without violating the stochastic stability of gene regulatory network in (3) is called intrinsic robustness.

Since w_i(t) and v(t) are stochastic processes, x(t) is also a stochastic process. Let us denote the disturbance sensitivity level ρ of the stochastic gene regulatory network in (3) as the following average energy ratio

$$\frac{E{\int }_0^{\infty }{x}^T\left(t\right)x\left(t\right)dt}{E{\int }_0^{\infty }{v}^T\left(t\right)v\left(t\right)dt}\le {\rho }^2$$

(4) where E(●) denotes the expectation of ● over population of gene regulatory networks.

The physical meaning of disturbance sensitivity level ρ in (4) is that the effect of environmental disturbance v(t) on the network state x(t) of the stochastic gene regulatory network is less than ρ from the average energy point of view. In general,^2,3,19 it is not easy to measure the system sensitivity to all possible environmental disturbances directly, and it is always to give an upper bound of the system sensitivity. Then we will decrease the upper bound to approach its minimum to obtain the system sensitivity. The disturbance sensitivity level ρ can be taken as the upper bound of the network sensitivity of a gene regulatory network to the environmental disturbance. If 0 ≤ ρ < 1, then the effect of environmental disturbance v(t) is attenuated by the gene regulatory network, meaning that it is less sensitive to environmental disturbance. If ρ > 1, then the effect of environmental disturbance is amplified by the gene regulatory network, meaning that it is more sensitive to environmental disturbance. The minimum disturbance sensitivity level ρ₀ in (4) is denoted as the network sensitivity of the gene regulatory network,^7,8 ie, ρ₀ is the lowest bound of ρ in (4).

Definition 3: the ability to resist the effect of environmental disturbance on the network state x(t) of gene regulatory network in (3) is called the environmental robustness, which is always inversely proportional to network sensitivity ρ₀.

Definition 4: The phenotype of the stochastic gene regulatory network in (3) is called robust stabilization (ie, phenotype robustness) if the gene regulatory network has enough network robustness to tolerate the intrinsic parameter fluctuations $\sum _{i=1}^L\Delta N_{i}x\left(t\right)d{w}_i\left(t\right)$ (this ability is called intrinsic robustness) and can attenuate the effect of environmental disturbance v(t) by (4) (this ability is called environmental robustness). In other words, the phenotype of the gene regulatory network can be maintained despite of intrinsic parameter fluctuations and environmental disturbances, if the network robustness can provide enough intrinsic robustness and environmental robustness.

In the linear stochastic gene regulatory network (3), let us denote the Lyapunov function of gene regulatory network in (2) as V (x(t)) = x^T (t)Px(t) for some positive definite symmetric matrix P = P^T > 0. Then the condition for robust stabilization of phenotype under intrinsic parameter fluctuations and environmental disturbances is given below.^7,8

Proposition 1 : The phenotype of the linear gene regulatory network with intrinsic parameter fluctuations and environmental disturbances is robustly stable if there exists some symmetric positive definite matrix P such that the following phenotype robustness criterion holds

$$PN+N^{T}P+\sum _{i=1}^L\Delta N_i^{T}P\Delta N_i+\frac{1}{{\rho }^2}{P}^{T}P+I\le 0$$

(5) where I is identity matrix. That is, if the above quadratic inequality holds, then the intrinsic parameter fluctuations could be tolerated, ie, the intrinsic parameter fluctuations $\sum _{i=1}^L\Delta N_{i}x\left(t\right)d{w}_i\left(t\right)$ would not destroy the network robustness, and the effect of environmental disturbances could be attenuated to a level ρ, ie, the effect of environmental disturbance v(t) on network state x(t) of the gene regulatory network is bounded below ρ in (4) from the energy point of view.

According to the Proposition 1, the network sensitivity ρ₀ of the stochastic gene regulatory network in (3) could be measured by solving the following constrained optimization problem

$$\begin{array}{l}{\rho }_0=\underset{p>0}{\min }\rho \hfill \\ \text{subject to}\left(5\right)\hfill \end{array}$$

(6)

This could be solved by decreasing ρ until no positive solution P > 0 in (5) is obtained.^7,8 Replacing ρ by ρ₀ in (5), we obtain the following phenotype robustness criterion of the gene regulatory network in (2) or (3)

$$PN+N^{T}P+\sum _{i=1}^L\Delta N_i^{T}P\Delta N_i+\frac{1}{{\rho }_0^2}{P}^{T}P+I\le 0$$

(7)

For the gene regulatory network (2) with intrinsic parameter fluctuations and environmental disturbances, if the Eigenvalues of interaction matrix N are more negative in the far left-hand side of s-domain (more robust) so that the phenotype robustness criterion (7) holds, then the gene regulatory network could tolerate more intrinsic parameter fluctuations and resist more environmental disturbances. In phenotype robustness criterion (7), -(PN + N^TP) in (7) can be taken as a measure of network robustness (a gene network with more negative Eigenvalues in N is with more network robustness, see Fig. 1); $\sum _{i=1}^L\Delta N_i^{T}P\Delta N_i$ in (7) due to the parametric fluctuation can be taken as a measure of intrinsic robustness; 1/ρ ² ₀ P^TP + I due to environmental disturbance can be taken as a measure of environmental robustness with network sensitivity ρ₀. In other words, environmental robustness is inversely proportional to network sensitivity to the environment. Therefore, the phenotype robustness criterion in (7) could be changed as

$$\underset{\mathrm{int}rinsic robustness}{\underbrace{\sum _{i=1}^L\Delta N_i^{T}P\Delta N_i}}+\underset{environmental robustness}{\underbrace{\frac{1}{{\rho }_0^2}{P}^{T}P+I}}\le \underset{network robustness}{\underbrace{-\left(PN+N^{T}P\right)}}$$

(8)

Figure 1.

The smaller distance between the locations of Eigenvalues of N and the image axis can be taken as the measure of network robustness for the linear stochastic gene network in (2).

The physical meaning of phenotype robustness criterion in (8) is that if the network robustness can confer enough intrinsic robustness for tolerating intrinsic parameter fluctuations and environmental robustness for resisting environmental disturbances, then the phenotype of the gene network is maintained. In general, feedback loop could shift the Eigenvalues of N to the more negative left-hand side in s-domain to increase network robustness in (8) for tolerating more intrinsic parameter fluctuation and for resisting more environmental disturbance (ie, lower network sensitivity with small ρ₀). Further, parallel loop and modular structure through redundancy duplication scheme could tolerate intrinsic parameter fluctuations in gene networks to reduce the influence of ΔN_i from random fluctuation source n_i(t). Feedback loop, parallel loop, and modular and redundant structure are popular in gene regulatory networks as they contribute to phenotype robustness and win natural selection's favor. ¹³

If the phenotype robustness criterion is violated, then network robustness cannot simultaneously provide enough intrinsic robustness for tolerating intrinsic parameter fluctuations or enough environmental robustness for resisting environmental disturbances simultaneously. In this situation, the phenotype of a gene network may change. In this case, we could improve network robustness by the gene circuit design method as follows^7,20

$$dx\left(t\right)=\left(N+K\right)x\left(t\right)dt+v\left(t\right)dt+\sum _{i=1}^L\Delta N_{i}x\left(t\right)d{w}_i\left(t\right)$$

(9) where K is the kinetic parameter matrix of gene circuits to be designed, in which k_ij denotes the gene circuit parameter to be specified for the gene circuit between gene j and gene i via transfection and transformation biotechnologies.^21,22 The gene circuit from gene j to gene i can be implemented by inserting the motif binding site of gene product j (ie, the protein of gene j) into the promoter region of gene i. This would enable the protein of gene j to bind to this inserted site and to act as a transcription factor (TF) regulating the gene expression of gene i. By inserting long or short motif-binding sites, we can get a large or small gene circuit parameter k_ij. The insertion of a TF-binding site into the promoter region can easily be done by using a highly efficient phage-based homologous recombination system, called recombineering.^23,24 This powerful biotechnology has been employed to engineer large segments of genomic DNA for generating transgenic and knockout constructs. However, for simplicity and feasibility, only a few practical gene circuits in Kx(t) are considered for gene circuit design in the gene network to improve network robustness; ie, the gene circuit parameter K has only a few elements k_ij. After gene circuit design as (9), the phenotype robustness criterion in (8) can be modified as follows

$$\underset{\mathrm{int}rinsic robustness}{\underbrace{\sum _{i=1}^L\Delta N_i^{T}P\Delta N_i}}+\underset{environmental robustness}{\underbrace{\frac{1}{{\rho }_0^2}{P}^{T}P+I}}\le \underset{network robustness}{\underbrace{-\left(P\left(N+K\right)+{\left(N+K\right)}^{T}P\right)}}$$

(10)

That is, the robust gene circuit designed in (9) is meant to specify K in (10) such that the improved network robustness be large enough to provide intrinsic robustness and environmental robustness.^8,25 This network robustness improvement biotechnology for genetic networks can be applied to genetic therapy. ²⁵

Remark 1: If the engineered circuit is also perturbed by $\sum _{j=1}^M\Delta K_j{n}_j\left(t\right)$ for some random noises n_j(t), then the phenotype robustness criterion in (10) should be modified as follows

$$\underset{\mathrm{int}rinsic robustness}{\underbrace{\sum _{i=1}^L\Delta N_i^{T}P\Delta N_i+\sum _{j=1}^M\Delta K_j^{T}P\Delta K_j}}+\underset{environmental robustness}{\underbrace{\frac{1}{{\rho }_0^2}{P}^{T}P+I}}\le \underset{network robustness}{\underbrace{-\left(P\left(N+K\right)+{\left(N+K\right)}^{T}P\right)}}$$

(11)

Nonlinear gene regulatory network

In real biological systems, the gene regulatory networks are always nonlinear. In this situation, the n-gene stochastic gene regulatory network in (3) should be modified as

$$dx\left(t\right)=N\left(x\left(t\right)\right)dt+v\left(t\right)dt+\sum _{i=1}^L\Delta N_i\left(x\left(t\right)\right)d{w}_i\left(t\right)$$

(12) where N(x(t)) denotes the nonlinear regulation vector and ΔN_i(x(t)) denotes the nonlinear intrinsic parameter fluctuation due to the i-th random fluctuation source w_i(t).

Consider the nonlinear stochastic system in (12). There exist many equilibrium points (phenotypes). Suppose a phenotype near a stable equilibrium point x_e is of interest (see Fig. 2). For convenience of analysis, the origin of the nonlinear stochastic gene regulatory network is shifted to the equilibrium point (phenotype) x_e to simplify the robust analysis procedure. Let us denote $\tilde{x}$ (t) = x(t) - x_e, then the following shifted nonlinear stochastic gene regulatory network can be obtained as follows^4,5

$$d\tilde{x}\left(t\right)=N\left(\tilde{x}\left(t\right)+x_e\right)dt+v\left(t\right)dt+\sum _{i=1}^L\Delta N_i\left(\tilde{x}\left(t\right)+x_e\right)d{w}_i\left(t\right)$$

(13) where N( $\tilde{x}$ (t)) and ΔN( $\tilde{x}$ (t)) denote N( $\tilde{x}$ (t) + x_e) and ΔN( $\tilde{x}$ (t) + x_e), respectively, for the convenience of notation because x_e is a constant. That is, the origin $\tilde{x}$ (t) = 0 of the nonlinear stochastic gene network in (13) is at the equilibrium point x_e of the original nonlinear stochastic gene regulatory network in (12). Then, let us consider the robust stabilization of phenotype at x_e of the nonlinear stochastic gene regulatory network in (12) or (13). According to the stochastic Lyapunov stability theory,^4,5 we get the following phenotype robustness for the nonlinear stochastic gene network. ⁷

Figure 2.

The stochastic nonlinear gene network has many local stable equilibrium points (phenotypes).

Proposition 2: The phenotype x_e of the stochastic nonlinear gene network under intrinsic parameter fluctuations and environmental disturbances in (13) is robustly stable if the following Hamilton Jacobi Inequality (HJI) holds for some positive function V( $\tilde{x}$ (t)) > 0

$$\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}N\left(\tilde{x}\left(t\right)\right)+\frac{1}{2}\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right)\frac{{\partial }^{2}V\left(\tilde{x}\left(t\right)\right)}{\partial {\tilde{x}}^2\left(t\right)}\Delta N_i\left(\tilde{x}\left(t\right)\right)+\frac{1}{4{\rho }^2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+{\tilde{x}}^T\left(t\right)\tilde{x}\left(t\right)\le 0$$

(14)

That is, if the HJI in (14) holds for some V( $\tilde{x}$ (t)) > 0, then the intrinsic parameter fluctuations $\sum _{i=1}^L\Delta N_i\left(\tilde{x}\left(t\right)\right)d{w}_i\left(t\right)$ could be tolerated by the nonlinear stochastic gene network and the effect of environmental disturbance v(t) could be attenuated to a level ρ in (4).

The network sensitivity ρ₀ of the nonlinear stochastic gene network in (13) can then be obtained by solving the following constrained optimization problem ⁸

$$\begin{array}{l}{\rho }_0=\underset{P>0}{\min }\rho \hfill \\ \text{subject to HJI in}\left(14\right)\hfill \end{array}$$

(15)

Substituting ρ₀ into ρ in (14), we get the following phenotype robustness criterion from the nonlinear stochastic gene regulatory network in (12)

$$\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}N\left(\tilde{x}\left(t\right)\right)+\frac{1}{2}\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right)\frac{{\partial }^{2}V\left(\tilde{x}\left(t\right)\right)}{\partial {\tilde{x}}^2\left(t\right)}\Delta N_i\left(\tilde{x}\left(t\right)\right)+\frac{1}{4{\rho }_0^2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+{\tilde{x}}^T\left(t\right)\tilde{x}\left(t\right)\le 0$$

(16)

From the phenotype robustness criterion in (16), it can be seen that if $\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)N\left(\tilde{x}\left(t\right)\right)$ is more negative, the nonlinear gene regulatory network will be more robustly stable and more able to tolerate the term 1/2 $\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right){\partial }^{2}V\left(\tilde{x}\left(t\right)\right)/\partial {\tilde{x}}^2\left(t\right)\Delta N_i\left(\tilde{x}\left(t\right)\right)$ in (16) due to the effect of intrinsic parameter fluctuations and the term 1/4 ${\rho }_0^2{\left(\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)\right)}^T\left(\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)\right)+{\tilde{x}}^T\left(t\right)\tilde{x}\left(t\right)$ and due to the effect of environmental disturbances. The phenotype robustness criterion in (16) can also be rewritten as

$$\underset{\mathrm{int}rinsic robustness}{\underbrace{\frac{1}{2}\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right)\frac{{\partial }^{2}V\left(\tilde{x}\left(t\right)\right)}{\partial {\tilde{x}}^2\left(t\right)}\Delta N_i\left(\tilde{x}\left(t\right)\right)}}+\underset{environmental robustness}{\underbrace{\frac{1}{4{\rho }_0^2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+{\tilde{x}}^T\left(t\right)\tilde{x}\left(t\right)}}\le \underset{network robustness}{\underbrace{-\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}N\left(\tilde{x}\left(t\right)\right)\right)}}$$

(17)

The phenotype robustness criterion in (17) for the nonlinear stochastic gene regulatory network in (12) can be denoted as “intrinsic robustness + environmental robustness ≤ network robustness.” In other words, if the sum of intrinsic robustness and environmental robustness is still less than the network robustness, then the phenotype of the nonlinear gene network remains robust under the influence of intrinsic parameter fluctuations and environmental disturbances.

In order to maintain phenotype robustness in (17), the network structure needs to make the term $\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)N\left(\tilde{x}\left(t\right)\right)$ more negative (with a deeper basin at the corresponding equilibrium point x_e in Fig. 2) through feedback loop. With this, the network robustness will be large enough to tolerate more intrinsic parameter fluctuations and the network sensitivity ρ₀ small enough to resist environmental disturbances. In addition, modular and redundant structures, which could reduce the effect of variations ΔN_i( $\tilde{x}$ (t)), are both favored by natural selection in the evolutionary process of the nonlinear gene regulatory network. Furthermore, the trade-off between intrinsic robustness and environmental robustness can be seen from the phenotype robustness criterion in (17). The nonlinear gene regulatory network cannot confer simultaneously large intrinsic robustness and large environmental robustness. The total amount of intrinsic robustness plus environmental robustness cannot exceed network robustness, and therefore there should be a trade-off between them.

If the gene regulatory network in (12) does not have sufficient network robustness such that the phenotype robustness in (17) cannot be guaranteed, then robust stabilization cannot be maintained at the equilibrium point x_e and the phenotype of the gene regulatory network may be shifted by intrinsic parameter fluctuations and environmental disturbances from x_e to another equilibrium point (another phenotype, see Fig. 2). In this situation, in order to maintain the phenotype at x_e, some gene circuits could be designed to improve the network robustness as follows

$$d\tilde{x}\left(t\right)=\left(N\left(\tilde{x}\left(t\right)\right)+K\left(\tilde{x}\left(t\right)\right)\right)dt+v\left(t\right)dt+\sum _{i=1}^L\Delta N_i\left(\tilde{x}\left(t\right)\right)d{w}_i\left(t\right)$$

(18) where K( $\tilde{x}$ (t)) is the designed gene circuit for improving network robustness of the nonlinear stochastic gene network at the equilibrium point x_e by transfection and transformation biotechnologies.^21,22

After the gene circuit is designed as (18), the phenotype robustness criterion in (17) should be modified as

$$\underset{\mathrm{int}rinsic robustness}{\underbrace{\frac{1}{2}\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right)\frac{{\partial }^{2}V\left(\tilde{x}\left(t\right)\right)}{\partial {\tilde{x}}^2\left(t\right)}\Delta N_i\left(\tilde{x}\left(t\right)\right)}}+\underset{environmental robustness}{\underbrace{\frac{1}{4{\rho }_0^2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+{\tilde{x}}^T\left(t\right)\tilde{x}\left(t\right)}}\le \underset{network robustness}{\underbrace{-\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\left(N\left(\tilde{x}\left(t\right)\right)+K\left(\tilde{x}\left(t\right)\right)\right)\right)}}$$

(19)

That is, the robust gene network design in (18) is to specify the designed gene circuit K( $\tilde{x}$ (t)) for making network robustness in the right-hand side of (19) as large a value as possible (ie, with a very deep basin at the equilibrium point x_e in Fig. 2) to override both intrinsic robustness and environmental robustness in the left-hand side of (19). In this situation, the phenotype of the nonlinear stochastic gene regulatory network is maintained in the basin near the equilibrium point x_e despite of intrinsic parameter fluctuations and environmental disturbances.

In general, it is very difficult to solve the HJI in (14) or (17) for the phenotype robustness criterion of the nonlinear stochastic gene network in (12) so as to gain more insight into the trade-off between intrinsic robustness, environmental robustness, and network robustness from the systematic point of view, such as the linear stochastic gene network. At present there is no good method for solving the nonlinear partial differential HJI either analytically or numerically. In this situation, some interpolation methods like global linearization ²⁶ and fuzzy techniques ² are employed to interpolate several local linearized gene networks in order to approximate the nonlinear stochastic gene regulatory network and to simplify the systematic analysis of phenotype robustness criterion in (14) or (17). When using the global linearization technique, all the global linearizations of the nonlinear stochastic gene regulatory network in (12) are bounded by a polytope consisting of M vertices as ²⁶

$$\left[\begin{array}{c}\frac{\partial N\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\\ \frac{\partial \Delta N_j\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\\ \frac{\partial K\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\end{array}\right]\in Co\left[\left[\begin{array}{c}{N}_1\\ \Delta N_{j1}\\ K_1\end{array}\right]\cdots \left[\begin{array}{c}{N}_i\\ \Delta N_{ji}\\ K_i\end{array}\right]\cdots \left[\begin{array}{c}{N}_m\\ \Delta N_{jm}\\ K_m\end{array}\right]\right],\begin{array}{c}i=1,2\cdots M\\ j=1,2\dots L\end{array}$$

(20) where Co denotes the convex hull of polytope with M vertices; ie, if all the linearized systems are inside the convex hull Co, then the state trajectories $\tilde{x}$ (t) of the nonlinear shifted gene network in (13) will belong to the convex combination of the state trajectories of the following M linearized gene regulatory networks, derived from the vertices of polytope ²⁶

$$d\tilde{x}\left(t\right)=N_i\tilde{x}\left(t\right)dt+v\left(t\right)dt+\sum _{j=1}^L\Delta N_{ji}\tilde{x}\left(t\right)d{w}_j\left(t\right),i=1,2\dots M$$

(21)

Remark 2: If we considered only the linearization at the origin $\tilde{x}$ (t) ≡ 0, then the global linearization system of the gene regulatory network in (21) is reduced as the linear gene regulatory network in (2).

According to the global linearization theory, ²⁶ if (20) holds, then every trajectory of the nonlinear stochastic gene network in (12) or (13) can be represented by a convex combination of M linearized stochastic gene networks in (21). Therefore, if we can prove that the convex interpolation of M linearized stochastic gene networks in (21) has phenotype robustness under intrinsic parameter fluctuations and environmental disturbances, then the original nonlinear stochastic gene regulatory network in (12) or (13) will have the same phenotype robustness. The convex combination of M linearized stochastic gene networks in (21) could be represented by

$$d\tilde{x}\left(t\right)=\sum _{i=1}^M{\alpha }_i\left(\tilde{x}\left(t\right)\right)\left[N_i\tilde{x}\left(t\right)dt+v\left(t\right)dt+\sum _{j=1}^L\Delta N_{ji}\tilde{x}\left(t\right)d{w}_j\left(t\right)\right]$$

(22) where the interpolation function α_i( $\tilde{x}$ (t)) satisfies 0 ≤ α_i( $\tilde{x}$ (t)) ≤ 1 and $\sum _{i=1}^M{\alpha }_i\left(\tilde{x}\left(t\right)\right)=1$ . In other words, the trajectory of the nonlinear stochastic gene network in (12) or (13) can be represented by the trajectory of the interpolated gene network in (22), which is the convex interpolation of M linearized gene regulatory networks at M vertices in (20). There are several kinds of interpolation techniques for interpolating some local linear systems to approximate a nonlinear system including T-S fuzzy interpolation method through fuzzy bases.^7,27,28

After the nonlinear stochastic gene network in (12) or (13) is represented by the interpolation of M linearized gene networks in (22), we then get the following result of the phenotype robustness criterion for the nonlinear stochastic gene regulatory network.^7,18

Proposition 3: If the following quadratic inequality holds with P > 0 for the phenotype robustness of the globally linearized gene network in (22)

$$P{N}_i+N_i^{T}P+\sum _{j=1}^L\Delta N_{ji}^{T}P\Delta N_{ji}+\frac{1}{{\rho }^2}{P}^{T}P+I\le 0,i=1,2\dots M$$

(23) then the nonlinear stochastic gene network in (22) or (13) is robustly stable at the equilibrium point x_e under intrinsic parameter fluctuations and environmental disturbances. In other words, if the phenotype robustness criterion in (23) holds, then the intrinsic parameter fluctuations can be tolerated and the effect of environmental disturbances can be attenuated to a level ρ.

Similarly, the network sensitivity of the global linearized gene network in (22) can be obtained by solving the following constrained optimization problem

$$\begin{array}{l}{\rho }_0=\underset{P>0}{\min }\rho \hfill \\ \text{subject to}\left(23\right)\hfill \end{array}$$

(24) which could be easily solved with the help of linear matrix inequality (LMI) toolbox in Matlab. Following the global linearization technique in (20)–(22), the phenotype robustness criterion in (17) for nonlinear gene regulatory network can be modified as

$$\underset{local Intrinsic robustness}{\underbrace{\sum _{j=1}^L\Delta N_{ji}^{T}P\Delta N_{ji}}}+\underset{local environmental robustness}{\underbrace{\frac{1}{{\rho }_0^2}{P}^{T}P+I}}\le \underset{local network robustness}{\underbrace{-\left(P{N}_i+N_i^{T}P\right)}},i=1,2\dots M$$

(25)

From the phenotype robust criterion in (25), it is seen that if the local network robustness is large enough for every local linearized gene network to simultaneously provide enough local intrinsic robustness for tolerating local random parameter fluctuations and enough local environmental robustness for resisting local environmental disturbances, then the equilibrium point x_e (phenotype) is robustly stable in spite of intrinsic parameter fluctuations and environmental disturbances. The trade-offs between the local intrinsic robustness, local environmental robustness, and local network robustness can also be seen in the phenotype robustness criterion (25) of local linearized gene regulatory networks of the nonlinear stochastic gene network.

In the case of gene network design for engineering some gene circuits to improve the network robustness in (18), it can also be approximated by the global linearization representation as follows

$$d\tilde{x}\left(t\right)=\sum _{i=1}^M{\alpha }_i\left(\tilde{x}\left(t\right)\right)\times \left[\left(N_i+K_i\right)\tilde{x}\left(t\right)dt+v\left(t\right)dt+\sum _{j=1}^L\Delta N_{ji}\tilde{x}\left(t\right)d{w}_j\left(t\right)\right]$$

(26)

In this situation, the phenotype robustness criterion in (25) should be modified as follows

$$\underset{local Intrinsic robustness}{\underbrace{\sum _{j=1}^L\Delta N_{ji}^{T}P\Delta N_{ji}}}+\underset{local environmental robustness}{\underbrace{\frac{1}{{\rho }_0^2}{P}^{T}P+I}}\le \underset{local network robustness}{\underbrace{-\left[\left(P\left(N_i+K_i\right)+{\left(N_i+K_i\right)}^{T}P\right)\right]}},i=1,2\dots M$$

(27)

The gene circuit design is then able to specify adequate K_i, i = 1, … M in order for the local network robustness in the right-hand side of (27) to be large enough to provide sufficient local intrinsic robustness for tolerating local intrinsic parameter fluctuations as well as enough local environmental robustness for resisting environmental disturbances.

Trade-off between Genetic Robustness, Environmental Robustness and Network Robustness in Evolutionary Biology

Linear stochastic evolutionary gene regulatory network

In the evolutionary process, a gene regulatory network will suffer from genetic variations and environmental changes or stresses such as temperature or salinity, which may perturb its phenotype. Consider the following stochastic gene regulatory network in the evolutionary process

$$\begin{array}{c}dx\left(t\right)=\left(N+\sum _{i=1}^L\Delta N_{i}P\left(t-t_i\right)\right)x\left(t\right)dt+v\left(t\right)dt\\ =Nx\left(t\right)dt+v\left(t\right)dt+\sum _{i=1}^L\Delta N_{i}x\left(t\right)P\left(t-t_i\right)dt\end{array}$$

(28) where $\sum _{i=1}^L\Delta N_{i}P\left(t-t_i\right)=\sum _{i=1}^L\Delta N_i$ is a weighted Poisson point process with intensity to denote the effect of discontinuous genetic mutations on the phenotype variation, ²⁹ in which $L=\sum _{i=1}^{L}P\left(t-t_i\right)=\sum _{i=1}^{L}1$ denotes a Poisson point (counting) process with mean λt and variance λt. ΔN_i denotes the random phenotypic variation due to the point genetic variation P(t - t_i) occurring at t = t_i. v(t) denotes the environmental stimuli. The stochastic gene regulatory network in (28) can be considered as a population of gene regulatory networks over all possible random genetic variations.

Remark 3: Unlike the parameter variations in (3), which are mainly due to random molecular fluctuations in the conventional gene regulatory network, that are always modeled as continuous Brownian motion (Wiener process), the parameter variations of the gene regulatory network in the evolutionary process are mainly due to discontinuous genetic mutations. These mutations are hereditable and will be accumulated to affect the phenotype of offspring in the evolutionary process. Therefore, phenotype variations due to genetic mutations are more likely random point processes described by the weighted Poisson processes. ²⁹ Additionally, the time scale of stochastic Poisson point process of evolutionary biology in (28) is much longer than the stochastic Wiener process of systems biology in (3).

Remark 4: From the population point of view, the stochastic gene evolutionary network in (28) with weighted Poisson point process can be considered as a population of gene networks over all possible random genetic variations $\sum _{i=1}^L\Delta N_{i}x\left(t\right)P\left(t-t_i\right)$ in the evolutionary process. The network population in (28) will be selected by natural selection to meet fitness function in the evolutionary process. Therefore, if the robust stability results hold from the stochastic network point of view, they could also hold from network population point of view.

Let us denote the phenotype variation of the evolutionary gene regulatory network as $\tilde{x}$ (t) = x(t) - x_e. The effectiveness of environmental disturbances or stimuli on phenotype variation in the evolutionary process (ie, the response of the gene regulatory network to environmental changes in the evolutionary process) is denoted as the following evolutionary level

$$\frac{E{\int }_0^{t_p}{\tilde{x}}^T\left(t\right)\tilde{x}\left(t\right)dt}{E{\int }_0^{t_p}{v}^T\left(t\right)v\left(t\right)dt}\le {\rho }^2$$

(29) where E(●) denotes the expectation operation on • over all random genetic variations (or population of genetic networks) and is similar to the disturbance sensitivity level in (4). That is, the evolution level ρ denotes the upper boundary of network evolvability ρ₀ of the stochastic gene regulatory network in (28) from the average energy viewpoint of network population. If evolution level ρ in (29) is large, then environmental stimuli will exert more influence on the phenotype variation, and vice versa. A small evolution level ρ denotes that the gene regulatory network has a reasonable ability to resist the effect of environmental stimuli in the evolutionary process. The numerator in (29) denotes the average energy of all possible phenotype variations in network population till the present time t_p. The evolution level ρ denotes the upper boundary of the normalized influence of environmental stimuli on the phenotype variations in the evolutionary process from the total energy point of view. The fitness function f for the phenotype x_e of the stochastic linear gene regulatory network (28) in the evolutionary process is inversely proportional to the evolution level ρ on the gene regulatory network; ie,

$$f=\frac{1}{\rho }$$

(30)

From (30), it is seen that a gene regulatory network with less evolution level to environmental stimuli exhibits more fitness function to the phenotype $\tilde{x}$ _e = 0, and vice versa. Therefore, the maximization of fitness function for a phenotype ˜_e = 0 in evolution is equivalent to the minimization of the evolution level ρ in (29), ie,

$$\max f=\frac{1}{\min \rho }$$

(31)

Let us denote the network fitness f₀ = max f and the network evolvability ρ₀ = min ρ of the gene regulatory network. Then we get

$$f_0=\frac{1}{{\rho }_0}$$

(32) ie, the network fitness f₀ of gene regulatory network at the phenotype x_e (ie, at the equilibrium point $\tilde{x}$ _e = 0) occurs at the case of minimum evolution level (network evolvability) ρ₀.

Remark 5: The fitness function in (30) is chosen according to the robust stabilization of phenotype $\tilde{x}$ _e = 0, ie, the capacity for resistance of the stochastic gene regulatory network against environmental stimuli v(t) in the evolutionary process.

However, it is still very difficult to solve the problem of network fitness f₀ or network evolvability ρ₀ of the stochastic gene regulatory network in (31) or (32) directly. A suboptimal method is employed to solve by minimizing its upper bound ρ in (29). In the following, the upper bound ρ in (29) is decreased, until its minimum value approaches the network evolvability ρ₀ in (31) or (32), ie, solving ρ₀ by the so-called suboptimal method. Then we get the following result for evolution level ρ of the stochastic gene regulatory network in (28).

Proposition 4: For the stochastic gene regulatory network in (28) with genetic variations and environmental stimuli in the evolutionary process, if the following phenotype robustness criterion holds for P > 0.

$$N^{T}P+PN+PP+\sum _{i=1}^L\Delta N_i^T\Delta N_i+\sum _{i=1}^L\lambda \Delta N_i^{T}P\Delta N_i+\frac{1}{{\rho }^2}PP+I\le 0$$

(33)

where λ is the intensity (mean) of genetic variations described by Poisson point process P(t - t_i) in (28), then the stochastic gene regulatory network has enough network robustness to tolerate random genetic variations. It also has an evolution level ρ in (29) capable of responding to the environmental stimuli. In other words, the stochastic gene regulatory network in (28) is robustly stable with an evolution level ρ in the evolutionary process.

Proof: see appendix A

Remark 6: The genetic variations of the evolutionary gene regulatory network used by Chen and Lin ¹⁴ are considered with the Wiener process. In this study, the random genetic variations of the gene regulatory network are considered with Poisson point process to mimic discontinuous genetic mutations, which are heritable in the evolutionary process.

Since the evolution level ρ is the upper bound of the evolvability of stochastic gene regulatory network, the network evolvability ρ₀ of the gene regulatory network is obtained as follows

$$\begin{array}{l}{\rho }_0=\underset{P>0}{\min }\rho \hfill \\ \text{subject to}\left(33\right)\hfill \end{array}$$

(34)

The network evolvability ρ₀ in (34) can be taken as the smallest upper bound of network evolution to resist environmental stimuli without violation of network stability. The constrained optimization problem for solving the network evolvability ρ₀ for the stochastic gene regulatory network in (28) can be achieved by decreasing ρ until no solution P > 0 exists in (33). By using the Schur complement, ²⁶ the quadratic inequality in (33) can be transformed into the following equivalent LMI:

$$\left[\begin{array}{cccc}{N}^{T}P+& PN& +I+\sum _{i=1}^L\Delta N_i^T\Delta N_i+\sum _{i=1}^L\lambda \Delta N_i^{T}P\Delta N_i& \hspace{1em}P\\ & P& & -\frac{{\rho }^2}{1+{\rho }^2}I\end{array}\right]\le 0$$

(35)

Therefore, the network evolvability problem in (34) is equivalent to solving the following constrained optimization problem

$$\begin{array}{l}{\rho }_0=\underset{P<0}{\min }\rho \hfill \\ \text{subject to LMI in}\left(35\right)\hfill \end{array}$$

(36)

From the analysis above, the network evolvability ρ₀ in (36) for the stochastic gene regulatory network in (28) could be measured by decreasing ρ until the LMI in (35) has no solution P > 0. The solution P > 0 of LMI in (35) could be easily obtained with the help of the LMI toolbox in Matlab. ²⁶

After obtaining the network evolvability ρ₀ by solving the constrained optimization problem in (36), according to the fact that the network fitness f₀ = 1/ρ₀ (ie, a small network evolvability implies a large network fitness and vice versa), the phenotype robustness criterion in (33) can then be modified as

$$\underset{genetic robustness}{\underbrace{PP+\sum _{i=1}^L\Delta N_i^T\Delta N_i+\sum _{i=1}^L\lambda \Delta N_i^{T}P\Delta N_i}}+\underset{environmental robustness}{\underbrace{f_0^{2}PP+I}}\le -\underset{network robustness}{\underbrace{\left(N^{T}P+PN\right)}}$$

(37) in which the first three terms are due to genetic variations and the second two terms are due to environmental disturbance (see the proof in Appendix A). From (37), it can be seen that the phenotype of the genetic network can be maintained in the evolutionary process if the network robustness can confer genetic robustness for tolerating genetic variations and environmental robustness for resisting environmental stimuli in the evolutionary process. The phenotype robustness criterion for evolutionary biology in (37) is quite similar to the phenotype robustness criterion for conventional systems biology except that the conventional Wiener random molecular fluctuations are replaced by the Poisson genetic variations in evolution. It is seen that the effect of Poisson genetic variations is much larger than that of Weiner random molecular fluctuations because of the accumulated heredity of Poisson genetic variations. Therefore, network robustness must demonstrate a higher level of genetic robustness than intrinsic robustness in (8) for systems biology to maintain phenotype robustness in (37) for evolutionary biology. The trade-off between evolvability and network robustness of gene regulatory networks in the evolutionary process is discussed in the following.

For the stochastic gene network in (28) with genetic variations and environmental stimuli, if the network robustness in the right-hand side of phenotype robustness criterion in (37) is fixed, the sum of genetic robustness and the environmental robustness should be also fixed. That is, large genetic robustness will lead to small environmental robustness, and vice versa. This relationship implies that a genetic gene network cannot tolerate a large amount of genetic variations and resist a large amount of environmental stimuli simultaneously. There should be a trade-off between genetic robustness and environmental robustness in the evolutionary process. In Lenski et al's study, ³⁰ the evolutionary cases of genetic robustness were discussed using different evolutionary scenarios. The correlation with environmental robustness is considered to be the most probable cause of genetic robustness in evolution. According to the congruence scenario, the genetic robustness of the gene regulatory network in evolution is a by-product of environmental robustness of the gene regulatory network to resist environmental disturbances in the evolutionary process. This is caused by the fact that environmental disturbances are more frequent than genetic perturbations in evolution, which have been confirmed in RNA folding and heat-shock protein. ³⁰ However, this correlation between genetic robustness and environmental robustness is obvious in the phenotype robustness criterion in (37). In order to provide a buffer against environmental disturbances (ie, the term f ² _o PP + I in (37)), the gene regulatory network needs feedback loops as well as modular and redundant mechanisms to improve its network robustness, resulting in a large - (N^T P + PN) in (37), which can also tolerate genetic variations (ie, the term $PP+\sum _{i=1}^L\Delta N_i^T\Delta N_i+\sum _{i=1}^L\lambda \Delta N_i^{T}P\Delta N$ simultaneously. In order to resist the effect of a large amount of environmental disturbances in evolution, the gene regulatory network requires greater network robustness to maintain its phenotype. However, this greater network robustness can also provide a buffer against large amount of neutral genetic variations simultaneously in the evolutionary process. This large amount of neutral genetic variations may provide raw materials for new evolutionary possibilities. Therefore, the phenotype robustness criterion in (37) could give more insight into the trade-off between evolvability and network robustness.

Nonlinear stochastic evolutionary gene regulatory network

In the real nonlinear stochastic gene regulatory network in evolution, the stochastic dynamic model in (28) should be modified as

$$dx\left(t\right)=\left(N\left(x\left(t\right)\right)+v\left(t\right)\right)dt+\sum _{i=1}^L\Delta N_i\left(x\left(t\right)\right)P\left(t-t_i\right)dt$$

(38) where N(x(t)) denotes the nonlinear regulatory function and ΔN_i (x(t)) denotes the nonlinear stochastic parameter variation due to the i-th Poisson genetic variation P(t - t_i) at t = t_i in the evolutionary process. Sources of these random genetic variations include DNA mutation, deletion, duplication, inversion, and translation of chromosomes in the evolutionary process. The nonlinear stochastic gene regulatory network in (38) represents the phenotype heterogeneity in a gene network population, due to random genetic variations and environmental disturbances in the evolutionary process.

Let us consider the nonlinear stochastic gene regulatory network in evolution. Many equilibrium points (phenotypes) exist (see Fig. 2). Suppose the phenotype of the nonlinear stochastic gene regulatory network is near the equilibrium point x_e, ie, the phenotype at x_e is favored by natural selection. For convenience of systematic analysis, the origin of the nonlinear stochastic gene network in (38) is shifted to the equilibrium point (phenotype) x_e. In this situation, if the shifted nonlinear stochastic gene regulatory network is robustly stable at the origin in spite of genetic variations and environmental disturbances in the evolutionary process, then the equilibrium point x_e of phenotype is also robustly stable. Thus, the systematic analysis procedure of the network evolvability and network robustness of nonlinear stochastic evolutionary gene network in (38) is simplified. Let us denote $\tilde{x}$ (t) = x(t) - x_e, then we get the following shifted nonlinear stochastic evolutionary gene network

$$d\tilde{x}\left(t\right)=\left(N\left(\tilde{x}\left(t\right)\right)+v\left(t\right)\right)dt+\sum _{i=1}^L\Delta N_i\left(\tilde{x}\left(t\right)\right)P\left(t-t_i\right)dt$$

(39) ie, the origin $\tilde{x}$ (t) ≡ 0 of the nonlinear stochastic evolutionary gene network in (39) is at the equilibrium point x_e of the nonlinear stochastic evolutionary gene network in (38). Let us consider the robust phenotype stabilization problem of the phenotype at $\tilde{x}$ (t) ≡ 0 for the gene regulatory network in (39) in the evolutionary process. According to the nonlinear stochastic gene regulatory network in (39) and the evolution level ρ defined in (29), the following robust stabilization result is obtained for the evolutionary gene network.

Proposition 5: For a nonlinear stochastic evolutionary gene network in (39), if the following phenotype robustness criterion holds for some positive function V( $\tilde{x}$ ) > 0

$${\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^{T}N\left(\tilde{x}\left(t\right)\right)+\frac{1}{2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+\frac{1}{2}\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right)\Delta N_i\left(\tilde{x}\left(t\right)\right)+\sum _{i=1}^L\left(\begin{array}{l}V\left(\tilde{x}\left(t\right)+\Delta N_i\left(\tilde{x}\left(t\right)\right)\right)-V\left(\tilde{x}\left(t\right)\right)\\ -{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\Delta N_i\left(\tilde{x}\left(t\right)\right)\end{array}\right)\lambda +\frac{1}{4{\rho }^2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+\tilde{x}{\left(t\right)}^T\tilde{x}\left(t\right)\le 0$$

(40) then the stochastic gene regulatory network in (38) or (39) is robustly stable at the equilibrium point x_e (or $\tilde{x}$ _e = 0) with an evolution level ρ in the evolutionary process.

Proof: see appendix B

Remark 7: (i) If the nonlinear stochastic evolutionary gene network in (38) or (39) is free of environmental disturbance v(t), ie, v(t) ≡ 0 for all time and the phenotype robustness criterion in (40) holds, then the network state x(t) will stochastically achieve the origin, ie, $\tilde{x}$ (t) → 0 (or x(t) → x_e). (ii) If the phenotype criterion in (40) holds for the nonlinear stochastic evolutionary gene regulatory network in (39) with v(t) ≠ 0, then the network state x(t) will be maintained in the basin near the equilibrium point x_e. In this situation, the network state of gene regulatory network cannot asymptotically converge to x_e but will fluctuate in the basin of x_e with an evolution level ρ to respond to environmental disturbances. (iii) As both the accumulation of genetic variations (ie, the second to the fourth terms in (40)) and the evolutionary response to environmental disturbances (ie, the last two terms in (40)) are so large that the network robustness is not sufficiently significant (the first term is not negative enough or the cliff of the equilibrium point x_e in Fig. 2 is not steep enough) to override the effects of genetic variations and environmental disturbances on the phenotype simultaneously. Thus, the network state x(t) of the evolutionary gene regulatory network cannot be maintained in the basin of one stable equilibrium point x_e and will be pushed to another equilibrium point (another phenotype).

In the case of phenotype robustness criterion in (40) for the nonlinear stochastic gene regulatory network in (38) or (39) in evolution, the network evolvability can be measured as follows

$$\begin{array}{l}{\rho }_0=\underset{P>0}{\min }\rho \hfill \\ \text{subject to HJI in}\left(40\right)\hfill \end{array}$$

(41)

After solving the constrained optimization problem in (41) for the network evolvability of the nonlinear stochastic gene regulatory network in (38) or (39), the phenotype robustness criterion in (40) be modified with network fitness f₀=1/ρ₀ in (38) as follows

$$\underset{genetic robustness}{\underbrace{\left(\begin{array}{l}\frac{1}{2}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+\frac{1}{2}\sum _{i=1}^L\Delta N_i^T\left(\tilde{x}\left(t\right)\right)\Delta N_i\left(\tilde{x}\left(t\right)\right)\\ +\sum _{i=1}^L\left(V\left(\tilde{x}\left(t\right)+\Delta N_i\left(\tilde{x}\left(t\right)\right)\right)-V\left(\tilde{x}\left(t\right)\right)-{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\Delta N_i\left(\tilde{x}\left(t\right)\right)\right)\lambda \end{array}\right)}}+\underset{environmental robustness}{\underbrace{\frac{f_0^2}{4}{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^T\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)+\tilde{x}{\left(t\right)}^T\tilde{x}\left(t\right)}}\le \underset{network robustness}{\underbrace{{\left(\frac{\partial V\left(\tilde{x}\left(t\right)\right)}{\partial \tilde{x}\left(t\right)}\right)}^{T}f\left(\tilde{x}\left(t\right)\right)}}$$

(42)

In (42), the network robustness should be large enough when the nonlinear stochastic evolutionary gene network in (38) or (39) suffers from simultaneous genetic variations and environmental stimuli in evolution. That is, ${\left(\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)\right)}^{T}f\left(\tilde{x}\left(t\right)\right)$ on the right-hand side of (42) must be more negative to enable network robustness. This grants enough genetic robustness in the first three terms for tolerating genetic variations while simultaneously providing enough environmental robustness in the last two terms to resisting environmental stimuli so that the phenotype of the gene regulatory network is resilient in the evolutionary process.

The phenotype robustness criterion in (42) for the nonlinear stochastic evolutionary gene regulatory network in (38) or (39) can be denoted as “genetic robustness + environmental robustness ≤ network robustness.” In other words, if the sum of genetic robustness and environmental robustness is still less than the network robustness in evolution, then the phenotype of the evolutionary gene regulatory network can be maintained under genetic variations and environmental disturbances in the evolutionary process. The phenotype robustness criterion of the nonlinear stochastic evolutionary gene network in (38) or (39) shows that if greater network robustness is evolved to allow the gene network to provide a buffer against more environmental disturbances in evolution, it can also provide a buffer against heritable genetic variations. Obviously, the correlation between genetic robustness and environmental robustness affects the network robustness of gene network in the evolutionary process. ¹⁴

As the phenotype robustness criterion in (42) holds, the accumulated genetic variations $\sum _{i=1}^L\Delta N_i\left(x\left(t\right)\right)P\left(t-t_i\right)$ in (38) and the effect due to environmental disturbances v(t) cannot impel the nonlinear stochastic gene network to move from one stable basin to another in order to change the phenotype. In this case, the evolution is still in stasis state under the influence of genetic variations and environmental disturbances. However, if these Poisson genetic variations and stress-induced genetic variations persist and accumulate for a long time, the phenotype robustness criterion in (42) will be violated. The phenotype transition of this nonlinear stochastic evolutionary gene network from the basin of one equilibrium point to another basin of other equilibrium point may then occur to change the phenotype.

In general, the cost of evolutionary strategy to resist environmental disturbances and to tolerate genetic variations is much higher. ³¹ In this case, if stress-avoidance strategy is imposed in the evolutionary process, ³¹ the basin of equilibrium point x_e in Fig. 2 becomes shallow; ie, the network robustness term ${\left(\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)\right)}^{T}f\left(\tilde{x}\left(t\right)\right)$ in (42) is less negative. In this case, the network evolvability ρ₀ must be large and the environmental robustness becomes small in evolution, so that the nonlinear stochastic gene network is sufficiently flexible to respond to environmental stimuli with rapid evolutionary changes. If the network robustness is very large, ie, ${\left(\partial V\left(\tilde{x}\left(t\right)\right)/\partial \tilde{x}\left(t\right)\right)}^{T}f\left(\tilde{x}\left(t\right)\right)$ is very negative, the basin of the equilibrium point x_e in Figure 2 becomes deep. In this situation, large network robustness allows genetic robustness to tolerate a large amount of neutral genetic mutations in the evolutionary process. At the first sight, one may expect large network robustness to slow or even stop the evolution of a phenotype in the nonlinear stochastic evolutionary gene network. However, because a robust phenotype (deep basin) can harbor a large amount of neutral genetic variations, the robust phenotypes of the nonlinear stochastic evolutionary gene networks might show increased rather than decreased evolutionary potential in the long term. The reasons for increasing the evolvability of robust phenotypes include the accumulation of hidden neutral genetic variations that may be useful for later evolutions, the buffering of the pleiotropic side effect of evolution and the increased potential for a neutral exploration of genotype space.^30,32 This is why network robustness is intrinsic to evolution and can improve evolution.^30,32 This view extends these evolutionary results from a genetic level to a gene network level, using the nonlinear stochastic evolution perspective.

Generally, it is still difficult to solve HJI in (40) for the phenotype robustness criterion of the nonlinear stochastic evolutionary gene regulatory network in (38) or (39). In order to determine the extent to which Poisson genetic variations can be tolerated by the nonlinear stochastic evolutionary gene network, or to solve the HJI-constrained optimization in (41) (so as to measure the network evolvability of the nonlinear stochastic gene network so as to gain greater insight into the systematic molecular mechanism of the nonlinear stochastic gene regulatory network in the evolutionary process) the global linearization techniques in (20)–(22) are employed to interpolate several local linear stochastic evolutionary gene networks in (38) or (39) as follows

$$d\tilde{x}\left(t\right)=\sum _{i=1}^M{\alpha }_i\left(\tilde{x}\left(t\right)\right)\left[N_i\tilde{x}\left(t\right)dt+v\left(t\right)dt+\sum _{j=1}^L\Delta N_{ji}\tilde{x}\left(t\right)P\left(t-t_i\right)dt\right]$$

(43) where the interpolation functions α_i ( $\tilde{x}$ (t)) are defined in (22), ie, we use the globally linearized stochastic system in (43) to replace the nonlinear stochastic gene network in (39).

After the nonlinear stochastic evolutionary gene regulatory network in (39) is represented by the interpolation of M linearized gene networks in (43), we obtain the following result of the phenotype robustness criterion for the nonlinear stochastic gene regulatory network in evolution.

Proposition 6: If the following phenotype robustness criterion holds for some P > 0.

$$N_i^{T}P+P{N}_i+PP+\sum _{j=1}^L\Delta N_{ji}^T\Delta N_{ji}+\sum _{j=1}^L\lambda \Delta N_{ji}^{T}P\Delta N_{ji}+\frac{1}{{\rho }^2}{P}^{T}P+I\le 0,\hspace{0.17em}i=1,\dots ,M$$

(44) then the nonlinear stochastic evolutionary gene network (39) is stochastically robustly stable at the equilibrium point $\tilde{x}$ (t) ≡ 0 or x(t) = x_e in spite of Poisson genetic variations and environmental disturbances. In other words, if the phenotype robustness criterion in (44) holds, the Poisson genetic variations could be tolerated and the response to environmental disturbance can be attenuated below an evolution level ρ in the evolutionary process.

Proof: see appendix C

Similarly, the network evolvability of nonlinear stochastic evolutionary gene network can also be measured by solving the following constrained optimization problem

$$\begin{array}{l}{\rho }_0=\underset{P>0}{\min }\rho \hfill \\ \text{subject to}\left(44\right)\hfill \end{array}$$

(45)

The network evolvability ρ₀ of the nonlinear stochastic evolutionary gene network in (45) could be measured by decreasing ρ in the inequalities in (44) until these inequalities have no positive solution P > 0. With the fact that the network fitness f₀ = 1/ρ₀, then the phenotype robustness criterion in (44) can be modified as

$$\underset{local genetic robustness}{\underbrace{PP+\sum _{j=1}^L\Delta N_{ji}^T\Delta N_{ji}+\sum _{j=1}^L{\lambda }_j\Delta N_{ji}^{T}P\Delta N_{ji}}}+\underset{local environmental robustness}{\underbrace{f_0^2{P}^{T}P+I}}\le \underset{local network robustness}{\underbrace{-\left(N_i^{T}P+P{N}_i\right)}},\hspace{1em}i=1,\dots ,M$$

(46)

From the phenotype robustness criterion of the nonlinear stochastic evolutionary gene regulatory network in (46), it can be seen that the phenotype of the nonlinear stochastic evolutionary gene network can be maintained if each local network robustness can confer each local genetic robustness and local environmental robustness so that genetic variations can be tolerated and the environmental disturbances can be resisted in evolution. If the local stochastic evolutionary gene networks in (43) are more robust, ie, the Eigenvalues of local evolutionary gene regulatory networks are farther in the left-half complex s-domain as shown in Figure 1, the nonlinear stochastic evolutionary gene network can provide more network robustness so that it can also provide more genetic robustness and environmental robustness in the evolutionary process. However, the trade-off between local genetic robustness and local environmental robustness is that their total sum cannot be more than the local network robustness. If the phenotype robustness criterion in (46) is violated, the network robustness may not provide enough genetic robustness to simultaneously tolerate genetic variations as well as enough environmental robustness to resist environmental stimuli. Large network evolvability with significant response to environmental disturbances will impel phenotype transition from the basin of equilibrium point x_e to another basin.

Computer Simulation Example

To confirm the validity of the stability robustness and the noise attenuation schemes in gene regulatory network, a computational example in systems biology is shown in the following section. Consider a typical genetic regulatory network, as shown in Figure 3.^18,21,33 This is a typical gene interaction system describing the gene, mRNA, and protein interactions. x₁(t) is an mRNA produced from gene 1, x₂(t) is an enzyme protein produced from x₁(t), and x₃(t) is an inducer protein catalyzed by x₂(t). In addition, x₄(t) is an mRNA produced from gene 4, and x₅(t) is a regulator protein produced from x₄(t). Positive feedback from the inducer protein x₃(t) and negative feedback from the regulator protein x₅(t) are assumed in the mRNA production processes of gene 1 and gene 4. ³⁴ Suppose the genetic regulatory network suffered from some stochastic parameter perturbations and environmental noises can be represented as follows

$$\left[\begin{array}{c}{\dot{x}}_1\left(t\right)\\ {\dot{x}}_2\left(t\right)\\ {\dot{x}}_3\left(t\right)\\ {\dot{x}}_4\left(t\right)\\ {\dot{x}}_5\left(t\right)\end{array}\right]=\left[\begin{array}{l}\left(5+n\left(t\right)\right)x_3\left(t\right)x_5^{-1}\left(t\right)-\left(10+2n\left(t\right)\right)x_1^2\left(t\right)+v_1\left(t\right)\hfill \\ \left(10-2n\left(t\right)\right)x_1^2\left(t\right)-\left(10-0.3n\left(t\right)\right)x_2^2\left(t\right)+v_2\left(t\right)\hfill \\ \left(10+0.8n\left(t\right)\right)x_2^{-1}\left(t\right)-\left(10+1.4\right)x_2^{-1}\left(t\right)x_3^2\left(t\right)+v_3\left(t\right)\hfill \\ \left(8+1.2n\left(t\right)\right)x_3^2\left(t\right)x_5^{-1}\left(t\right)-\left(10+2n\left(t\right)\right)x_4^2\left(t\right)+v_4\left(t\right)\hfill \\ 10x_4^2\left(t\right)-\left(10+1.8n\left(t\right)\right)x_5^2\left(t\right)+v_5\left(t\right)\hfill \end{array}\right]$$

(47) where n(t) denotes the standard white Gaussian noise with zero mean and unit variance. The extrinsic noise vector v(t) = [v₁(t) v₂(t) v₃(t) v₄(t) v₅(t)]^T is the white Gaussian noise with mean = 2 and variance = 1. The initial value x₀ = [0.5 0.5 0.5 0.5 0.5] ^T , and the nominal equilibrium point x_e = [0.7339 0.7339 1 0.9283 0.9283] ^T without n(t) and v(t). The genetic regulatory network under such perturbations and environmental noises then becomes the following nonlinear stochastic system ie, $\tilde{x}$ _e = [0 0 0 0 0] ^T is the stable equilibrium point.

Figure 3.

A typical genetic regulatory network describing the gene, mRNA and protein interactions. ²¹

In order to discuss the robust stability of the stochastic regulatory network at the equilibrium point x_e, following (20), we can globally linearize the system and solve the inequalities in (23) to see if the perturbed stochastic system is stable or not under these stochastic intrinsic and extrinsic noises. If there exists a positive definite P of the phenotype robustness criterion (25), then the stochastic system (49) is stochastically stable at $\tilde{x}$ _e = [0 0 0 0 0] ^T under these stochastic

$$d\left[\begin{array}{l}{x}_1\left(t\right)\hfill \\ x_2\left(t\right)\hfill \\ x_3\left(t\right)\hfill \\ x_4\left(t\right)\hfill \\ x_5\left(t\right)\hfill \end{array}\right]=\left[\begin{array}{l}5x_3\left(t\right)x_5^{-1}\left(t\right)-10x_1^2\left(t\right)+v_1\left(t\right)\hfill \\ 10x_1^2\left(t\right)-10x_2^2\left(t\right)+v_2\left(t\right)\hfill \\ 10x_2^{-1}\left(t\right)-10x_2^{-1}\left(t\right)x_3^2\left(t\right)+v_3\left(t\right)\hfill \\ 8x_3^2\left(t\right)x_5^{-1}\left(t\right)10x_4^2\left(t\right)+v_4\left(t\right)\hfill \\ 10x_4^2\left(t\right)-10x_5^2\left(t\right)+v_5\left(t\right)\hfill \end{array}\right]dt+\left[\begin{array}{l}{x}_3\left(t\right)x_5^{-1}\left(t\right)-2x_1^2\left(t\right)\hfill \\ -2x_1^2\left(t\right)+0.3x_2^2\left(t\right)\hfill \\ 0.8x_2^{-1}\left(t\right)-1.4x_2^{-1}\left(t\right)x_3^2\left(t\right)\hfill \\ 1.2x_3^2\left(t\right)x_5^{-1}\left(t\right)-2x_4^2\left(t\right)\hfill \\ -1.8x_5^2\left(t\right)\hfill \end{array}\right]d\omega \left(t\right)$$

(48) where dω(t) = n(t)dt, and ω(t) is the standard Wiener process (or Brownian motion).

Since we are interested in the robust stability of the equilibrium x_e under stochastic parameter perturbation and environmental noise, the origin should be shifted to the equilibrium point x_e, ie, $\tilde{x}$ ₀ = x₀–x_e = [–0.2339 − 0.2339 − 0.5 − 0.4283 − 0.4283] ^T The genetic regulatory network under coordinate shift should then be rewritten to intrinsic and extrinsic noises. To validate the network robust stability under stochastic intrinsic fluctuations and extrinsic noises, we also want to confirm the network sensitivity ρ₀ in (25) about the estimation of disturbance attenuation of the gene regulatory network.

By solving the constrained optimization in (24), we can find the network sensitivity ρ₀ = 0.7 and a positive P as follows

$$d\left[\begin{array}{l}{\tilde{x}}_1\left(t\right)\hfill \\ {\tilde{x}}_2\left(t\right)\hfill \\ {\tilde{x}}_3\left(t\right)\hfill \\ {\tilde{x}}_4\left(t\right)\hfill \\ {\tilde{x}}_5\left(t\right)\hfill \end{array}\right]=\left[\begin{array}{l}5\left({\tilde{x}}_3\left(t\right)+1\right){\left({\tilde{x}}_5\left(t\right)+0.9283\right)}^{-1}-10{\left({\tilde{x}}_1\left(t\right)+0.7339\right)}^2+v_1\left(t\right)\hfill \\ 10{\left({\tilde{x}}_1\left(t\right)+0.7339\right)}^2-10{\left({\tilde{x}}_2\left(t\right)+0.7339\right)}^2+v_2\left(t\right)\hfill \\ 10{\left({\tilde{x}}_2\left(t\right)+0.7339\right)}^{-1}-10{\left({\tilde{x}}_2\left(t\right)+0.339\right)}^{-1}{\left({\tilde{x}}_3\left(t\right)+1\right)}^2+v_3\left(t\right)\hfill \\ 8{\left({\tilde{x}}_3\left(t\right)+1\right)}^2{\left({\tilde{x}}_5\left(t\right)+0.9283\right)}^{-1}-10{\left({\tilde{x}}_4\left(t\right)+0.9283\right)}^2+v_4\left(t\right)\hfill \\ 10{\left({\tilde{x}}_4\left(t\right)+0.9283\right)}^2-10{\left({\tilde{x}}_5\left(t\right)+0.9283\right)}^2+v_5\hfill \end{array}\right]dt+\left[\begin{array}{l}\left({\tilde{x}}_3\left(t\right)+1\right){\left({\tilde{x}}_5\left(t\right)+0.9283\right)}^{-1}-2{\left({\tilde{x}}_1\left(t\right)+0.7339\right)}^2\hfill \\ -2{\left({\tilde{x}}_1\left(t\right)+0.7339\right)}^2+0.3{\left({\tilde{x}}_2\left(t\right)+0.7339\right)}^2\hfill \\ 0.8{\left({\tilde{x}}_2\left(t\right)+0.7339\right)}^{-1}-1.4{\left({\tilde{x}}_2\left(t\right)+0.7339\right)}^{-1}{\left({\tilde{x}}_3\left(t\right)+1\right)}^2\hfill \\ 1.2{\left({\tilde{x}}_3\left(t\right)+1\right)}^2{\left({\tilde{x}}_5\left(t\right)+0.9283\right)}^{-1}-2{\left({\tilde{x}}_4\left(t\right)+0.9283\right)}^2\hfill \\ -1.8{\left({\tilde{x}}_5\left(t\right)+0.9283\right)}^2\hfill \end{array}\right]d\omega \left(t\right)$$

(49)

$$P=\left[\begin{array}{ccccc}0.38196& 0.076187& -0.024633& -0.1377& -0.033869\\ 0.076187& 0.21315& -0.01921& -0.026898& -0.033226\\ -0.024633& -0.01921& 0.12503& 0.0021901& -0.065039\\ -0.1377& -0.026898& 0.0021901& 0.25956& 0.025702\\ -0.033869& -0.033226& -0.065039& 0.025702& 0.18759\end{array}\right]$$

(50)

This means that the effect of environmental disturbance on the gene regulatory network cannot exceed this value. We can compute the energy ratio of x(t) and v(t) in (4) to verify the network sensitivity ρ₀ = 0.7 obtained from the above computational result. Thus, by Monte Carlo simulation results with 100 runs for the stochastic gene network, we estimate network sensitivity as follows

$$\frac{E{\int }_0^{150}x{\left(t\right)}^{T}x\left(t\right)dt}{E{\int }_0^{150}v{\left(t\right)}^{T}v\left(t\right)dt}\approx {0.6554}^2<{0.7}^2$$

(51)

In this genetic regulatory network we found that the network sensitivity of gene regulatory network can be estimated by the proposed systems biology method and validated by system simulation.

Discussion

One of the most important features of biology is the ability of organisms to persist in face of changing conditions. To achieve this consistently, organisms must have a balance between robustness and evolvability, that is, between resisting and allowing change in their own internal states. ³¹ Moreover, they must achieve this balance on multiple time scales, such as physiological responses to changes over an individual life span as well as evolutionary responses, in which a population of genomes continually update its encoded information about past environments and how future generations should respond given that record. ³¹ There are many examples of robust biological systems found at many scales, from biochemical to ecological. At each scale, robustness may reflect the properties of individual elements or, alternatively, the dynamic feedbacks between interacting elements. The expression of some metabolic functions may be robust in face of temperature changes. For example, an enzyme maintains its shape and specificity across a range of temperatures or because an interconnected network of reactions sustains the supply of product, even when some enzymes fail. A genome may be robust, on the other hand, because it encodes proofreading and repairs systems that reduce replication errors or because it is organized such that many mutations have little effect on its phenotype.

One important question is whether there exists a unifying mathematical framework that can encompass such diverse examples of biological robustness under intrinsic perturbations and extrinsic disturbances. Might new insights come from such conceptual and mathematical unification, or will future understanding require detailed analyses of specific cases? Across the different biological system scales, recurring mechanisms for achieving robustness, which include redundancy, modularity and feedback, might serve as organizing principles of robust biological networks. Yet, similar robust mechanism could mask important differences in the evolutionary origins of those robust mechanisms. At the level of genes in genomes or of cells in multicellular organisms, it is reasonable to suggest that redundancy evolved by natural selection in order to maintain some functional capacity in face of intrinsic perturbations and extrinsic disturbances. However, while species redundancy could also be critical for robustness of ecosystem functions, differences in redundancy might be an emergent property rather than an ecosystem-level adaptation as selection generally acts at lower levels. ³¹ If robustness has evolved to maintain performance, what would prevent biological networks from becoming ever more robust? To answer these questions, in this study we have focused on genetic, environmental, and network robustness, and the interplay between them to discuss the phenotype robustness of biological networks at different scales from a unifying nonlinear stochastic system perspective. We found that if a biological network becomes more robust in the evolutionary process, a robust phenotype of biological network will harbor a large amount of neutral genetic variations, which might show increased rather than decreased evolutionary potential in the long run. This is the reason why network robustness is essential to the evolution of biological network and why it can improve the network evolution in biological systems.

In this study, we have developed a single unifying mathematical framework for encompassing diverse examples of stochastic biological networks to discuss intrinsic, genetic, environmental, and network robustness, and their trade-offs in systems and evolutionary biology.

According to our analyses, the phenotype robustness criteria of stochastic biological networks in systems and evolutionary biology have a similar mathematical framework. The biological networks in systems and evolutionary biology can be modeled as nonlinear stochastic systems with intrinsic parameter fluctuations, genetic variations and environmental disturbances, in which intrinsic parameter fluctuations are described by the Wiener (Brownian) process, environmental disturbances are described by the Gaussian white noise, and genetic variations are described by the Poisson point process. The interplay between these four different areas of robustness can therefore be analyzed by the nonlinear stochastic system theory. The linear stochastic system theory can then be applied when the global linearization technique is employed to interpolate several local linear stochastic systems to approximate the nonlinear stochastic system.

From the system theory perspective, the phenotype robustness of nonlinear stochastic gene networks in systems and evolutionary biology need to obey a similar phenotype robustness criterion. In order words, “intrinsic robustness + genetic robustness + environmental robustness ≤ network robustness.” This means network robustness needs to be strong enough to tolerate either heritable perturbations (genetic variations) or non-heritable perturbations (ie, random molecular fluctuations and environmental disturbances) in order for the phenotype of gene networks to be maintained in systems and evolutionary biology with a similar mathematical framework. The phenotype robustness of the stochastic gene network is completely consistent with the idea of canalization of development and inheritance of acquired characters as described by Waddington. ³⁵ According to these phenotype robustness criteria, the correlation among intrinsic, genetic, environmental, and network robustness by recent genomic experiments in yeast (genes conferring similar intrinsic, genetic and environmental robustness to maintain phenotypic robustness) can be rationally explained from the systematic perspective. ³⁶ In other words, if the network robustness of gene network is large enough, genetic perturbations or environmental disturbances can then be taken over respectively or simultaneously to maintain the functional phenotype in systems and evolutionary biology.

Genetic, environmental and phenotypic random variations are inevitably noisy processes in systems and evolutionary biology rather than desirable features of biological networks.^37,38 These stochastic processes arise from the complexity of the evolutionary process of biological gene networks. However, there is numerous evidence of high fidelity and minimal noise, including the proof editing of DNA replication and protein translation in systems and evolutionary biology. Enzymes have also evolved toward high specificity thereby increasing fidelity. Gene expression is also regulated by elaborate mechanisms, and random variations seem to have been minimized in systems, synthetic, and evolutionary biology. However, chemophysical constraints on the specificity and fidelity of biological networks are costly and there are generally trade-offs. If biological networks want to retain enough network robustness to give intrinsic robustness for tolerating intrinsic parameter fluctuations, genetic robustness for buffering genetic variations and environmental robustness for resisting environmental disturbances, so as to keep their proper function (ie, phenotype robustness), much effort has to be taken and a high cost must be paid.

In general, random genetic variations, phenotype perturbations, and heterogeneity are neither desired nor deliberate outcomes of systems and evolutionary biology. However, heterogeneity and diversity form the very basis of evolutionary biology, not only within genetically diverse populations but also within the same allele or genome. Thus, random genetic variations, environmental disturbances, and phenotypic perturbations are inherent features of biological systems and networks. Random perturbative biological networks may contain more connected and interconnected systems and may provide multifunctionality of the biological network. This multifunctionality may result in increased robustness and a capacity to cope with diverse challenges. However, multifunctionality also increases the complexity and the variations in the biological network, which may increase adaptive potential. Thus, behind the façade of perfection and optimality of systems and evolutionary biology lies the messy biology that originates from the genetic variations and environmental disturbances in evolution. ³⁹ There exist the trade-off among intrinsic, genetic, environmental, and network robustness in the phenotype robustness of stochastic biological networks. That is, if intrinsic robustness + genetic robustness + environmental robust ≤ network robustness, then the phenotype of biological network is maintained. This sheds light on the mechanisms that govern the exploitation and toleration of the messiness of biological networks in systems and evolutionary biology, from the systematic perspective. Obviously, network robustness needs to be strong enough to tolerate either heritable perturbations (genetic variations) or non-heritable perturbations (random molecular fluctuations and environmental disturbances) so that the phenotype can be maintained in biological network at different levels.

The interplay between evolvability and network robustness in evolutionary network has been discussed by Chen and Lin. ¹⁴ However, in this study, the Wiener (or Brownian) processes for modeling genetic variations in evolutionary gene networks have been replaced by the Poisson point processes to better mimic the discontinuous genetic mutations of gene networks in the evolutionary process. Further, some genetic algorithms (GAs) based on genetic mutations and natural selection in the evolutionary process have been widely applied to both control engineering design to select the most adequate controller^40,41 and genetic circuit designs in synthetic biology to select the most adequate circuit components^42,43 to satisfy the prescribed design specification respectively.

In general, it is very difficult to solve the HJI in (14), (40) for the phenotype robustness criteria in biological networks at different levels. With the global linearization technique, the HJI problem for robust stabilization of nonlinear stochastic biological network is reduced to solving an equivalent set of Riccati-like inequalities in (23), (44) for the robust stabilization of each local linearized biological network. We also found that if the network robustness of each local linearized biological network can take on the local intrinsic robustness, genetic robustness and environmental robustness of each local linearized biological network, then the phenotype of the nonlinear stochastic biological network could also be maintained.

For biological networks at different levels, two favored strategies can improve phenotype robustness in the evolutionary process. One is to improve network robustness to provide enough intrinsic robustness for tolerating intrinsic parameter fluctuations, genetic robustness for buffering genetic variations, and environmental robustness for resisting environmental disturbances so that phenotype robustness of the biological network can be maintained under these uncertain perturbations and environmental disturbances. Negative feedback is a mechanism that can improve network robustness (ie, it can make right-hand sides of (17), (19), (25), (27), (42) or (46) larger) and is favored by natural selection in biological networks at different levels in the evolutionary process. Another strategy is to reduce the effect of intrinsic parameter fluctuations, genetic variations and environmental disturbances on different biological networks (ie, it can make the left-side of (17), (19), (25), (27), (42) or (46) smaller). Redundancies and repairs are the mechanisms of this strategy and are favored by natural selection in evolution. This is the reason why there are so many redundancies from duplicated genes in networks and redundant pathways in biochemical networks.

Conclusion

This paper presents a unifying mathematical framework to describe different levels of stochastic biological networks under intrinsic parameter fluctuations, genetic variations, and environmental disturbances. The phenotype robustness criteria of biological networks in systems and evolutionary biology were also investigated, according to the unifying stochastic biological systems, from the robust stabilization and disturbance sensitivity perspective. It was found that if intrinsic robustness + genetic robustness + environmental robustness ≤ network robustness (ie, network robustness can confer intrinsic robustness for tolerating intrinsic parameter fluctuations, genetic robustness for buffering genetic variations and environmental robustness for resisting the environmental disturbances) then the phenotype will be robust in biological networks at different levels in systems and evolutionary biology. Using the global linearization method, we also found that if the network robustness of each local linearized system is greater than the total sum of intrinsic robustness, genetic robustness, and environmental robustness of each local linear system, then the phenotype of the biological network is also maintained despite intrinsic parameter fluctuations, genetic variations, and environmental disturbances. Finally, an example in silico is given to estimate the network sensitivity of a gene regulatory network, which can be also validated by Monte Carlo simulation.