Steady-State-Preserving Simulation of Genetic Regulatory Systems

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.


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][4][5][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][13][14][15][16][17] and the references therein).
A typical system of ODEs governing an -gene activation-inhibition system has the form (see Polynikis et al. [15]) Transcription:̇= ( ) − , where, for = 1, . . . , , 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 1 ( 1 ), . . . , ( ). If protein has no effect on gene , ( ) does not appear in . The partial derivative / > 0 if protein is an activator of gene and / < 0 if protein is an inhibitor of gene . The genetic regulatory system (1) can be written in matrix forṁ ( ) = −Γ ( ) + ( ( )) , 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 acetonebutanol-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 amplificationfitted 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 longterm 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.

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 ODEṡ= where : [0, +∞) → R and "" represents the derivative of with respect to time. We make the following assumptions on system (3): (A1) The origin * = 0 is a steady state of the system; that is, (0) = 0. (A2) The Jacobian ( / )(0) 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, lim →+∞ ( ) = 0.
(A3) The function : R → R is continuously differentiable and satisfies the Lipschitz condition; that is, there exists a constant (called the Lipschitz constant) such that for all , ∈ R .
Let us recall the general-purpose Runge-Kutta methods for IVP (3).
Computational and Mathematical Methods in Medicine 3 Scheme (5) can be expressed briefly by the Butcher tableau of its coefficients The 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 forṁ − Ω = ( ) , where the matrix Ω = ( / )(0), the function ( ) = ( )− Ω = O(‖ ‖ 2 ), ‖ ⋅ ‖ is the Euclidean norm, and Ω is called the rate matrix of system (3). The matrix form of the variationof-constants formula for system (3) is given by where and ℎ are real numbers and = ℎ, = 0, 1, . . .. 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.
It is convenient to express scheme (9) by the Butcher tableau where exp( ) = (exp( 1 ) , . . . , exp( ) ) . 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 * ̸ = 0, the exponential RK scheme (9) should take the form where = ℎΩ with Ω = ( / )( * ).

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 , , = 1, . . . , , are the constant coefficients of the RK method (5). Note that as Ω → 0 ( → 0), 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 = (1, . . . , 1) , exp(( − ) = (exp((1 − 1 ) ) , . . . , exp((1 − ) ) ) , and exp( ) exp(− ) = (exp(( − ) ) ) × . In Kronecker's notation, scheme (12) can be written as where is the × unit matrix and exp Among simple examples are the following: (a) The exponential Euler method (explicit), denoted by expEuler: 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 LE = O(ℎ +1 ). Theorem 3. The expRK method (12) has the same algebraic order as its prototype RK method.

Computational and Mathematical Methods in Medicine 5
Proof. The conclusion follows by expanding the exponential functions exp( ), exp(( − ) ), and exp((1 − ) ) 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). (3) satisfies the Lipschitz condition and the origin is a steady state of the system; that is, (0) = 0. Then the origin is also a fixed point of the expRK method (12) for small step size ℎ.

Theorem 4. Suppose that in system
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 ( ) = ( )− * , V( ) = ( ) − * to translate the steady state ( * , * ) of the system to the origin and yieldṡ where ( * ) is the Jacobian matrix of ( ) at point * and Then system (2) can be written in the form (7) and the function ( ) = ( (V) , 0) .

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 bẏ where ( ( )) = /(1 + ( ) 2 / 2 ) 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 = 1,2, 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 functions are the Hill functions of activation and repression, respectively. The parameters 1 , 2 are the expression thresholds. The integer value of ( = 1, 2), 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 2 and the repression function − is decreasing in 1 .
(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 ( ) = 0, has the forṁ= where 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, * ( * = , 0, 1) 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.  (24) can be determined by the cubic equation (1/ 2 ) * 3 + * + / = 0 and the relation * = ( / ) * . If the system is written in the form (7), the rate matrix and the function this system has a unique positive steady state ( * , * ) = (0.6, 2) where the rate matrix Ω has eigenvalues 1,2 = −1.2500 ± 2.9767 , where is the imaginary unit satisfying 2 = −1. Since the two eigenvalues both have negative real parts, the steady state is asymptotically stable. Figure 1 Putting in the form (7), system (25) has the rate matrix and the function ) . The characteristic equation of the rate matrix Ω is where = 1 2 1 2 For a certain value of = Hopf , the real part of one of the eigenvalues crosses zero, indicating a loss of stability through a Hopf bifurcation. For 1 > 1 and 2 > 1, [17] has calculated this value explicitly as In our experiment, we take the values of parameters as follows: The rate matrix (33) has the eigenvalues with negative real parts: Therefore the steady state * is asymptotically stable.

Computational and Mathematical Methods in Medicine
In Tables 3 and 4, average errors are compared for differential step sizes.
(42) and the function We use the parameter values (see [33]) as follows: For simplicity, we take the small function ( ) ≡ 0. The system has a unique steady state ( * , * , * , * ) = (9.42094, 0.0372868, 3.49529, 0). The rate matrix Ω has the eigenvalues 1 = −38.4766, 2,3 = −0.0028 ± 0.0220 , Since all the eigenvalues have negative real parts, the steady state is asymptotically stable. The problem is solved on the interval [0, 50] for differential step sizes and the average errors are presented in Tables 5  and 6.

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:  (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 [0, 100] with the step sizes ℎ = 1/2 , = 1, 2, 3, 4. 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 ℎ = 1/2, expRK4 cannot produce smaller error, just increasing the computation effort.

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 ℎ = 1.32. 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 Table 2: One-gene system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 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  Step size RK3/8 expRK3/8 RK4 expRK4   the error of estimation and the testing frequency | | ≪ 1, 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 ( = 1), expRK methods become zerodistortive, 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 + 1 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.

B. Analysis of Linear Stability, Distortion, Dispersion, and Dissipation
In this section, in order to examine the behavior of an RK or expRK method Φ ℎ : → +1 for the dynamical system defined by ODE (3) with a stable steady state * = 0, we analyze its linear stability and estimate the orders of distortion, dispersion, and dissipation.
Let us consider the linear scalar test equatioṅ where 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)  where ( , ) is called the stability function of the method.
For an -stage RK method, where is the × identity matrix and the -dimensional vector = (1, . . . , 1) . For an explicit -stage RK method, since the matrix is an × lower-triangular, then = 0 and On the other hand, an -stage expRK method (12)  is called the stability region of the method. It 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, | | ≪ 1, then the coefficient of the leading term of the distortion Dist( ) for the expRK4 method is much smaller than that for the RK method. Moreover, if the frequency can be estimated exactly ( = 1), 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, Ω = , = , 2 = −1. Therefore, 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 | | ≪ 1, then the coefficient of the leading term of the dissipation Diss( ) for the expRK4 (expRK3/8) method is much smaller than that for RK4 (RK3/8). Furthermore, if = 1, 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.