Splitting strategy for simulating genetic regulatory networks.

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 [3-6]. 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 [7-10]. 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 [11-13]. 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 N-gene regulatory network can be described by the following system of ordinary differential equations: where r(t) = (r 1(t),…, r (t)) and p(t) = (p 1(t), p 2(t),…, p (t)) are N-dimensional vectors representing the concentrations of mRNAs and proteins at time t, respectively, and F(p(t)) = (F 1(p(t)),…, F (p(t))), Γ = diag⁡(γ 1,…, γ ), M = diag⁡(μ 1,…, μ ), and K = diag⁡(κ 1,…, κ ) are diagonal matrices. The system (1) can be written in components as where, for i = 1,2,…, N, r (t) and p (t) ∈ ℝ+ are the concentrations of mRNA i and the corresponding protein i at time t, respectively; γ and μ , positive constants, are the degradation rates of mRNA i and protein i, respectively; κ is a positive constant; and F (·), the regulatory function of gene i, is a nonlinear and monotonic function in each of its variables. If gene i is activated by protein j, ∂F /∂p > 0, and if gene i is inhibited by protein j, ∂F /∂p < 0. If protein j has no action on gene i, p will not appear in the expression of F .

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 f(p(t)) = α/(1 + p(t)2/θ 2) 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 r 1 and r 2 are the concentrations of mRNA 1 and mRNA 2, respectively, p 1 and p 2 are the concentrations of their corresponding products protein 1 and protein 2, respectively, λ 1, and λ 2 represent the maximal transcription rates of gene 1 and gene 2, respectively, γ 1 and γ 2 are the degradation rates of mRNA 1 and mRNA 2, respectively, δ 1 and δ 2 are the degradation rates of protein 1 and protein 2, respectively, are the Hill functions for activation and repression, respectively, n 1 and n 2 are the Hill coefficients, and θ 1 and θ 2 are the thresholds. It is easy to see that the activation function h + is increasing in p 2 and the repression function h − is decreasing in p 1.

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 P represents the concentration of the p53 tumour suppressor, M (mdm2) is the concentration of the p53's main negative regulator, C is the concentration of the p53-mdm2 complex, P is the concentration of an active form of p53 that is resistant against mdm2-mediated degradation, S(t) is a transient stress stimulus which has the form S(t) = −e ,  c = γk , s ∗(∗ = p, m0, m1) are de novo synthesis rates, k ∗(∗ = a, c, u) are production rates, j ∗(∗ = a, c) are reverse reactions (e.g., dephosphorylation), d is the degradation rate of active p53, and K 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 z = z(t) ∈ R and the function g : ℝ → ℝ is smooth enough as required. The system (7) together with initial value z(0) = z 0 is called an initial value problem (IVP). Throughout this paper we make the following assumptions.

The system (7) has a unique positive steady state z*; that is, there is a unique z* = (z 1*, z 2*,…, z *) with z * > 0 for i = 1,2,…, d such that g(z*) = 0.

The steady state z* is asymptotically stable; that is, for any solution z(t) of the system (7) through a positive initial point z 0, lim⁡⁡z(t) = z*.

