Abstract
A novel family of exponential Runge-Kutta (expRK) methods are designed incorporating the stable steady-state structure of genetic regulatory systems. A natural and convenient approach to constructing new expRK methods on the base of traditional RK methods is provided. In the numerical integration of the one-gene, two-gene, and p53-mdm2 regulatory systems, the new expRK methods are shown to be more accurate than their prototype RK methods. Moreover, for nonstiff genetic regulatory systems, the expRK methods are more efficient than some traditional exponential RK integrators in the scientific literature.
1. Introduction
One of the challenges in systems biology is to understand how biochemical molecules, such as DNAs, mRNAs, and proteins, interact to form harmonic and uniform cellular systems which give rise to life (see [1, 2]). The synthetic genetic regulatory networks (GRNs) play an important role in the investigation of protein regulation processes in living organisms (see [3–6]). By introducing ordinary differential equations (ODEs) to describe the rates of change in the concentrations of biochemical molecules, such as mRNAs and proteins, more detailed understanding and insights of the dynamic behavior exhibited by biological systems can be achieved (see [7]). The first attempt to model the oscillation in genetic regulation in terms of ODEs was made by Goodwin [8]. A standard presentation of general regulatory dynamics can be found in the monographs by Thomas and D’Ari [9] and by Fall et al. [10]. Iwamoto et al. [11] presented a dynamical model of the DNA damage signaling pathway that includes p53 and whole cell cycle regulation and explored the relationship between p53 oscillation and cell fate selection. ODEs models admit mathematically qualitative and quantitative analysis to reveal the profound properties from steady states with stability, bistability, oscillation, and limit cycles to chaos (see [12–17] and the references therein).
A typical system of ODEs governing an -gene activation-inhibition system has the form (see Polynikis et al. [15]) where, for , is the concentration of mRNA produced by gene , is the concentration of protein translated from mRNA , is the degradation rate of , and is the degradation rate of . Function is the translation function. Function is the regulation function, typically taking the form of a sum of products of functions . If protein has no effect on gene , does not appear in . The partial derivative if protein is an activator of gene and if protein is an inhibitor of gene . The genetic regulatory system (1) can be written in matrix formwhere and are -dimensional vectors representing the concentrations of mRNAs and proteins at time , respectively, and , , , and are diagonal matrices.
The analytical solution of system (1) is in general not acquirable due to the nonlinearity of the functions and . Therefore, in order to explore the dynamics of the gene regulatory system (1), one usually resorts to numerical solution. For example, Shinto et al. [18] proposed a kinetic simulation model of metabolic pathways that describes the dynamic behaviors of metabolites in acetone-butanol-ethanol (ABE) production by Clostridium saccharoperbutylacetonicum N1-4 using a simulator WinBEST-KIT. So far, differential equations for genetic regulation are mostly solved by the classical four-stage Runge-Kutta (RK) method or by the Runge-Kutta-Fehlberg adaptive method (see Hairer et al. [19]). However, general-purpose RK methods have not taken into account the special structure of system (1) and fail to capture the dynamical features of the system effectively, especially in the long time simulation (see Figure 11). Thanks to new advances in the last two decades, new approaches have been developed aiming at preserving the intrinsic geometric or physical structures of the true solution. A comprehensive account of structure-preserving algorithms can be found in the monographs by Stuart and Humphries [20], Hairer et al. [21], and Wu et al. [22]. Recently Hochbruck and Ostermann [23] investigated exponential Runge-Kutta methods for initial value problems of parabolic differential equations. This type of methods simulates exactly the linear structure of the differential equations. Defterli et al. [24] and Weber et al. [25] considered discretizing and optimizing the so-called gene-environment networks based on usually finite data series.
From the dynamics point of view, there are two basic categories of genetic regulatory systems: Category 1 consists of systems having sustained oscillation, such as limit cycles; Category 2 consists of systems having steady states. For the genetic regulatory system (1) with a limit-cycle structure, You [26] proposed a new class of phase-fitted and amplification-fitted Runge-Kutta type methods which were shown to be more effective and more efficient than the traditional Runge-Kutta methods of the same order. Very recently, You et al. [27] developed a splitting approach for genetic regulatory systems with a stable steady state. In the numerical simulation, the new splitting methods constructed in that paper are shown to be remarkably more effective and more suitable for long-term computation with large steps than the general-purpose Runge-Kutta methods. In order to respect the oscillatory feature of the solution of some genetic regulatory systems, Chen et al. [28] developed a new type of exponentially fitted TDRK (EFTDRK) methods. Zhang et al. [29] constructed a family of phase-fitted symmetric splitting methods of order two and order four. The result of the numerical experiment on the Lotka–Volterra system shows that the new phase-fitted symmetric splitting methods are more effective than their prototype splitting methods and can preserve the invariant of the system in the long term compared with the classical Runge-Kutta method of order four.
The purpose of this paper is to develop a novel type of exponential RK methods for the simulation of genetic regulatory systems which have an asymptotically stable steady state. In Section 2 we present the general scheme of exponential Runge-Kutta (expRK) methods for solving initial value problems of ODEs based on a matrix form of the variation-of-constants formula. A convenient approach of transiting traditional RK methods into a special type of expRK methods is given. In Section 3 we integrate the above three regulatory systems by the new expRK methods as well as their prototype RK methods for comparison. Section 4 is devoted to conclusive remarks. In Appendix, we analyze the linear stability and phase properties of the expRK methods.
2. Exponential Runge-Kutta Methods
2.1. Formulation of Exponential RK Methods for Systems with a Stable Steady-State Structure
Prior to dealing with the genetic regulatory system (2) numerically, we first consider the general initial value problem (IVP) of the autonomous system of ODEs where and “” represents the derivative of with respect to time. We make the following assumptions on system (3):(A1)The origin is a steady state of the system; that is, .(A2)The Jacobian has eigenvalues of negative real parts in a neighborhood of the origin. Then, according to the theorem in Section 8.5 of Hirsch et al. [30], the origin is an asymptotically stable steady state; that is, for every solution of system (3) through a point in the neighborhood of the steady state, .(A3)The function is continuously differentiable and satisfies the Lipschitz condition; that is, there exists a constant (called the Lipschitz constant) such that for all .
Let us recall the general-purpose Runge-Kutta methods for IVP (3).
Definition 1. An -stage Runge-Kutta (RK) method for system (3) has the schemewhere , are real numbers.
Scheme (5) can be expressed briefly by the Butcher tableau of its coefficientsThe order conditions for the RK method (5) can be found in Hairer et al. [19]. Note that the general RK scheme (5) does not take into account the special structure of the equilibrium structure of the system so that the computational results are usually not satisfactory. In order to simulate more effectively system (3) with a steady state at the origin, we rewrite problem (3) in an equivalent formwhere the matrix , the function , is the Euclidean norm, and is called the rate matrix of system (3). The matrix form of the variation-of-constants formula for system (3) is given bywhere and are real numbers and . Approximating the integral on the right-hand side of (8) by some effective quadrature formulas leads to a numerical integrator. In general we have the following definition of the so-called exponential Runge-Kutta methods.
Definition 2. An -stage exponential Runge-Kutta method for system (3) has the schemewhere , are real numbers and and , are real matrix-valued functions of matrix .
It is convenient to express scheme (9) by the Butcher tableauwhere . For a comprehensive review of exponential integrators with the construction, analysis of convergence and error bounds, order conditions, example integrators, and applications, the reader is referred to Hochbruck and Ostermann [23].
It is noted that if we need to integrate a system near a stable steady state , the exponential RK scheme (9) should take the formwhere with .
2.2. A Special Class of Exponential RK Methods
Based on an RK method (5), we can formulate, as a special case of the exponential RK method (9), the following scheme for system (3): where , , and , are the constant coefficients of the RK method (5). Note that as (), scheme (12) reduces to the RK method (5). The latter is called the prototype RK method of the former. In the sequel of this paper, scheme (12) will be referred to as an expRK method. Its Butcher tableau can be written as follows:where , , and .
In Kronecker’s notation, scheme (12) can be written as where is the unit matrix and .
Among simple examples are the following:(a)The exponential Euler method (explicit), denoted by expEuler:(b)The exponential backward Euler method (implicit), denoted by expImEuler:(c)The exponential trapezoidal rule (implicit), denoted by expTrap:(d)The exponential Heun rule (explicit), denoted by expHeun:(e)The exponential midpoint rule (implicit), denoted by expMid:
Two typical expRK methods, denoted by expRK3/8 and expRK4, have the prototype RK methods of order four, denoted by RK3/8 and RK4, respectively, whose respective Butcher tableaux are given in Page 138 of [19]
The (algebraic) order is a measure of the accuracy of numerical method. A method is said to have order if its local error .
Theorem 3. The expRK method (12) has the same algebraic order as its prototype RK method.
Proof. The conclusion follows by expanding the exponential functions , , and in scheme (12) in series of , hence in series of , and comparing this series with that of the true solution.
The next theorem asserts that exponential RK method (12) preserves the steady state of system (3).
Theorem 4. Suppose that in system (3) satisfies the Lipschitz condition and the origin is a steady state of the system; that is, . Then the origin is also a fixed point of the expRK method (12) for small step size .
The proof is given in Appendix A.
In order to apply the expRK method (9) or (12) to system (2), we first use a coordinate transform to translate the steady state of the system to the origin and yieldswhere is the Jacobian matrix of at point and . Then system (2) can be written in the form (7) with where the rate matrix and the function .
3. Numerical Illustrations
From Theorem 4, contrast to traditional RK methods, exponential RK methods, especially expRK methods, retain the rate of growth, phase, and amplification of the exact solution of the test equation (B.1) without error. Then it is reasonable to expect expRK (12) to be more effective than their prototype RK methods. In this section, in order to compare their effectiveness, we apply them to three test systems—one-gene self-regulation system, two-gene cross-regulation system, and the p53-mdm2 system.
(I) A One-Gene System of Self-Regulation. The first model we consider is the one-gene system with self-regulation given bywhere represents the action of an inhibitory protein that acts as a dimer and are positive constants. For a similar model with delays see Xiao and Cao [31].
(II) A Two-Gene System with Cross-Regulation. The second model is a two-gene activation-inhibition system with cross-regulation (studied by Polynikis et al. [15], Widder et al. [17], Chen et al. [32], You [26], and You et al. [27])where, for , is the concentration of mRNA produced by gene , is the concentration of protein , is the maximal transcription rate of gene , is the translation rate of mRNA , and and are the degradation rates of mRNA and protein , respectively. The functionsare the Hill functions of activation and repression, respectively. The parameters are the expression thresholds. The integer value of , called the Hill coefficient, stands for the number of protein monomers required for saturation of binding to DNA. It is easy to see that the activation function is increasing in and the repression function is decreasing in .
(III) The p53-mdm2 System. The third model is for the damped oscillation of the p53-mdm2 regulatory pathway. Strictly speaking, this system is not of the form (1). We adopt this model since its solutions also have a stable steady-state structure of interest. The system, given by van Leeuwen et al. [33] with the small transient stress stimulus , has the formwhere represents the concentration of the p53 tumor 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, 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.1. Accuracy Test
3.1.1. The One-Gene System
Steady states of system (24) can be determined by the cubic equation and the relation . If the system is written in the form (7), the rate matrixand the function where . With the parameter values (provided by [31])this system has a unique positive steady state where the rate matrix has eigenvalues , where is the imaginary unit satisfying . Since the two eigenvalues both have negative real parts, the steady state is asymptotically stable. Figure 1(a) shows three solution trajectories on the phase plane starting at , respectively. Figure 1(b) shows the time evolution of concentrations of mRNA and protein starting at .
(a)
(b)
With the above values of parameters and initial data, we integrate system (25) on the time interval by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods. We plot the error growth of the protein on the time interval . The numerical results are presented in Figures 2 and 3.
(a)
(b)
(a)
(b)
The system is solved on the time interval with initial values of mRNA and protein and and with different step sizes. The numerical results are presented in Tables 1 and 2.
3.1.2. The Two-Gene System
The steady states of system (25) are determined by the equations
Putting in the form (7), system (25) has the rate matrixand the functionwith . The characteristic equation of the rate matrix iswhereFor a certain value of , the real part of one of the eigenvalues crosses zero, indicating a loss of stability through a Hopf bifurcation. For and , [17] has calculated this value explicitly as
In our experiment, we take the values of parameters as follows:Figure 4(a) shows three solution trajectories projected on the mRNA 1-protein 1 plane starting at , , and , respectively. Figure 4(b) shows the time evolution of concentrations of protein 1 and protein 2.
(a)
(b)
We solve (31) for by Newton’s iteration and then substitute it into (32) obtainingThe rate matrix (33) has the eigenvalues with negative real parts: Therefore the steady state is asymptotically stable.
With the values of parameters (38) and initial data and initial values , , , and , we integrate system (25) on the time interval by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step sizes , respectively. In Figures 5 and 6, we plot the global error growth of protein 1 on the time interval .
(a)
(b)
(a)
(b)
In Tables 3 and 4, average errors are compared for differential step sizes.
3.1.3. The p53-mdm2 System
The steady state is determined by the following equations:Putting in the form (7), system (27) has the rate matrixand the functionwhere , andWe use the parameter values (see [33]) as follows:For simplicity, we take the small function . The system has a unique steady state . The rate matrix has the eigenvalues Since all the eigenvalues have negative real parts, the steady state is asymptotically stable.
Figure 7(a) shows three solution trajectories projected on the inactive p53-complex plane starting at , respectively. Figure 7(b) shows the time evolution of concentrations of p53 and mdm2.
(a)
(b)
With the values of parameters in (45) and initial data and initial values , , , and , we integrate system (25) on the time interval by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step size . In Figures 8 and 9, we plot the global error growth of the inactive p53 on the time interval .
(a)
(b)
(a)
(b)
The problem is solved on the interval for differential step sizes and the average errors are presented in Tables 5 and 6.
From Figures 2, 3, 5, 6, 8, and 9 and Tables 1–6, we can see that the new expRK methods expEuler, expHeun, expRK3/8, and expRK4 are much more accurate than their corresponding prototype methods.
3.2. Efficiency Test
In this subsection we will compare the simulation efficiency of the newly constructed exponential RK methods with some famous exponential integrators. The integrators we choose for comparison are listed as follows:(i)expRK3/8: the expRK method defined by (12) whose prototype RK method is given by (20)(ii)expRK4: the expRK method defined by (12) whose prototype RK method is given by (21)(iii)COX-MATTHEWS: the exponential RK method given by Cox and Matthews [34](iv)KROGSTAD: the exponential RK method given by Krogstad [35](v)STREHMEL-WEINER: the exponential RK method given by Strehmel [36] (Example 4.5.5)(vi)HOCHBRUCK-OSTERMANN: the exponential RK method given by Hochbruck and Ostermann [23]
The criterion for the efficiency is the digital logarithm of the global error against the CPU-time consumed. System (25) with parameters (38) is solved on the time interval with the step sizes . The numerical results are displayed in Figure 10, where we can see that the new exponential RK methods expEuler, expHeun, expRK3/8, and expRK4 are much more efficient than the other exponential RK methods we select from the literature. For the two-gene system, among all the exponent integrators we consider, expEuler, though the simplest, turns out to be the most efficient. It is also interesting to observe that as integration step size decreases from , expRK4 cannot produce smaller error, just increasing the computation effort.
4. Conclusions and Discussions
Most genetic regulatory systems carry their own structures, such as stable steady states and sustained oscillation, bistability. Traditional Runge-Kutta (RK) methods have not taken into account these characteristic structures and may give misleading information. To see this, one only needs to integrate the two-gene system (25) by RK4 with step size . The simulation result, as presented in Figure 11, is qualitatively wrong. The exponential RK methods for system (3) originate from the discretization of the matrix form of variation-of-constants formula (8). This type of integrators has the property that they can integrate exactly linear systems of ODEs and thus can preserve the steady-state structure of system (3). From the numerical results presented in Section 3, an expRK method is more accurate than its traditional prototype RK method for long-term simulation with large step sizes. On the other hand, despite the simple form, the new expRK methods considered in this paper are tested to be more efficient than those prominent exponential RK methods when they are applied to genetic regulatory systems. The following advantages contribute to the high accuracy and high efficiency of the expRK methods compared to the classical RK methods.
(a) The scheme of the expRK methods recovers by the exponential functions the principal oscillatory structure of the true solution, which is contained in linear part (the Jacobian) of the system.
(b) The construction of the scheme is very simple and the coefficients are immediately obtained from a classical Runge-Kutta (RK) method.
(c) As shown in Appendix, the expRK methods (RK3/8 and RK4) have distortion and dissipation of the same order as their prototype RK methods but have dispersion of one order higher than their prototype RK methods. As the principal frequency is estimated accurately enough, that is, the ratio of the error of estimation and the testing frequency , the coefficients of the leading terms of distortion and dissipation of the new expRK methods are much less than those of their prototype RK methods. Moreover, if is known to be the exact frequency (), expRK methods become zero-distortive, zero-dispersive, and zero-dissipative.
Although there may exist other methods of higher algebraic order that have higher accuracy, the expRK methods (12) are the most natural ones and are the most convenient to use. Note also that, similar to the approach of ode45, we can control the step by embedding two expRK methods of orders and into a pair to achieve higher efficiency.
Indeed, the test problem (25) is assumed to be nonstiff and this may be why our expRK methods outperform the existing exponential Runge-Kutta methods of Cox and Matthews, Hochbruck and Ostermann, Krogstad, and Stehmel and Weiner. The consideration of expRK methods for stiff genetic regulatory system (such as the p53-mdm2 systems, whose Jacobian possesses eigenvalues with large negative real parts or with purely imaginary eigenvalues of large modulus) is an interesting theme for future work. These systems contain different time scales due to coexistence of fast reactions and slow reactions which are frequently encountered in real cell processes. The expRK methods adapted to stiff systems will have a more delicate scheme to incorporate the stiff structure.
Appendix
A. Proof of Theorem 4
A fixed point of scheme (12) is a point such that implies . We have assumed that . Suppose that in scheme (12). The internal stages of scheme (12) define a nonlinear operator on . For , the internal stages of system (12) arewhere if . Therefore, is a contraction on . By the Banach contraction theorem, the operator has a unique fixed point ; that is, for which implies that . We have proven that the origin is the unique fixed point of the exponential RK method (12).
B. Analysis of Linear Stability, Distortion, Dispersion, and Dissipation
In this section, in order to examine the behavior of an RK or expRK method for the dynamical system defined by ODE (3) with a stable steady state , we analyze its linear stability and estimate the orders of distortion, dispersion, and dissipation.
Let us consider the linear scalar test equationwhere the complex number is an estimate of the principal rate and is the error of the estimation. Applying an RK method (5) or an expRK method (12) to (B.1) yieldswhere is called the stability function of the method.
For an -stage RK method,where is the identity matrix and the -dimensional vector . For an explicit -stage RK method, since the matrix is an lower-triangular, then and
On the other hand, an -stage expRK method (12) has a distortionIf the expRK method is explicit,
Definition B.1. For an -stage method with the stability function , the region in the planeis called the stability region of the method.
Figure 12 presents the stability regions of RK4 (left) and expRK4, respectively.
(a)
(b)
Definition B.2. For an -stage expRK method (5) with the stability function , the quantityis called the distortion (or rate error) of the method. If as , then the method is called distortive of order . If , then the method is called zero-distortive.
The distortion of RK4 method is On the other hand, if we denote the ratio , then and and the distortion of the expRK4 method isIt can be seen from (B.9) and (B.10) that RK4 and expRK4 are both distortive of order four. However, if an estimate of the true rate is accurate enough, that is, , then the coefficient of the leading term of the distortion for the expRK4 method is much smaller than that for the RK method. Moreover, if the frequency can be estimated exactly (), then all the expRK methods are zero-distortive.
It has been observed that most genetic regulatory systems have oscillatory solutions. Therefore it is of importance to measure to what extent a numerical method can preserve the oscillation. In order to compare the accuracy RK methods and that of expRK methods in preserving the oscillatory properties of the exact solution, we assume that and in the test equation (B.1) are purely imaginary; that is, .
Definition B.3. For a method with stability function , denote the real and imaginary parts of the stability function by and , respectively, where and . The quantitiesare called the dispersion (or phase lag) and dissipation (or amplification factor error) of the method, respectively. If and as , then the method is called dispersive of order and dissipative of order , respectively. If and , then the method is called zero-dispersive and zero-dissipative, respectively.
For both RK4 and RK3/8, the dispersion and dissipation satisfyTherefore, RK4 and RK3/8 both are dispersive of order four and dissipative of order five.
Let . Then the dispersion and dissipation for both expRK4 and expRK3/8 areTherefore, expRK4 and expRK3/8 both are dispersive of order five and dissipative of order five. Then expRK4 (expRK3/8) has dispersion of one order higher than RK4 (RK3/8). If , then the coefficient of the leading term of the dissipation for the expRK4 (expRK3/8) method is much smaller than that for RK4 (RK3/8). Furthermore, if , then all the expRK methods are zero-dispersive and zero-dissipative.
The above analysis explains why, when the rate of the problem is estimated accurately enough, an expRK method is more accurate than its prototype RK method in preserving the distortion, dispersion, and dissipation of the exact solution.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was partially supported by the National Natural Science Foundation of China (NSFC) (Grant no. 11171155), the State Key Program of Natural Science Foundation of China (Grant no. 31330067), and the Fundamental Research Fund for the Central Universities (Grant no. KYZ201424).