The splitting approach is developed for the numerical simulation of genetic regulatory networks with a stable steady-state structure. The numerical results of the simulation of a one-gene network, a two-gene network, and a p53-mdm2 network show that the new splitting methods constructed in this paper are remarkably more effective and more suitable for long-term computation with large steps than the traditional general-purpose Runge-Kutta methods. The new methods have no restriction on the choice of stepsize due to their infinitely large stability regions.

1. Introduction

The exploration of mechanisms of gene expression and regulation has become one of the central themes in medicine and biological sciences such as cell biology, molecular biology, and systems biology [1, 2]. For example, it has been acknowledged that the p53 tumor suppressor plays key regulatory roles in various fundamental biological processes, including development, ageing, and cell differentiation. It can regulate its downstream genes through their signal pathways and further implement cell cycle arrest and cell apoptosis [36]. The qualitative analysis as well as numerical simulation has become an important route in the investigation of differential equations of genetic regulatory networks (GRNs) in the past few years [710]. Up till now, algorithms used in the simulation of GRNs have primarily been classical Runge-Kutta (RK) methods (typically of order four) or Runge-Kutta-Fehlberg embedded pairs as employed in the scientific computing software MATLAB [1113]. However, if we are required to achieve a very high accuracy, we have to take very small stepsize. Moreover, the traditional Runge-Kutta type methods often fail to retain some important qualitative properties of the system of interest. This prevents us from acquiring correct knowledge of the dynamics of genetic regulatory networks.

Geometric numerical integration aims at solving differential equations effectively while preserving the geometric properties of the exact flow [14]. Recently, You et al. [15] develop a family of trigonometrically fitted Scheifele two-step (TFSTS) methods, derive a set of necessary and sufficient conditions for TFSTS methods to be of up to order five based on the linear operator theory, and construct two practical methods of algebraic four and five, respectively. Very recently, You [16] develops a new family of phase-fitted and amplification methods of Runge-Kutta type which have been proved very effective for genetic regulatory networks with a limit-cycle structure.

Splitting is one of the effective techniques in geometric integration. For example, Blanes and Moan [17] construct a symmetric fourth- and sixth-order symplectic partitioned Runge-Kutta and Runge-Kutta-Nyström methods and show that these methods have an optimized efficiency. For a systematic presentation of the splitting technique, the reader is referred to Hairer et al. [14]. The purpose of this paper is to develop the splitting methods for GRNs. In Section 2 we present the system of differential equations governing the GRNs and basic assumptions for the system. In Section 3 we describe the idea and formation of the approach of splitting strategy which intends to simulate exactly the characteristic part of the system. Section 4 gives the simulation results of the new splitting methods and the traditional Runge-Kutta methods when they are applied to a one-gene network, a two-gene network, and a p53-mdm2 network. We compare their accuracy and computational efficiency. Section 5 is devoted to conclusive remarks. Section 6 is for discussions. In Appendix, the linear stability of the new splitting methods is analyzed.

2. Materials

2.1. mRNA-Protein Networks

An -gene regulatory network can be described by the following system of ordinary differential equations: where and are -dimensional vectors representing the concentrations of mRNAs and proteins at time , respectively, and , , , and are diagonal matrices. The system (1) can be written in components as where, for , and are the concentrations of mRNA and the corresponding protein at time , respectively; and , positive constants, are the degradation rates of mRNA and protein , respectively; is a positive constant; and , the regulatory function of gene , is a nonlinear and monotonic function in each of its variables. If gene is activated by protein , , and if gene is inhibited by protein , . If protein has no action on gene , will not appear in the expression of .

In particular, we are concerned in this paper with the following two simple models.

(I) The first model is a one-gene regulatory network which can be written as where represents the action of an inhibitory protein that acts as a dimer and , and are positive constants. This model with delays can be found in Xiao and Cao [18].

(II) The second model is a two-gene cross-regulatory network [7, 19]: where and are the concentrations of mRNA and mRNA , respectively, and are the concentrations of their corresponding products protein and protein , respectively, , and represent the maximal transcription rates of gene and gene , respectively, and are the degradation rates of mRNA and mRNA , respectively, and are the degradation rates of protein and protein , respectively, are the Hill functions for activation and repression, respectively, and are the Hill coefficients, and and are the thresholds. It is easy to see that the activation function is increasing in and the repression function is decreasing in .

2.2. A p53-mdm2 Regulatory Pathway

Another model we are interested in is for the damped oscillation of the p53-mdm2 regulatory pathway which is given by (see [20]) where represents the concentration of the p53 tumour suppressor, (mdm2) is the concentration of the p53's main negative regulator, is the concentration of the p53-mdm2 complex, is the concentration of an active form of p53 that is resistant against mdm2-mediated degradation, is a transient stress stimulus which has the form ,  , are de novo synthesis rates, are production rates, are reverse reactions (e.g., dephosphorylation), is the degradation rate of active p53, and is the saturation coefficient.