The most frequently used algorithms for the system (7) are the so-called Runge-Kutta methods which read where z is an approximation of the solution z(t) at t ,  n = 0,1,…, a , b ,  i, j = 1,…, s, are real numbers, s is the number of internal stages Z , and h is the step size. The scheme (8) can be represented by the Butcher tableau: or simply by (c, A, b(ν)), where c = ∑ a for i = 1,…, s. Two of the most famous fourth-order RK methods have the tableaux as follows (see [13]):

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 g of the system (7) has a split structure Assume also that both systems z′ = g [1](z) and z′ = g [2](z) can be solved in closed form or are accurately integrated and their exact flows are denoted by φ [1] and φ [2], 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 a 1, b 1, a 2, b 2,…, a , b are positive constants satisfying some appropriate conditions. See, for example, [22-25].

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

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

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 f 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 z(t) = (r(t), p(t)). Then the N-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 x(t) = r(t) − r*, y(t) = p(t) − p* translates the steady state (r*, p*) of the system (1) to the origin and yields where F′(p*) is the Jacobian matrix of F(p) at point p* and G(y(t)) = F(p* + y(t)) − F′(p*)y(t)  −  F(p*). We employ this special structure of the system (16) to reach the decomposition of the vector field: where z(t) = (x(t), y(t)). The system here is called the linearization of the system (1) at the steady state (r*, p*). g [1] in (17) is the linear principal part of the vector field g which has the exact flow . However, it is not easy or impossible to obtain the exact solution of g [2](x (t), y (t)) due to its nonlinearity. So we have to use an approximation flow ψ [2] and form the splitting method: When ψ [2] 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 ψ [2] 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 (r*, p*) = (0.6,2) where the Jacobian matrix has eigenvalues ω 1 = −1.2500 + 2.9767i,  ω 2 = −1.2500 − 2.9767i, where i is the imaginary unit satisfying i 2 = −1. Since the two eigenvalues both have negative real parts, the steady state is asymptotically stable.

Table 1

Parameter values for the one-gene network.

α = 3

γ = 1

μ = 1.5

κ = 5

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 [0,100] with initial values of mRNA and protein r(0) = 0.6,  p(0) = 0.8 and with different stepsizes. The numerical results are presented in Table 2.

Table 2

One-gene network: average errors for different stepsizes.

Stepsize	RK4	RK3/8	Split(Exact:RK4)	Split(Exact:RK3/8)
1.2	1.7733 × 100	6.2404 × 10−1	1.3910 × 10−3	1.3910 × 10−3
1.5	3.2171 × 1017	1.2337 × 101	1.1862 × 10−3	1.1862 × 10−3
2.0	3.3619 × 1059	6.0455 × 1047	5.4651 × 10−4	5.4651 × 10−4
10.0	8.7073 × 1062	7.6862 × 1062	1.1451 × 10−7	1.1451 × 10−7

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

Table 3

One-gene network: average errors for fixed stepsize h = 2 on different time intervals.

Time interval	RK4	RK3/8	Split(Exact:RK4)	Split(Exact:RK3/8)
[0,100]	3.3619 × 1059	6.0455 × 1047	5.4917 × 10−4	5.4917 × 10−4
[0,500]	4.5549 × 10299	2.0782 × 10240	1.1158 × 10−4	1.1158 × 10−4
[0,1000]	NaN	NaN	5.5903 × 10−5	5.5903 × 10−5
[0,1500]	NaN	NaN	3.7294 × 10−5	3.7294 × 10−5

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 (r 1*, r 2*, p 1*, p 2*) = (0.814713,0.614032,0.814713,0.614032). Since the four eigenvalues ω 1 = −1.9611 + 0.9611i,  ω 2 = −1.9611 − 0.9611i,  ω 3 = −0.0389 + 0.9611i, and ω 4 = −0.0389 − 0.9611i of the Jacobian matrix at the steady state (r*, p*) all have negative real parts, the steady state is asymptotically stable.

Table 4

Parameter values for the two-gene network.

λ 1 = 1.8	γ 1 = 1	κ 1 = 1	δ 1 = 1	θ 1 = 0.6542	n 1 = 3
λ 2 = 1.8	γ 2 = 1	κ 2 = 1	δ 2 = 1	θ 1 = 0.6542	n 2 = 3

For this system, the decomposition (16) becomes

The system is solved on the time interval [0,100] with the initial values r 1(0) = 0.6,  r 2(0) = 0.8,  p 1(0) = 0.6, and p 2(0) = 0.8 and with different stepsizes. The numerical results are presented in Table 5.

Table 5

Two-gene network: average errors for different stepsizes.

