Splitting Strategy for Simulating Genetic Regulatory Networks

You, Xiong; Liu, Xueping; Musa, Ibrahim Hussein

doi:https://doi.org/10.1155/2014/683235

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Materials Methods Results Conclusions Appendix Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2014 | Article ID 683235 | https://doi.org/10.1155/2014/683235

Splitting Strategy for Simulating Genetic Regulatory Networks

Xiong You,¹Xueping Liu,¹and Ibrahim Hussein Musa¹

Academic Editor: Damien Hall

Received30 May 2013

Accepted24 Oct 2013

Published02 Feb 2014

Abstract

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 -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]):(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, [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 .

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.

Appendix

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.

(a)

(b)

(a)

(b)

Conflict of Interests

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

Acknowledgments

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

References

J. Cao and F. Ren, “Exponential stability of discrete-time genetic regulatory networks with delays,” IEEE Transactions on Neural Networks, vol. 19, no. 3, pp. 520–523, 2008.
View at: Publisher Site | Google Scholar
F. Ren and J. Cao, “Asymptotic and robust stability of genetic regulatory networks with time-varying delays,” Neurocomputing, vol. 71, no. 4–6, pp. 834–842, 2008.
View at: Publisher Site | Google Scholar
J. J. Chou, H. Li, G. S. Salvesen, J. Yuan, and G. Wagner, “Solution structure of BID, an intracellular amplifier of apoptotic signaling,” Cell, vol. 96, no. 5, pp. 615–624, 1999.
View at: Google Scholar
C. A. Perez and J. A. Purdy, “Treatment planning in radiation oncology and impact on outcome of therapy,” Rays, vol. 23, no. 3, pp. 385–426, 1998.
View at: Google Scholar
B. Vogelstein, D. Lane, and A. J. Levine, “Surfing the p53 network,” Nature, vol. 408, no. 6810, pp. 307–310, 2000.
View at: Publisher Site | Google Scholar
J. P. Qi, S. H. Shao, D. D. Li, and G. P. Zhou, “A dynamic model for the p53 stress response networks under ion radiation,” Amino Acids, vol. 33, no. 1, pp. 75–83, 2007.
View at: Publisher Site | Google Scholar
A. Polynikis, S. J. Hogan, and M. di Bernardo, “Comparing different ODE modelling approaches for gene regulatory networks,” Journal of Theoretical Biology, vol. 261, no. 4, pp. 511–530, 2009.
View at: Publisher Site | Google Scholar
A. Goldbeter, “A model for circadian oscillations in the Drosophila period protein (PER),” Proceedings of the Royal Society B, vol. 261, no. 1362, pp. 319–324, 1995.
View at: Publisher Site | Google Scholar
C. Gérard, D. Gonze, and A. Goldbeter, “Dependence of the period on the rate of protein degradation in minimal models for circadian oscillations,” Philosophical Transactions of the Royal Society A, vol. 367, no. 1908, pp. 4665–4683, 2009.
View at: Publisher Site | Google Scholar
J. C. Leloup and A. Goldbeter, “Toward a detailed computational model for the mammalian circadian clock,” Proceedings of the National Academy of Sciences of the United States of America, vol. 100, no. 12, pp. 7051–7056, 2003.
View at: Publisher Site | Google Scholar
J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, 2nd edition, 2008.
J. C. Butcher and G. Wanner, “Runge-Kutta methods: some historical notes,” Applied Numerical Mathematics, vol. 22, no. 1–3, pp. 113–151, 1996.
View at: Google Scholar
E. Hairer, S. P. Nørsett, and G. Wanner, Ordinary Differential Equations I: Nonstiff Problems, Springer, Berlin, Germany, 1993.
E. Hairer, C. Lubich, and G. Wanner, Solving Geometric Numerical Integration: Structure-Preserving Algorithms, Springer, 2nd edition, 2006.
X. You, Y. Zhang, and J. Zhao, “Trigonometrically-fitted Scheifele two-step methods for perturbed oscillators,” Computer Physics Communications, vol. 182, no. 7, pp. 1481–1490, 2011.
View at: Publisher Site | Google Scholar
X. You, “Limit-cycle-preserving simulation of gene regulatory oscillators,” Discrete Dynamics in Nature and Society, vol. 2012, Article ID 673296, 22 pages, 2012.
View at: Publisher Site | Google Scholar
S. Blanes and P. C. Moan, “Practical symplectic partitioned Runge-Kutta and Runge-Kutta-Nyström methods,” Journal of Computational and Applied Mathematics, vol. 142, no. 2, pp. 313–330, 2002.
View at: Publisher Site | Google Scholar
M. Xiao and J. Cao, “Genetic oscillation deduced from Hopf bifurcation in a genetic regulatory network with delays,” Mathematical Biosciences, vol. 215, no. 1, pp. 55–63, 2008.
View at: Publisher Site | Google Scholar
S. Widder, J. Schicho, and P. Schuster, “Dynamic patterns of gene regulation I: simple two-gene systems,” Journal of Theoretical Biology, vol. 246, no. 3, pp. 395–419, 2007.
View at: Publisher Site | Google Scholar
I. M. M. van Leeuwen, I. Sanders, O. Staples, S. Lain, and A. J. Munro, “Numerical and experimental analysis of the p53-mdm2 regulatory pathway,” in Digital Ecosystems, F. A. Basile Colugnati, L. C. Rodrigues Lopes, and S. F. Almeida Barretto, Eds., vol. 67 of Lecture Notes of the Institute for Computer Sciences, Social Informatics and Telecommunications Engineering, pp. 266–284, 2010.
View at: Google Scholar
G. Strang, “On the construction and comparison of difference schemes,” SIAM Journal on Numerical Analysis, vol. 5, pp. 506–517, 1968.
View at: Google Scholar
R. D. Ruth, “Canonical integration technique,” IEEE Transactions on Nuclear Science, vol. NS-30, no. 4, pp. 2669–2671, 1983.
View at: Google Scholar
M. Suzuki, “Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations,” Physics Letters A, vol. 146, no. 6, pp. 319–323, 1990.
View at: Google Scholar
M. Suzuki, “General theory of higher-order decomposition of exponential operators and symplectic integrators,” Physics Letters A, vol. 165, no. 5-6, pp. 387–395, 1992.
View at: Google Scholar
H. Yoshida, “Construction of higher order symplectic integrators,” Physics Letters A, vol. 150, no. 5–7, pp. 262–268, 1990.
View at: Google Scholar
H. Yang, X. Wu, X. You, and Y. Fang, “Extended RKN-type methods for numerical integration of perturbed oscillators,” Computer Physics Communications, vol. 180, no. 10, pp. 1777–1794, 2009.
View at: Publisher Site | Google Scholar
Z. Chen, X. You, W. Shi, and Z. Liu, “Symmetric and symplectic ERKN methods for oscillatory Hamiltonian systems,” Computer Physics Communications, vol. 183, no. 1, pp. 86–98, 2012.
View at: Publisher Site | Google Scholar
J. Li, B. Wang, X. You, and X. Wu, “Two-step extended RKN methods for oscillatory systems,” Computer Physics Communications, vol. 182, no. 12, pp. 2486–2507, 2011.
View at: Publisher Site | Google Scholar
X. You, Z. Chen, and Y. Fang, “New explicit adapted Numerov methods for second-order oscillatory differential equations,” Applied Mathematics and Computation, vol. 219, pp. 6241–6255, 2013.
View at: Google Scholar

Copyright

Copyright © 2014 Xiong You et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1127

Downloads

1138

Citations