3. Methods

3.1. Runge-Kutta Methods

Either the mRNA-protein network (1) or the p53-mdm2 regulatory pathway (6) can be regarded as a special form of a system of ordinary differential equations (ODEs): where and the function is smooth enough as required. The system (7) together with initial value is called an initial value problem (IVP). Throughout this paper we make the following assumptions. (i)The system (7) has a unique positive steady state ; that is, there is a unique with for such that .(ii)The steady state is asymptotically stable; that is, for any solution of the system (7) through a positive initial point , .

The most frequently used algorithms for the system (7) are the so-called Runge-Kutta methods which read where is an approximation of the solution at ,  , ,  , are real numbers, is the number of internal stages , and is the step size. The scheme (8) can be represented by the Butcher tableau: or simply by , where for . Two of the most famous fourth-order RK methods have the tableaux as follows (see [13]):xy(10)

which we denote as RK4 and RK3/8, respectively.

3.2. Splitting Methods

Splitting methods have been proved to be an effective approach to solve ODEs. The main idea is to split the vector field into two or more integrable parts and treat them separately. For a concise account of splitting methods, see Chapter II of Hairer et al. [14].

Suppose that the vector field of the system (7) has a split structure Assume also that both systems and can be solved in closed form or are accurately integrated and their exact flows are denoted by and , respectively.

Definition 1. (i) The method defined by is the simplest splitting method for the system (7) based on the decomposition (11) (see [13]).
(ii) The Strang splitting is the following symmetric version (see [21]):
(iii) The general splitting method has the form where are positive constants satisfying some appropriate conditions. See, for example, [2225].

Theorem 5.6 in Chapter II of Hairer et al. [14] gives the conditions for the splitting method (14) to be of order .

However, in most occasions, the exact flows and for and in Definition 1 are not available. Hence, we have to use instead some approximations or numerical flows which are denoted by and .

3.3. Splitting Methods for Genetic Regulatory Networks Based on Their Characteristic Structure

For a given genetic regulatory network, different ways of decomposition of the vector field may produce different results of computation. Thus a question arises as follows: which decomposition is more appropriate or more effective. In the following we take the system (1) for example. The analysis of the p53-mdm2 pathway (6) is similar. Denote . Then the -gene regulatory network (1) has a natural form of decomposition: Unfortunately, it has been checked through practical test that the splitting methods based on this decomposition cannot lead to effective results. To find a way out, we first observe that a coordinate transform translates the steady state of the system (1) to the origin and yields where is the Jacobian matrix of at point and   −  . We employ this special structure of the system (16) to reach the decomposition of the vector field: where . The system here is called the linearization of the system (1) at the steady state . in (17) is the linear principal part of the vector field which has the exact flow . However, it is not easy or impossible to obtain the exact solution of due to its nonlinearity. So we have to use an approximation flow and form the splitting method: When is taken as an RK method, then the resulting splitting method is denoted by Split(Exact:RK). Hence we write Split(Exact:RK4) and Split(Exact:RK3/8) for the splitting methods with taken as RK4 and RK3/8, respectively.

4. Results

In order to examine the numerical behavior of the new splitting methods Split(Exact:RK4) and Split(Exact:RK3/8), we apply them to the three models presented in Section 2. Their corresponding RK methods RK4 and RK3/8 are also used for comparison. We will carry out two observations: effectiveness and efficiency. For effectiveness, we first find the errors produced by each method with different values of stepsize. We also solve each problem with a fixed stepsize on different lengths of time intervals.

4.1. The One-Gene Network

Table 1 gives the parameter values which are provided by Xiao and Cao [18]. This system has a unique steady state where the Jacobian matrix has eigenvalues ,  , where is the imaginary unit satisfying . Since the two eigenvalues both have negative real parts, the steady state is asymptotically stable.

In order to apply the splitting methods Split(Exact:RK4) and Split(Exact:RK3/8), the vector field of the system (3) is decomposed in the way of (16) as

The system is solved on the time interval with initial values of mRNA and protein ,   and with different stepsizes. The numerical results are presented in Table 2.

Then we solve the problem with a fixed stepsize on several lengths of time intervals. The numerical results are given in Table 3.

4.2. The Two-Gene Network

Table 4 gives the parameter values which can be found in Widder et al. [19]. This system has a unique steady state . Since the four eigenvalues ,  ,  , and of the Jacobian matrix at the steady state all have negative real parts, the steady state is asymptotically stable.

For this system, the decomposition (16) becomes

The system is solved on the time interval with the initial values ,  ,  , and and with different stepsizes. The numerical results are presented in Table 5.

Then we solve the problem with a fixed stepsize on several lengths of time intervals. The numerical results are given in Table 6.

4.3. The p53-mdm2 Network

Table 7 gives the parameter values which are used by van Leeuwen et al. [20]. For simplicity, we take the small function . This system has a unique steady state   = . Since the eigenvalues ,  ,  , and of the Jacobian matrix at the steady state all have negative real parts, the steady state is asymptotically stable.