Stepsize	RK4	RK3/8	Split(Exact:RK4)	Split(Exact:RK3/8)
0.1	9.8188 × 10−2	9.3019 × 10−2	9.9338 × 10−4	9.9338 × 10−4
1.2	2.5157 × 10−1	5.5798 × 10−2	7.6803 × 10−3	7.6803 × 10−3
1.5	4.3953 × 100	2.9958 × 100	9.5832 × 10−3	9.5832 × 10−3
2.0	4.1821 × 1024	7.5238 × 1017	1.6686 × 10−2	1.6686 × 10−2
5.0	1.9692 × 1069	3.2225 × 1068	3.1227 × 10−2	3.1227 × 10−2

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

Table 6

Two-gene network: average errors for fixed stepsize h = 2 on different time intervals.

Time interval	RK4	RK3/8	Split(Exact:RK4)	Split(Exact:RK3/8)
[0,500]	4.4414 × 100	1.5342 × 100	1.9352 × 10−3	1.9352 × 10−3
[0,1000]	4.4396 × 100	1.0172 × 100	9.6936 × 10−4	9.6936 × 10−4
[0,1500]	NaN	2.3844 × 1095	1.1409 × 10−3	1.1409 × 10−3
[0,2000]	NaN	NaN	8.5613 × 10−4	8.5613 × 10−4
[0,2500]	NaN	NaN	6.8520 × 10−4	6.8520 × 10−4

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 S(t) ≡ 0. This system has a unique steady state (P *, M*, C*, P *  = (9.42094,0.0372868,3.49529,0). Since the eigenvalues ω 1 = −38.4766,  ω 2 = −0.0028 + 0.0220i,  ω 3 = −0.0028 − 0.0220i, and ω 4 = −0.2002 of the Jacobian matrix at the steady state all have negative real parts, the steady state is asymptotically stable.

Table 7

Parameter values for the p53-mdm2 pathway.

s m0 = 2 × 10−3 nM min−1	k a = 20 min−1	j a = 0.2 min−1	γ = 2.5
s m1 = 0.15 nM min−1	k c = 4 min−1 nM−1	j c = 2 × 10−3 min−1
s m2 = 0.2 nM min−1	k u = 0.4 min  −1	d m = 0.4 min−1
s p = 1.4 nM min−1	K m = 100 nM	d p = 2 × 10−4 min−1

For this system, decomposition (16) becomes where

The system is solved on the time interval [0,100] with initial values P (0) = 9.42 nM, C(0) = 0.037 nM, M(0) = 3.49 nM, and P (0) = 0 nM and with different stepsizes. The numerical results are presented in Table 8.

Table 8

p53-mdm2 network: average errors for different stepsizes.

Stepsize	RK4	RK3/8	Split(Exact:RK4)	Split(Exact:RK3/8)
0.05	3.7686 × 10−5	3.7046 × 10−5	1.2091 × 10−6	1.2065 × 10−6
0.08	4.6118 × 100	4.5057 × 100	2.1383 × 10−6	2.1290 × 10−6
0.10	NaN	5.3550 × 100	2.8098 × 10−6	2.7945 × 10−6
0.12	NaN	NaN	3.4986 × 10−6	3.4762 × 10−6
5.00	NaN	NaN	7.3001 × 10−4	3.8208 × 10−4

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

Table 9

p53-mdm2 network: average errors for fixed stepsize h = 10 on different time intervals.

Time interval	RK4	RK3/8	Split(Exact:RK4)	Split(Exact:RK3/8)
[0,100]	NaN	NaN	6.8810 × 10−2	6.6940 × 10−2
[0,500]	NaN	NaN	2.1634 × 10−2	2.0911 × 10−2
[0,1000]	NaN	NaN	1.1625 × 10−2	1.1214 × 10−2
[0,1500]	NaN	NaN	7.8314 × 10−3	7.5529 × 10−3
[0,2000]	NaN	NaN	5.8881 × 10−3	5.6786 × 10−3

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.

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.

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.