For this system, decomposition (16) becomes where

The system is solved on the time interval with initial values  nM,  nM,  nM, and  nM and with different stepsizes. The numerical results are presented in Table 8.

Then we solve the problem with a fixed stepsize on several lengths of time intervals. The numerical results are given in Table 9.

5. Conclusions

In this paper we have developed a new type of splitting algorithms for the simulation of genetic regulatory networks. The splitting technique has taken into account the special structure of the linearizing decomposition of the vector field. From the results of numerical simulation of Tables 2, 5, and 8, we can see that the new splitting methods Split(Exact:RK4) and Split(Exact:RK3/8) are much more accurate than the traditional Runge-Kutta methods RK4 and RK3/8. For large steps when RK4 and RK3/8 completely lose effect, Split(Exact:RK4) and Split(Exact:RK3/8) continue to work very well. On the other hand, Tables 3, 6, and 9 show that for comparatively large steps, RK4 and RK3/8 can solve the problem only on short time intervals while Split(Exact:RK4) and Split(Exact:RK3/8) work for very long time intervals.

We conclude that, for genetic regulatory networks with an asymptotically stable steady state, compared with the traditional Runge-Kutta, the new splitting methods have two advantages.(a)They are extremely accurate for large steps. This promises high efficiency for solving large-scale systems (complex networks containing a large number of distinct proteins) in a long-term simulation.(b)They work effectively for very long time intervals. This makes it possible for us to explore the long-run behavior of genetic regulatory network which is important in the research of gene repair and gene therapy.

The special structure of the new splitting methods and their remarkable stability property (see Appendix) are responsible for the previous two advantages.

6. Discussions

The splitting methods designated in this paper have opened a novel approach to effective simulation of the complex dynamical behaviors of genetic regulatory network with a characteristic structure. It is still possible to enhance the effectiveness of the new splitting methods. For example, higher-order splitting methods can be obtained by recursive composition (14) or by employing higher order Runge-Kutta methods; see II.5 of [13]. Another possibility is to consider embedded pairs of two splitting methods which can improve the efficiency; see II.4 of [13].

The genetic regulatory networks considered in this paper are nonstiff. For stiff systems (whose Jacobian possesses eigenvalues with large negative real parts or with purely imaginary eigenvalues of large modulus), the previous techniques suggested by the reviewer are applicable. Moreover, the error control technique which can increase the efficiency of the methods is an interesting theme for future work.

There are more qualitative properties of the genetic regulatory networks that can be taken into account in the designation of simulation algorithms. For example, oscillation in protein levels is observed in most regulatory networks. Symmetric and symplectic methods have been shown to have excellent numerical behavior in the long-term integration of oscillatory systems even if they are not Hamiltonian systems. A brief account of symmetric and symplectic extended Runge-Kutta-Nyström (ERKN) methods for oscillatory Hamiltonian systems and two-step ERKN methods can be found, for instance, in Yang et al. [26], Chen et al. [27], Li et al. [28], and You et al. [29].

Finally, a problem related to this work remains open. We observe that, in Tables 3 and 9 for the p53-mdm2 pathway, as the time interval extends, the error produced by Split(Exact:RK4) and Split(Exact:RK3/8) becomes even smaller. This phenomenon is yet to be explained.


Stability Analysis of Runge-Kutta Methods and Splitting Methods

Stability analysis is a necessary step for a new numerical method before it is put into practice. Numerically unstable methods are completely useless. In this appendix, we examine the linear stability of the new splitting method constructed in Section 3. To this end, we consider the linear scalar test equation: where is a test rate which is an estimate of the principal rate of a scalar problem, and is the error of the estimation. Applying the splitting method (18) to the test equation (A.1), we obtain the difference equation where is called the stability function with and .

Definition A.1. The region in the plane is called the stability region of the method and the curve defined by is call the boundary of stability region.

By simple computation, we obtain the stability function of an RK method where , is the unit matrix and the stability function a splitting method Split(Exact:RK)

The stability regions of RK4 and RK3/8 are presented in Figure 1 and those of Split(Exact:RK4) and Split(Exact:RK3/8) are presented in Figure 2. It is seen that Split(Exact:RK4) and Split(Exact:RK3/8) have much larger stability regions than RK4 and RK3/8. Moreover, the infinity area of the stability regions of Split(Exact:RK4) and Split(Exact:RK3/8) means that these two methods have no limitation to the stepsize for the stability reason but for the consideration of accuracy, while RK4 and RK3/8 will completely lose effect when the stepsize becomes large to some extent. This has been verified in Section 4.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The authors are grateful to the anonymous referees for their invaluable comments and constructive suggestions which help greatly to improve this paper. This work was partially supported by NSFC (Grant no. 11171155) and the Fundamental Research Fund for the Central Universities (Grant no. Y0201100265).