# Steady-State-Preserving Simulation of Genetic Regulatory Systems

Authors: Ruqiang Zhang; Julius Osato Ehigie; Xilin Hou; Xiong You; Chunlu Yuan
Journal: Computational and Mathematical Methods in Medicine
 (2017)

---

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


---

## Body

## 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 anN-gene activation-inhibition system has the form (see Polynikis et al. [15]) (1)Transcription:r˙i=fiRp-γiri,Translation:p˙i=fiPri-μipi,where, for i=1,…,N, ri is the concentration of mRNA Ri produced by gene gi, pi is the concentration of protein Pi translated from mRNA Ri, γi is the degradation rate of Ri, and μi is the degradation rate of Pi. Function fiP(ri) is the translation function. Function fiR(p) is the regulation function, typically taking the form of a sum of products of functions fi1R(p1),…,fiNR(pN). If protein Pj has no effect on gene gi, fjR(pj) does not appear in fiR. The partial derivative ∂fiR/∂pj>0 if protein Pj is anactivator of gene gi and ∂fiR/∂pj<0 if protein Pj is aninhibitor of gene gi. The genetic regulatory system (1) can be written in matrix form(2)r˙t=-Γrt+Fpt,p˙t=Krt-Mpt,where r(t)=(r1(t),…,rN(t))T and p(t)=(p1(t),p2(t),…,pN(t))T are N-dimensional vectors representing the concentrations of mRNAs and proteins at time t, respectively, and F(p(t))=(F1(p(t)),…,FN(p(t)))T, Γ=diag⁡γ1,…,γN, M=diag⁡μ1,…,μN, and K=diag⁡κ1,…,κN are diagonal matrices.The analytical solution of system (1) is in general not acquirable due to the nonlinearity of the functions fiP(ri) and fijR(pj). 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 byClostridium 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 Section2 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 (3)y˙=fy,t>0,y0=y0,where y:[0,+∞)→Rd and “y˙” represents the derivative of y with respect to time. We make the following assumptions on system (3):(A1)
The originy∗=0 is a steady state of the system; that is, f(0)=0.(A2)
The Jacobian∂f/∂y(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 y(t) of system (3) through a point in the neighborhood of the steady state, limt→+∞y(t)=0.(A3)
The functionf:Rd→Rd is continuously differentiable and satisfies the Lipschitz condition; that is, there exists a constant L (called theLipschitz constant) such that (4)fy-fz≤Ly-zfor ally,z∈Rd.Let us recall the general-purpose Runge-Kutta methods for IVP (3).Definition 1.
Ans-stage Runge-Kutta (RK) method for system (3) has the scheme(5)Yi=yn+h∑j=1saijfYj,i=1,…,s,yn+1=yn+h∑i=1sbifYi,where ci,aij,bi,i,j=1,…,s, are real numbers.Scheme (5) can be expressed briefly by the Butcher tableau of its coefficients(6)
​
c
A
b
T
=
c
1
a
11
⋯
a
1
s
⋮
⋮
⋱
⋮
c
s
a
s
1
⋯
a
s
s
b
1
⋯
b
s
.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 form(7)y˙-Ωy=gy,y0=y0,where the matrix Ω=∂f/∂y(0), the function g(y)=f(y)-Ωy=Oy2, · is the Euclidean norm, and Ω is called therate matrix of system (3). The matrix form of the variation-of-constants formula for system (3) is given by(8)ytn+μh=exp⁡μhΩytn+exp⁡μhΩ∫tntn+μhexp⁡-ξ-tnΩgyξdξ,where μ and h are real numbers and tn=nh,n=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.Definition 2.
Ans-stage exponential Runge-Kutta method for system (3) has the scheme(9)Yi=exp⁡ciVyn+∑j=1saijVhfYj-VYj,i=1,…,s,yn+1=exp⁡Vyn+∑i=1sbiVhfYi-VYi,where ci,i=1,…,s, are real numbers and aij(V) and bi(V),i,j=1,…,s, are real d×d matrix-valued functions of matrix V=hΩ.It is convenient to express scheme (9) by the Butcher tableau(10)
​
c
exp
⁡
c
V
A
V
exp
⁡
V
b
T
V
=
c
1
e
x
p
c
1
V
a
11
V
⋯
a
11
V
⋮
⋮
⋮
⋱
⋮
c
s
e
x
p
c
1
V
a
s
1
V
⋯
a
s
s
V
exp
⁡
V
b
1
V
⋯
b
s
V
,where exp⁡cV=(exp⁡c1VT,…,exp⁡(csV)T)T. 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 systemz′=ψ(z) near a stable steady state z∗≠0, the exponential RK scheme (9) should take the form(11)yn=zn-z∗,Yi=exp⁡ciVyn+∑j=1saijVhψYj+z∗-VYj,i=1,…,s,yn+1=exp⁡Vyn+∑i=1sbiVhψYi+z∗-VYi,zn+1=yn+1+z∗,where V=hΩ with Ω=∂ψ/∂z(z∗).


### 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): (12)Yi=exp⁡ciVyn+∑j=1saijexp⁡ci-cjVhfYj-VYj,i=1,…,s,yn+1=exp⁡Vyn+∑i=1sbiexp⁡1-ciVhfYi-VYi,where ci, aij, and bi,i,j=1,…,s, are the constant coefficients of the RK method (5). Note that as Ω→0 (V→0), scheme (12) reduces to the RK method (5). The latter is called theprototype RK method of the former. In the sequel of this paper, scheme (12) will be referred to as anexpRK method. Its Butcher tableau can be written as follows:(13)
​
c
exp
⁡
c
V
e
x
p
c
V
T
A
exp
⁡
-
c
V
exp
⁡
V
b
T
exp
⁡
e
-
c
V
=
c
1
exp
⁡
c
1
V
a
11
a
12
exp
⁡
c
1
-
c
2
V
⋯
a
1
s
exp
⁡
c
1
-
c
s
V
⋮
⋮
⋮
⋮
⋱
⋮
c
s
exp
⁡
c
1
V
a
s
1
exp
⁡
c
s
-
c
1
V
a
s
2
exp
⁡
c
s
-
c
2
V
⋯
a
s
s
exp
⁡
V
b
1
exp
⁡
1
-
c
1
V
b
2
exp
⁡
1
-
c
2
V
⋯
b
s
exp
⁡
1
-
c
s
V
,where e=(1,…,1)T, exp⁡((e-c)V=(exp⁡((1-c1)V)T,…,exp⁡((1-cs)V)T)T, and exp⁡(cV)TAexp⁡(-cV)=(exp⁡((ci-cj)V)aij)s×s.In Kronecker’s notation, scheme (12) can be written as (14)Y=exp⁡cVyn⊗e+exp⁡cVA⊗Idexp⁡-cVhfY-VY,yn+1=exp⁡Vyn+exp⁡VbT⊗Idexp⁡-cVhfY-VY,where Id is the d×d unit matrix and exp⁡-cV(hf(Y)-VY)=(exp⁡((-c1)V)T(hf(Y1)-VY1),…,exp⁡((-cs)V)T(hf(Ys)-VYs))T.Among simple examples are the following:(a)
Theexponential Euler method (explicit), denoted by expEuler:(15)yn+1=exp⁡Vyn+exp⁡Vhfyn-Vyn(b)
Theexponential backward Euler method (implicit), denoted by expImEuler:(16)yn+1=exp⁡Vyn+hfyn+1-Vyn+1(c)
Theexponential trapezoidal rule (implicit), denoted by expTrap:(17)yn+1=expVyn+expV2hfyn-Vyn+12hfyn+1-Vyn+1(d)
Theexponential Heun rule (explicit), denoted by expHeun:(18)yp=exp⁡Vyn+exp⁡Vhfyn-Vyn,yc=exp⁡Vyn+hfyp-Vyp,yn+1=yp+yc2(e)
Theexponential midpoint rule (implicit), denoted by expMid:(19)yn+1=exp⁡Vyn+exp⁡V2hfexp⁡V/2yn+exp⁡-V/2yn+12-12VexpV2yn+exp-V2yn+1Two 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](20)
​
0
1
3
1
3
2
3
-
1
3
1
1
1
-
1
1
1
8
3
8
3
8
1
8
,
(21)
​
0
1
2
1
2
1
2
0
1
2
1
0
0
1
1
6
2
6
2
6
1
6The (algebraic) order is a measure of the accuracy of numerical method. A method is said to haveorderp if its local error LE=O(hp+1).Theorem 3.
The expRK method (12) has the same algebraic order as its prototype RK method.Proof.
The conclusion follows by expanding the exponential functionsexp⁡(ciV), exp⁡((ci-cj)V), and exp⁡((1-ci)V) in scheme (12) in series of V=hΩ, hence in series of h, 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 thatf in system (3) satisfies the Lipschitz condition and the origin is a steady state of the system; that is, f(0)=0. Then the origin is also a fixed point of the expRK method (12) for small step size h.The proof is given in AppendixA.In order to apply the expRK method (9) or (12) to system (2), we first use a coordinate transform u(t)=r(t)-r∗,v(t)=p(t)-p∗ to translate the steady state (r∗,p∗) of the system to the origin and yields(22)u˙t=-Γut+F′p∗vt+Gvt,v˙t=Kut-Mvt,where F′(p∗) is the Jacobian matrix of F(p) at point p∗ and G(v(t))=F(p∗+v(t))-F′(p∗)v(t)-F(p∗). Then system (2) can be written in the form (7) with y(t)=(u(t),v(t))T where the rate matrix (23)Ω=-ΓF′p∗K-Mand the function g(y)=(G(v)T,0)T.


## 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 (3)y˙=fy,t>0,y0=y0,where y:[0,+∞)→Rd and “y˙” represents the derivative of y with respect to time. We make the following assumptions on system (3):(A1)
The originy∗=0 is a steady state of the system; that is, f(0)=0.(A2)
The Jacobian∂f/∂y(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 y(t) of system (3) through a point in the neighborhood of the steady state, limt→+∞y(t)=0.(A3)
The functionf:Rd→Rd is continuously differentiable and satisfies the Lipschitz condition; that is, there exists a constant L (called theLipschitz constant) such that (4)fy-fz≤Ly-zfor ally,z∈Rd.Let us recall the general-purpose Runge-Kutta methods for IVP (3).Definition 1.
Ans-stage Runge-Kutta (RK) method for system (3) has the scheme(5)Yi=yn+h∑j=1saijfYj,i=1,…,s,yn+1=yn+h∑i=1sbifYi,where ci,aij,bi,i,j=1,…,s, are real numbers.Scheme (5) can be expressed briefly by the Butcher tableau of its coefficients(6)
​
c
A
b
T
=
c
1
a
11
⋯
a
1
s
⋮
⋮
⋱
⋮
c
s
a
s
1
⋯
a
s
s
b
1
⋯
b
s
.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 form(7)y˙-Ωy=gy,y0=y0,where the matrix Ω=∂f/∂y(0), the function g(y)=f(y)-Ωy=Oy2, · is the Euclidean norm, and Ω is called therate matrix of system (3). The matrix form of the variation-of-constants formula for system (3) is given by(8)ytn+μh=exp⁡μhΩytn+exp⁡μhΩ∫tntn+μhexp⁡-ξ-tnΩgyξdξ,where μ and h are real numbers and tn=nh,n=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.Definition 2.
Ans-stage exponential Runge-Kutta method for system (3) has the scheme(9)Yi=exp⁡ciVyn+∑j=1saijVhfYj-VYj,i=1,…,s,yn+1=exp⁡Vyn+∑i=1sbiVhfYi-VYi,where ci,i=1,…,s, are real numbers and aij(V) and bi(V),i,j=1,…,s, are real d×d matrix-valued functions of matrix V=hΩ.It is convenient to express scheme (9) by the Butcher tableau(10)
​
c
exp
⁡
c
V
A
V
exp
⁡
V
b
T
V
=
c
1
e
x
p
c
1
V
a
11
V
⋯
a
11
V
⋮
⋮
⋮
⋱
⋮
c
s
e
x
p
c
1
V
a
s
1
V
⋯
a
s
s
V
exp
⁡
V
b
1
V
⋯
b
s
V
,where exp⁡cV=(exp⁡c1VT,…,exp⁡(csV)T)T. 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 systemz′=ψ(z) near a stable steady state z∗≠0, the exponential RK scheme (9) should take the form(11)yn=zn-z∗,Yi=exp⁡ciVyn+∑j=1saijVhψYj+z∗-VYj,i=1,…,s,yn+1=exp⁡Vyn+∑i=1sbiVhψYi+z∗-VYi,zn+1=yn+1+z∗,where V=hΩ with Ω=∂ψ/∂z(z∗).


## 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): (12)Yi=exp⁡ciVyn+∑j=1saijexp⁡ci-cjVhfYj-VYj,i=1,…,s,yn+1=exp⁡Vyn+∑i=1sbiexp⁡1-ciVhfYi-VYi,where ci, aij, and bi,i,j=1,…,s, are the constant coefficients of the RK method (5). Note that as Ω→0 (V→0), scheme (12) reduces to the RK method (5). The latter is called theprototype RK method of the former. In the sequel of this paper, scheme (12) will be referred to as anexpRK method. Its Butcher tableau can be written as follows:(13)
​
c
exp
⁡
c
V
e
x
p
c
V
T
A
exp
⁡
-
c
V
exp
⁡
V
b
T
exp
⁡
e
-
c
V
=
c
1
exp
⁡
c
1
V
a
11
a
12
exp
⁡
c
1
-
c
2
V
⋯
a
1
s
exp
⁡
c
1
-
c
s
V
⋮
⋮
⋮
⋮
⋱
⋮
c
s
exp
⁡
c
1
V
a
s
1
exp
⁡
c
s
-
c
1
V
a
s
2
exp
⁡
c
s
-
c
2
V
⋯
a
s
s
exp
⁡
V
b
1
exp
⁡
1
-
c
1
V
b
2
exp
⁡
1
-
c
2
V
⋯
b
s
exp
⁡
1
-
c
s
V
,where e=(1,…,1)T, exp⁡((e-c)V=(exp⁡((1-c1)V)T,…,exp⁡((1-cs)V)T)T, and exp⁡(cV)TAexp⁡(-cV)=(exp⁡((ci-cj)V)aij)s×s.In Kronecker’s notation, scheme (12) can be written as (14)Y=exp⁡cVyn⊗e+exp⁡cVA⊗Idexp⁡-cVhfY-VY,yn+1=exp⁡Vyn+exp⁡VbT⊗Idexp⁡-cVhfY-VY,where Id is the d×d unit matrix and exp⁡-cV(hf(Y)-VY)=(exp⁡((-c1)V)T(hf(Y1)-VY1),…,exp⁡((-cs)V)T(hf(Ys)-VYs))T.Among simple examples are the following:(a)
Theexponential Euler method (explicit), denoted by expEuler:(15)yn+1=exp⁡Vyn+exp⁡Vhfyn-Vyn(b)
Theexponential backward Euler method (implicit), denoted by expImEuler:(16)yn+1=exp⁡Vyn+hfyn+1-Vyn+1(c)
Theexponential trapezoidal rule (implicit), denoted by expTrap:(17)yn+1=expVyn+expV2hfyn-Vyn+12hfyn+1-Vyn+1(d)
Theexponential Heun rule (explicit), denoted by expHeun:(18)yp=exp⁡Vyn+exp⁡Vhfyn-Vyn,yc=exp⁡Vyn+hfyp-Vyp,yn+1=yp+yc2(e)
Theexponential midpoint rule (implicit), denoted by expMid:(19)yn+1=exp⁡Vyn+exp⁡V2hfexp⁡V/2yn+exp⁡-V/2yn+12-12VexpV2yn+exp-V2yn+1Two 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](20)
​
0
1
3
1
3
2
3
-
1
3
1
1
1
-
1
1
1
8
3
8
3
8
1
8
,
(21)
​
0
1
2
1
2
1
2
0
1
2
1
0
0
1
1
6
2
6
2
6
1
6The (algebraic) order is a measure of the accuracy of numerical method. A method is said to haveorderp if its local error LE=O(hp+1).Theorem 3.
The expRK method (12) has the same algebraic order as its prototype RK method.Proof.
The conclusion follows by expanding the exponential functionsexp⁡(ciV), exp⁡((ci-cj)V), and exp⁡((1-ci)V) in scheme (12) in series of V=hΩ, hence in series of h, 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 thatf in system (3) satisfies the Lipschitz condition and the origin is a steady state of the system; that is, f(0)=0. Then the origin is also a fixed point of the expRK method (12) for small step size h.The proof is given in AppendixA.In order to apply the expRK method (9) or (12) to system (2), we first use a coordinate transform u(t)=r(t)-r∗,v(t)=p(t)-p∗ to translate the steady state (r∗,p∗) of the system to the origin and yields(22)u˙t=-Γut+F′p∗vt+Gvt,v˙t=Kut-Mvt,where F′(p∗) is the Jacobian matrix of F(p) at point p∗ and G(v(t))=F(p∗+v(t))-F′(p∗)v(t)-F(p∗). Then system (2) can be written in the form (7) with y(t)=(u(t),v(t))T where the rate matrix (23)Ω=-ΓF′p∗K-Mand the function g(y)=(G(v)T,0)T.


## 3. Numerical Illustrations

From Theorem4, 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 by(24)r˙t=-γrt+Fpt,p˙t=-μpt+κrt,where F(p(t))=α/1+p(t)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])(25)r˙1=m1H+p2;θ2,n2-γ1r1,r˙2=m2H-p1;θ1,n1-γ2r2,p˙1=κ1r1-μ1p1,p˙2=κ2r2-μ2p2,where, for i=1,2, ri is the concentration of mRNA Ri produced by gene gi, pi is the concentration of protein Pi, mi is the maximal transcription rate of gene gi, κi is the translation rate of mRNA Ri, and γi and μi are the degradation rates of mRNA Ri and protein Pi, respectively. The functions(26)H+p2;θ2,n2=p2n2θ2n2+p2n2,H-p1;θ1,n1=θ1n1θ1n1+p1n1are theHill functions of activation and repression, respectively. The parameters θ1,θ2 are theexpression thresholds. The integer value of ni(i=1,2), called theHill coefficient, stands for the number of protein monomers required for saturation of binding to DNA. It is easy to see that the activation function H+ is increasing in p2 and the repression function H- is decreasing in p1.(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 S(t)=0, has the form(27)P˙I=sp+jaPA-dpPI-kcPIM+jcC,M˙=sm0+sm1PI+sm2PAPI+PA+Km+ku+jcC-dm+kcPIM,C˙=kcPIM-jc+kuC,P˙A=-ja+dpPA,where PI represents the concentration of the p53 tumor suppressor, M (mdm2) is the concentration of the p53’s main negative regulator, C is the concentration of the p53-mdm2 complex, PA is the concentration of an active form of p53 that is resistant against mdm2-mediated degradation, 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), dp is the degradation rate of active p53, and Km 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 1/θ2p∗3+p∗+ακ/γμ=0 and the relation r∗=μ/κp∗. If the system is written in the form (7), the rate matrix(28)Ω=-γ-2αp∗/θ21+p∗2/θ22κ-μ,and the function (29)gy=α1+y2+p∗2/θ2+2αp∗/θ21+p∗2/θ22y20,where y=(y1,y2)T=(r-r∗,p-p∗). With the parameter values (provided by [31])(30)α=3,γ=1,μ=1.5,κ=5,θ=1,this system has a unique positive steady state (r∗,p∗)=(0.6,2) where the rate matrix Ω has eigenvalues λ1,2=-1.2500±2.9767i, where i is the imaginary unit satisfying i2=-1. 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 (r(0),p(0))=(0.2,0.3),(0.2,0.8),(0.2,2.3), respectively. Figure 1(b) shows the time evolution of concentrations of mRNA and protein starting at (r(0),p(0))=(1.2,0.3).Figure 1
One-gene system: (a) solution trajectories; (b) time evolution of concentrations of mRNA and protein.
(a) (b)With the above values of parameters and initial data, we integrate system (25) on the time interval [0,100] 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 [50,100]. The numerical results are presented in Figures 2 and 3.Figure 2
Accuracy comparison for the one-gene system: (a) Euler and expEuler; (b) Heun and expHeun with step sizeh=1.2.
(a) (b)Figure 3
Accuracy comparison for the one-gene system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4 with step sizeh=1.6.
(a) (b)The system is solved on the time interval[0,100] with initial values of mRNA and protein r(0)=0.6 and p(0)=0.8 and with different step sizes. The numerical results are presented in Tables 1 and 2.Table 1
One-gene system: comparison of average errors for Euler, expEuler, Heun, and expEuler methods.
Step size Euler expEuler Heun expHeun 1 / 4 3.5849 × 1 0 - 2 1.1696 × 1 0 - 3 2.1384 × 1 0 - 2 1.5011 × 1 0 - 4 1 / 8 2.6231 × 1 0 - 2 5.2015 × 1 0 - 4 2.0930 × 1 0 - 2 5.7342 × 1 0 - 5 1 / 16 2.2715 × 1 0 - 2 2.3383 × 1 0 - 4 2.0539 × 1 0 - 2 3.8767 × 1 0 - 5 1 / 32 2.1259 × 1 0 - 2 9.8536 × 1 0 - 5 2.0284 × 1 0 - 2 3.5349 × 1 0 - 5Table 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 1 / 4 3.0743 × 1 0 - 2 3.1270 × 1 0 - 5 3.1713 × 1 0 - 2 2.2006 × 1 0 - 4 1 / 8 2.5309 × 1 0 - 2 3.4131 × 1 0 - 5 2.5490 × 1 0 - 2 3.4655 × 1 0 - 5 1 / 16 2.2533 × 1 0 - 2 3.4373 × 1 0 - 5 2.2569 × 1 0 - 2 3.4403 × 1 0 - 5 1 / 32 2.1219 × 1 0 - 2 3.4407 × 1 0 - 5 2.1227 × 1 0 - 2 3.4408 × 1 0 - 5


#### 3.1.2. The Two-Gene System

The steady statesy∗=(r1∗,r2∗,p1∗,p2∗)T of system (25) are determined by the equations(31)θ2n2p1∗p1∗n1+θ1n1n2+p1∗-m1k1γ1μ1m2k2γ2μ2θ1n1n2=0,(32)r1∗=μ1κ1p1∗,r2∗=μ2κ2p2∗,p2∗=m2κ2θ1n1γ2μ2θ1n1+p1∗n1.Putting in the form (7), system (25) has the rate matrix(33)Ω=∂f∂yy∗=-γ100m1n2θ2n2p2∗n2-1θ2n2+p2∗n220-γ2-m2n1θ1n1p1∗n1-1θ1n1+p1∗n120κ10-μ100κ20-μ2and the function(34)gy=m1p2n2θ2n2+y4+p2∗n2-γ1r1∗-m1n2θ2n2p2∗n2-1θ2n2+p2∗n22y4m2θ1n1θ1n1+y3+p1∗n1-γ2r2∗-m2n1θ1n1p1∗n1-1θ1n1+p1∗n12y300with y=(y1,y2,y3,y4)T=(r1-r1∗,r2-r2∗,p1-p1∗,p4-p4∗)T. The characteristic equation of the rate matrix Ω is(35)λ+γ1λ+γ2λ+μ1λ+μ2+D=0,where(36)D=m1m2κ1κ2n1n2p1∗n1-1p2∗n2-1θ1n1θ2n2p1∗n1+θ1n12p2∗n2+θ2n22.For a certain value of D=DHopf, the real part of one of the eigenvalues crosses zero, indicating a loss of stability through a Hopf bifurcation. For n1>1 and n2>1, [17] has calculated this value explicitly as (37)DHopf=γ1+γ2γ1+μ1γ2+μ1γ1+μ2γ2+μ2μ1+μ2γ1+γ2+μ1+μ22.In our experiment, we take the values of parameters as follows:(38)n1=n2=3,m1=m2=1.8,k1=k2=1,γ1=γ2=1,μ1=μ2=1,θ1=θ2=0.6542.Figure 4(a) shows three solution trajectories projected on the mRNA 1-protein 1 plane starting at (r1(0),r2(0),p1(0),p2(0))=(0.1,0.5,0,0), (0.3,0.3,0,0), and (0.5,0.2,0,0), respectively. Figure 4(b) shows the time evolution of concentrations of protein 1 and protein 2.Figure 4
Two-gene system: (a) phase trajectories projected on mRNA 1-protein 1 plane; (b) time evolution of concentrations of proteins.
(a) (b)We solve (31) for p1∗ by Newton’s iteration and then substitute it into (32) obtaining(39)r1∗=0.81471271066221,r2∗=0.61403210214378,p1∗=0.81471271066221,p2∗=0.61403210214378.The rate matrix (33) has the eigenvalues with negative real parts: (40)λ1,2=-1.94911366876016±0.94911366876016i,λ3,4=-0.05088633123984±0.94911366876016i.Therefore the steady state y∗ is asymptotically stable.With the values of parameters (38) and initial data and initial values r1(0)=0.6, r2(0)=0.8, p1(0)=0, and p2(0)=0, we integrate system (25) on the time interval [0,100] by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step sizes h=1.8,h=1.7,h=1.5,h=1.678, respectively. In Figures 5 and 6, we plot the global error growth of protein 1 on the time interval [50,100].Figure 5
Accuracy comparison for the two-gene system: (a) Euler and expEuler; (b) Heun and expHeun.
(a) (b)Figure 6
Accuracy comparison for the two-gene system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4.
(a) (b)In Tables3 and 4, average errors are compared for differential step sizes.Table 3
Two-gene system: comparison of average errors for Euler, expEuler, Heun, and expHeun methods.
Step size Euler expEuler Heun expHeun 1 1.6832 4.4456 × 1 0 - 3 2.9451 × 1 0 - 1 2.1871 × 1 0 - 3 1 / 2 9.21627 × 1 0 - 1 2.4275 × 1 0 - 3 4.0230 × 1 0 - 2 1.0848 × 1 0 - 3 1 / 4 5.1922 × 1 0 - 1 1.4028 × 1 0 - 3 1.0096×10-2 9.5000 × 1 0 - 4 1 / 8 1.8775 × 1 0 - 1 1.0103 × 1 0 - 3 2.7588 × 1 0 - 3 9.4565 × 1 0 - 4 1 / 16 5.0636 × 1 0 - 2 9.0239 × 1 0 - 4 1.1937 × 1 0 - 3 9.4830 × 1 0 - 4 1 / 32 1.8282 × 1 0 - 2 8.9416 × 1 0 - 4 9.8137 × 1 0 - 4 9.4942 × 1 0 - 4Table 4
Two-gene system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 1 4.6619 × 1 0 - 2 8.7331 × 1 0 - 4 8.6240 × 1 0 - 3 9.7906 × 1 0 - 4 1 / 2 3.7028 × 1 0 - 1 9.4577 × 1 0 - 4 5.2468 × 1 0 - 1 9.4457 × 1 0 - 4 1 / 4 3.2666 × 1 0 - 1 9.4805 × 1 0 - 4 3.6973 × 1 0 - 1 9.4766 × 1 0 - 4 1 / 8 1.3398 × 1 0 - 1 9.4894 × 1 0 - 4 1.4412 × 1 0 - 1 9.4891 × 1 0 - 4 1 / 16 4.4359 × 1 0 - 2 9.4947 × 1 0 - 4 4.5572 × 1 0 - 2 9.4947 × 1 0 - 4 1 / 32 1.7338 × 1 0 - 2 9.4974 × 1 0 - 4 1.7525 × 1 0 - 2 9.4974 × 1 0 - 4


#### 3.1.3. The p53-mdm2 System

The steady statey∗=(PI∗,M∗,C∗,PA∗) is determined by the following equations:(41)sp-dpPI∗-kckudmjc+kuPI∗sm0+sm1PI∗+KmPI∗=0,M∗=sm0dm+sm1PI∗dmPI∗+Km,C∗=kcjc+kuPI∗M∗,PA∗=0.Putting in the form (7), system (27) has the rate matrix(42)Ω=∂f∂yy∗=-dp-kcM∗-kcPI∗jcjasm1Km+sm1-sm2PA∗PI∗+PA∗+Km2-kcM∗-dm+kcPI∗ku+jcsm2Km+sm2-sm1PI∗PI∗+PA∗+Km2kcM∗kcPI∗-ku+jc0000ja+dp.and the function(43)gy=-kcy1y2g2ykcy1y20,where y=(y1,y2,y3,y4)T=(PI-PI∗,M-M∗,C-C∗,PA-PA∗)T, and(44)g2y=-sm1-sm2PA∗+sm1Kmy1+sm2-sm1PI∗+sm2Kmy4y1+y4y1+PI∗+y4+PA∗+KmPI∗+PA∗+Km2-kcy1y2.We use the parameter values (see [33]) as follows:(45)sm0=2×10-3nMmin-1,ka=20min-1,ja=0.2min-1,sm1=0.15nMmin-1,kc=4min-1nM-1,jc=2×10-3min-1,sm2=0.2nMmin-1,ku=0.4min-1,dm=0.4min-1,sp=1.4nMmin-1,Km=100nM,dp=2×10-4min-1.For simplicity, we take the small function S(t)≡0. The system has a unique steady state (PI∗,M∗,C∗,PA∗)=(9.42094,0.0372868,3.49529,0). The rate matrix Ω has the eigenvalues (46)λ1=-38.4766,λ2,3=-0.0028±0.0220i,λ4=-0.2002.Since all the eigenvalues have negative real parts, the steady state is asymptotically stable.Figure7(a) shows three solution trajectories projected on the inactive p53-complex plane starting at (PI(0),M0(0),C0(0),PA(0))=(5,0.01,0,0),(30,0.01,0,0),(15,0.01,0,0), respectively. Figure 7(b) shows the time evolution of concentrations of p53 and mdm2.Figure 7
p53-mdm2 system: (a) phase curves projected on the inactive p53-complex plane; (b) time evolution of concentrations of p53 and mdm2.
(a) (b)With the values of parameters in (45) and initial data and initial values PI(0)=10, M(0)=0.1, C(0)=3, and PA(0)=0,8, we integrate system (25) on the time interval [0,100] by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step size h=1/16. In Figures 8 and 9, we plot the global error growth of the inactive p53 on the time interval [50,100].Figure 8
Accuracy comparison for the p53-mdm2 system: (a) Euler and expEuler; (b) Heun and expHeun.
(a) (b)Figure 9
Accuracy comparison for the p53-mdm2 system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4.
(a) (b)The problem is solved on the interval[0,50] for differential step sizes and the average errors are presented in Tables 5 and 6.Table 5
p53-mdm2 system: comparison of average errors for Euler, expEuler, Heun, and expHeun methods.
Step size Euler expEuler Heun expHeun 1 / 16 3.5455 9.7152 × 1 0 - 3 3.3784 2.47061 × 1 0 - 3 1 / 32 1.3141 × 1 0 - 3 5.3681 × 1 0 - 3 1.7446 × 1 0 - 4 1.2767 × 1 0 - 3 1 / 64 6.5825 × 1 0 - 4 2.7553 × 1 0 - 3 5.9564 × 1 0 - 5 4.383 × 1 0 - 4 1 / 128 3.0970 × 1 0 - 4 1.3769 × 1 0 - 3 4.7031 × 1 0 - 5 1.1482 × 1 0 - 4Table 6
p53-mdm2 system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 1 / 16 3.1651 2.2126 × 1 0 - 4 3.2096 7.8709 × 1 0 - 4 1 / 32 1.2676 × 1 0 - 3 2.5953 × 1 0 - 5 1.2770 × 1 0 - 3 9.4086 × 1 0 - 5 1 / 64 6.4492 × 1 0 - 4 4.2016 × 1 0 - 5 6.4759 × 1 0 - 4 4.8393 × 1 0 - 5 1 / 128 3.0628 × 1 0 - 4 4.4244 × 1 0 - 5 3.0697 × 1 0 - 4 4.4697 × 1 0 - 5From Figures2, 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 [0,100] with the step sizes h=1/2j,j=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 h=1/2, expRK4 cannot produce smaller error, just increasing the computation effort.Figure 10
Efficiency curves with step sizesh=1/2j,j=1,2,3,4.Figure 11
RK4 simulation with step sizeh=1.32 compared with the exact solution.


## 3.1. Accuracy Test

### 3.1.1. The One-Gene System

Steady states of system (24) can be determined by the cubic equation 1/θ2p∗3+p∗+ακ/γμ=0 and the relation r∗=μ/κp∗. If the system is written in the form (7), the rate matrix(28)Ω=-γ-2αp∗/θ21+p∗2/θ22κ-μ,and the function (29)gy=α1+y2+p∗2/θ2+2αp∗/θ21+p∗2/θ22y20,where y=(y1,y2)T=(r-r∗,p-p∗). With the parameter values (provided by [31])(30)α=3,γ=1,μ=1.5,κ=5,θ=1,this system has a unique positive steady state (r∗,p∗)=(0.6,2) where the rate matrix Ω has eigenvalues λ1,2=-1.2500±2.9767i, where i is the imaginary unit satisfying i2=-1. 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 (r(0),p(0))=(0.2,0.3),(0.2,0.8),(0.2,2.3), respectively. Figure 1(b) shows the time evolution of concentrations of mRNA and protein starting at (r(0),p(0))=(1.2,0.3).Figure 1
One-gene system: (a) solution trajectories; (b) time evolution of concentrations of mRNA and protein.
(a) (b)With the above values of parameters and initial data, we integrate system (25) on the time interval [0,100] 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 [50,100]. The numerical results are presented in Figures 2 and 3.Figure 2
Accuracy comparison for the one-gene system: (a) Euler and expEuler; (b) Heun and expHeun with step sizeh=1.2.
(a) (b)Figure 3
Accuracy comparison for the one-gene system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4 with step sizeh=1.6.
(a) (b)The system is solved on the time interval[0,100] with initial values of mRNA and protein r(0)=0.6 and p(0)=0.8 and with different step sizes. The numerical results are presented in Tables 1 and 2.Table 1
One-gene system: comparison of average errors for Euler, expEuler, Heun, and expEuler methods.
Step size Euler expEuler Heun expHeun 1 / 4 3.5849 × 1 0 - 2 1.1696 × 1 0 - 3 2.1384 × 1 0 - 2 1.5011 × 1 0 - 4 1 / 8 2.6231 × 1 0 - 2 5.2015 × 1 0 - 4 2.0930 × 1 0 - 2 5.7342 × 1 0 - 5 1 / 16 2.2715 × 1 0 - 2 2.3383 × 1 0 - 4 2.0539 × 1 0 - 2 3.8767 × 1 0 - 5 1 / 32 2.1259 × 1 0 - 2 9.8536 × 1 0 - 5 2.0284 × 1 0 - 2 3.5349 × 1 0 - 5Table 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 1 / 4 3.0743 × 1 0 - 2 3.1270 × 1 0 - 5 3.1713 × 1 0 - 2 2.2006 × 1 0 - 4 1 / 8 2.5309 × 1 0 - 2 3.4131 × 1 0 - 5 2.5490 × 1 0 - 2 3.4655 × 1 0 - 5 1 / 16 2.2533 × 1 0 - 2 3.4373 × 1 0 - 5 2.2569 × 1 0 - 2 3.4403 × 1 0 - 5 1 / 32 2.1219 × 1 0 - 2 3.4407 × 1 0 - 5 2.1227 × 1 0 - 2 3.4408 × 1 0 - 5


### 3.1.2. The Two-Gene System

The steady statesy∗=(r1∗,r2∗,p1∗,p2∗)T of system (25) are determined by the equations(31)θ2n2p1∗p1∗n1+θ1n1n2+p1∗-m1k1γ1μ1m2k2γ2μ2θ1n1n2=0,(32)r1∗=μ1κ1p1∗,r2∗=μ2κ2p2∗,p2∗=m2κ2θ1n1γ2μ2θ1n1+p1∗n1.Putting in the form (7), system (25) has the rate matrix(33)Ω=∂f∂yy∗=-γ100m1n2θ2n2p2∗n2-1θ2n2+p2∗n220-γ2-m2n1θ1n1p1∗n1-1θ1n1+p1∗n120κ10-μ100κ20-μ2and the function(34)gy=m1p2n2θ2n2+y4+p2∗n2-γ1r1∗-m1n2θ2n2p2∗n2-1θ2n2+p2∗n22y4m2θ1n1θ1n1+y3+p1∗n1-γ2r2∗-m2n1θ1n1p1∗n1-1θ1n1+p1∗n12y300with y=(y1,y2,y3,y4)T=(r1-r1∗,r2-r2∗,p1-p1∗,p4-p4∗)T. The characteristic equation of the rate matrix Ω is(35)λ+γ1λ+γ2λ+μ1λ+μ2+D=0,where(36)D=m1m2κ1κ2n1n2p1∗n1-1p2∗n2-1θ1n1θ2n2p1∗n1+θ1n12p2∗n2+θ2n22.For a certain value of D=DHopf, the real part of one of the eigenvalues crosses zero, indicating a loss of stability through a Hopf bifurcation. For n1>1 and n2>1, [17] has calculated this value explicitly as (37)DHopf=γ1+γ2γ1+μ1γ2+μ1γ1+μ2γ2+μ2μ1+μ2γ1+γ2+μ1+μ22.In our experiment, we take the values of parameters as follows:(38)n1=n2=3,m1=m2=1.8,k1=k2=1,γ1=γ2=1,μ1=μ2=1,θ1=θ2=0.6542.Figure 4(a) shows three solution trajectories projected on the mRNA 1-protein 1 plane starting at (r1(0),r2(0),p1(0),p2(0))=(0.1,0.5,0,0), (0.3,0.3,0,0), and (0.5,0.2,0,0), respectively. Figure 4(b) shows the time evolution of concentrations of protein 1 and protein 2.Figure 4
Two-gene system: (a) phase trajectories projected on mRNA 1-protein 1 plane; (b) time evolution of concentrations of proteins.
(a) (b)We solve (31) for p1∗ by Newton’s iteration and then substitute it into (32) obtaining(39)r1∗=0.81471271066221,r2∗=0.61403210214378,p1∗=0.81471271066221,p2∗=0.61403210214378.The rate matrix (33) has the eigenvalues with negative real parts: (40)λ1,2=-1.94911366876016±0.94911366876016i,λ3,4=-0.05088633123984±0.94911366876016i.Therefore the steady state y∗ is asymptotically stable.With the values of parameters (38) and initial data and initial values r1(0)=0.6, r2(0)=0.8, p1(0)=0, and p2(0)=0, we integrate system (25) on the time interval [0,100] by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step sizes h=1.8,h=1.7,h=1.5,h=1.678, respectively. In Figures 5 and 6, we plot the global error growth of protein 1 on the time interval [50,100].Figure 5
Accuracy comparison for the two-gene system: (a) Euler and expEuler; (b) Heun and expHeun.
(a) (b)Figure 6
Accuracy comparison for the two-gene system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4.
(a) (b)In Tables3 and 4, average errors are compared for differential step sizes.Table 3
Two-gene system: comparison of average errors for Euler, expEuler, Heun, and expHeun methods.
Step size Euler expEuler Heun expHeun 1 1.6832 4.4456 × 1 0 - 3 2.9451 × 1 0 - 1 2.1871 × 1 0 - 3 1 / 2 9.21627 × 1 0 - 1 2.4275 × 1 0 - 3 4.0230 × 1 0 - 2 1.0848 × 1 0 - 3 1 / 4 5.1922 × 1 0 - 1 1.4028 × 1 0 - 3 1.0096×10-2 9.5000 × 1 0 - 4 1 / 8 1.8775 × 1 0 - 1 1.0103 × 1 0 - 3 2.7588 × 1 0 - 3 9.4565 × 1 0 - 4 1 / 16 5.0636 × 1 0 - 2 9.0239 × 1 0 - 4 1.1937 × 1 0 - 3 9.4830 × 1 0 - 4 1 / 32 1.8282 × 1 0 - 2 8.9416 × 1 0 - 4 9.8137 × 1 0 - 4 9.4942 × 1 0 - 4Table 4
Two-gene system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 1 4.6619 × 1 0 - 2 8.7331 × 1 0 - 4 8.6240 × 1 0 - 3 9.7906 × 1 0 - 4 1 / 2 3.7028 × 1 0 - 1 9.4577 × 1 0 - 4 5.2468 × 1 0 - 1 9.4457 × 1 0 - 4 1 / 4 3.2666 × 1 0 - 1 9.4805 × 1 0 - 4 3.6973 × 1 0 - 1 9.4766 × 1 0 - 4 1 / 8 1.3398 × 1 0 - 1 9.4894 × 1 0 - 4 1.4412 × 1 0 - 1 9.4891 × 1 0 - 4 1 / 16 4.4359 × 1 0 - 2 9.4947 × 1 0 - 4 4.5572 × 1 0 - 2 9.4947 × 1 0 - 4 1 / 32 1.7338 × 1 0 - 2 9.4974 × 1 0 - 4 1.7525 × 1 0 - 2 9.4974 × 1 0 - 4


### 3.1.3. The p53-mdm2 System

The steady statey∗=(PI∗,M∗,C∗,PA∗) is determined by the following equations:(41)sp-dpPI∗-kckudmjc+kuPI∗sm0+sm1PI∗+KmPI∗=0,M∗=sm0dm+sm1PI∗dmPI∗+Km,C∗=kcjc+kuPI∗M∗,PA∗=0.Putting in the form (7), system (27) has the rate matrix(42)Ω=∂f∂yy∗=-dp-kcM∗-kcPI∗jcjasm1Km+sm1-sm2PA∗PI∗+PA∗+Km2-kcM∗-dm+kcPI∗ku+jcsm2Km+sm2-sm1PI∗PI∗+PA∗+Km2kcM∗kcPI∗-ku+jc0000ja+dp.and the function(43)gy=-kcy1y2g2ykcy1y20,where y=(y1,y2,y3,y4)T=(PI-PI∗,M-M∗,C-C∗,PA-PA∗)T, and(44)g2y=-sm1-sm2PA∗+sm1Kmy1+sm2-sm1PI∗+sm2Kmy4y1+y4y1+PI∗+y4+PA∗+KmPI∗+PA∗+Km2-kcy1y2.We use the parameter values (see [33]) as follows:(45)sm0=2×10-3nMmin-1,ka=20min-1,ja=0.2min-1,sm1=0.15nMmin-1,kc=4min-1nM-1,jc=2×10-3min-1,sm2=0.2nMmin-1,ku=0.4min-1,dm=0.4min-1,sp=1.4nMmin-1,Km=100nM,dp=2×10-4min-1.For simplicity, we take the small function S(t)≡0. The system has a unique steady state (PI∗,M∗,C∗,PA∗)=(9.42094,0.0372868,3.49529,0). The rate matrix Ω has the eigenvalues (46)λ1=-38.4766,λ2,3=-0.0028±0.0220i,λ4=-0.2002.Since all the eigenvalues have negative real parts, the steady state is asymptotically stable.Figure7(a) shows three solution trajectories projected on the inactive p53-complex plane starting at (PI(0),M0(0),C0(0),PA(0))=(5,0.01,0,0),(30,0.01,0,0),(15,0.01,0,0), respectively. Figure 7(b) shows the time evolution of concentrations of p53 and mdm2.Figure 7
p53-mdm2 system: (a) phase curves projected on the inactive p53-complex plane; (b) time evolution of concentrations of p53 and mdm2.
(a) (b)With the values of parameters in (45) and initial data and initial values PI(0)=10, M(0)=0.1, C(0)=3, and PA(0)=0,8, we integrate system (25) on the time interval [0,100] by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step size h=1/16. In Figures 8 and 9, we plot the global error growth of the inactive p53 on the time interval [50,100].Figure 8
Accuracy comparison for the p53-mdm2 system: (a) Euler and expEuler; (b) Heun and expHeun.
(a) (b)Figure 9
Accuracy comparison for the p53-mdm2 system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4.
(a) (b)The problem is solved on the interval[0,50] for differential step sizes and the average errors are presented in Tables 5 and 6.Table 5
p53-mdm2 system: comparison of average errors for Euler, expEuler, Heun, and expHeun methods.
Step size Euler expEuler Heun expHeun 1 / 16 3.5455 9.7152 × 1 0 - 3 3.3784 2.47061 × 1 0 - 3 1 / 32 1.3141 × 1 0 - 3 5.3681 × 1 0 - 3 1.7446 × 1 0 - 4 1.2767 × 1 0 - 3 1 / 64 6.5825 × 1 0 - 4 2.7553 × 1 0 - 3 5.9564 × 1 0 - 5 4.383 × 1 0 - 4 1 / 128 3.0970 × 1 0 - 4 1.3769 × 1 0 - 3 4.7031 × 1 0 - 5 1.1482 × 1 0 - 4Table 6
p53-mdm2 system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 1 / 16 3.1651 2.2126 × 1 0 - 4 3.2096 7.8709 × 1 0 - 4 1 / 32 1.2676 × 1 0 - 3 2.5953 × 1 0 - 5 1.2770 × 1 0 - 3 9.4086 × 1 0 - 5 1 / 64 6.4492 × 1 0 - 4 4.2016 × 1 0 - 5 6.4759 × 1 0 - 4 4.8393 × 1 0 - 5 1 / 128 3.0628 × 1 0 - 4 4.4244 × 1 0 - 5 3.0697 × 1 0 - 4 4.4697 × 1 0 - 5From Figures2, 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.1.1. The One-Gene System

Steady states of system (24) can be determined by the cubic equation 1/θ2p∗3+p∗+ακ/γμ=0 and the relation r∗=μ/κp∗. If the system is written in the form (7), the rate matrix(28)Ω=-γ-2αp∗/θ21+p∗2/θ22κ-μ,and the function (29)gy=α1+y2+p∗2/θ2+2αp∗/θ21+p∗2/θ22y20,where y=(y1,y2)T=(r-r∗,p-p∗). With the parameter values (provided by [31])(30)α=3,γ=1,μ=1.5,κ=5,θ=1,this system has a unique positive steady state (r∗,p∗)=(0.6,2) where the rate matrix Ω has eigenvalues λ1,2=-1.2500±2.9767i, where i is the imaginary unit satisfying i2=-1. 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 (r(0),p(0))=(0.2,0.3),(0.2,0.8),(0.2,2.3), respectively. Figure 1(b) shows the time evolution of concentrations of mRNA and protein starting at (r(0),p(0))=(1.2,0.3).Figure 1
One-gene system: (a) solution trajectories; (b) time evolution of concentrations of mRNA and protein.
(a) (b)With the above values of parameters and initial data, we integrate system (25) on the time interval [0,100] 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 [50,100]. The numerical results are presented in Figures 2 and 3.Figure 2
Accuracy comparison for the one-gene system: (a) Euler and expEuler; (b) Heun and expHeun with step sizeh=1.2.
(a) (b)Figure 3
Accuracy comparison for the one-gene system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4 with step sizeh=1.6.
(a) (b)The system is solved on the time interval[0,100] with initial values of mRNA and protein r(0)=0.6 and p(0)=0.8 and with different step sizes. The numerical results are presented in Tables 1 and 2.Table 1
One-gene system: comparison of average errors for Euler, expEuler, Heun, and expEuler methods.
Step size Euler expEuler Heun expHeun 1 / 4 3.5849 × 1 0 - 2 1.1696 × 1 0 - 3 2.1384 × 1 0 - 2 1.5011 × 1 0 - 4 1 / 8 2.6231 × 1 0 - 2 5.2015 × 1 0 - 4 2.0930 × 1 0 - 2 5.7342 × 1 0 - 5 1 / 16 2.2715 × 1 0 - 2 2.3383 × 1 0 - 4 2.0539 × 1 0 - 2 3.8767 × 1 0 - 5 1 / 32 2.1259 × 1 0 - 2 9.8536 × 1 0 - 5 2.0284 × 1 0 - 2 3.5349 × 1 0 - 5Table 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 1 / 4 3.0743 × 1 0 - 2 3.1270 × 1 0 - 5 3.1713 × 1 0 - 2 2.2006 × 1 0 - 4 1 / 8 2.5309 × 1 0 - 2 3.4131 × 1 0 - 5 2.5490 × 1 0 - 2 3.4655 × 1 0 - 5 1 / 16 2.2533 × 1 0 - 2 3.4373 × 1 0 - 5 2.2569 × 1 0 - 2 3.4403 × 1 0 - 5 1 / 32 2.1219 × 1 0 - 2 3.4407 × 1 0 - 5 2.1227 × 1 0 - 2 3.4408 × 1 0 - 5


## 3.1.2. The Two-Gene System

The steady statesy∗=(r1∗,r2∗,p1∗,p2∗)T of system (25) are determined by the equations(31)θ2n2p1∗p1∗n1+θ1n1n2+p1∗-m1k1γ1μ1m2k2γ2μ2θ1n1n2=0,(32)r1∗=μ1κ1p1∗,r2∗=μ2κ2p2∗,p2∗=m2κ2θ1n1γ2μ2θ1n1+p1∗n1.Putting in the form (7), system (25) has the rate matrix(33)Ω=∂f∂yy∗=-γ100m1n2θ2n2p2∗n2-1θ2n2+p2∗n220-γ2-m2n1θ1n1p1∗n1-1θ1n1+p1∗n120κ10-μ100κ20-μ2and the function(34)gy=m1p2n2θ2n2+y4+p2∗n2-γ1r1∗-m1n2θ2n2p2∗n2-1θ2n2+p2∗n22y4m2θ1n1θ1n1+y3+p1∗n1-γ2r2∗-m2n1θ1n1p1∗n1-1θ1n1+p1∗n12y300with y=(y1,y2,y3,y4)T=(r1-r1∗,r2-r2∗,p1-p1∗,p4-p4∗)T. The characteristic equation of the rate matrix Ω is(35)λ+γ1λ+γ2λ+μ1λ+μ2+D=0,where(36)D=m1m2κ1κ2n1n2p1∗n1-1p2∗n2-1θ1n1θ2n2p1∗n1+θ1n12p2∗n2+θ2n22.For a certain value of D=DHopf, the real part of one of the eigenvalues crosses zero, indicating a loss of stability through a Hopf bifurcation. For n1>1 and n2>1, [17] has calculated this value explicitly as (37)DHopf=γ1+γ2γ1+μ1γ2+μ1γ1+μ2γ2+μ2μ1+μ2γ1+γ2+μ1+μ22.In our experiment, we take the values of parameters as follows:(38)n1=n2=3,m1=m2=1.8,k1=k2=1,γ1=γ2=1,μ1=μ2=1,θ1=θ2=0.6542.Figure 4(a) shows three solution trajectories projected on the mRNA 1-protein 1 plane starting at (r1(0),r2(0),p1(0),p2(0))=(0.1,0.5,0,0), (0.3,0.3,0,0), and (0.5,0.2,0,0), respectively. Figure 4(b) shows the time evolution of concentrations of protein 1 and protein 2.Figure 4
Two-gene system: (a) phase trajectories projected on mRNA 1-protein 1 plane; (b) time evolution of concentrations of proteins.
(a) (b)We solve (31) for p1∗ by Newton’s iteration and then substitute it into (32) obtaining(39)r1∗=0.81471271066221,r2∗=0.61403210214378,p1∗=0.81471271066221,p2∗=0.61403210214378.The rate matrix (33) has the eigenvalues with negative real parts: (40)λ1,2=-1.94911366876016±0.94911366876016i,λ3,4=-0.05088633123984±0.94911366876016i.Therefore the steady state y∗ is asymptotically stable.With the values of parameters (38) and initial data and initial values r1(0)=0.6, r2(0)=0.8, p1(0)=0, and p2(0)=0, we integrate system (25) on the time interval [0,100] by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step sizes h=1.8,h=1.7,h=1.5,h=1.678, respectively. In Figures 5 and 6, we plot the global error growth of protein 1 on the time interval [50,100].Figure 5
Accuracy comparison for the two-gene system: (a) Euler and expEuler; (b) Heun and expHeun.
(a) (b)Figure 6
Accuracy comparison for the two-gene system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4.
(a) (b)In Tables3 and 4, average errors are compared for differential step sizes.Table 3
Two-gene system: comparison of average errors for Euler, expEuler, Heun, and expHeun methods.
Step size Euler expEuler Heun expHeun 1 1.6832 4.4456 × 1 0 - 3 2.9451 × 1 0 - 1 2.1871 × 1 0 - 3 1 / 2 9.21627 × 1 0 - 1 2.4275 × 1 0 - 3 4.0230 × 1 0 - 2 1.0848 × 1 0 - 3 1 / 4 5.1922 × 1 0 - 1 1.4028 × 1 0 - 3 1.0096×10-2 9.5000 × 1 0 - 4 1 / 8 1.8775 × 1 0 - 1 1.0103 × 1 0 - 3 2.7588 × 1 0 - 3 9.4565 × 1 0 - 4 1 / 16 5.0636 × 1 0 - 2 9.0239 × 1 0 - 4 1.1937 × 1 0 - 3 9.4830 × 1 0 - 4 1 / 32 1.8282 × 1 0 - 2 8.9416 × 1 0 - 4 9.8137 × 1 0 - 4 9.4942 × 1 0 - 4Table 4
Two-gene system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 1 4.6619 × 1 0 - 2 8.7331 × 1 0 - 4 8.6240 × 1 0 - 3 9.7906 × 1 0 - 4 1 / 2 3.7028 × 1 0 - 1 9.4577 × 1 0 - 4 5.2468 × 1 0 - 1 9.4457 × 1 0 - 4 1 / 4 3.2666 × 1 0 - 1 9.4805 × 1 0 - 4 3.6973 × 1 0 - 1 9.4766 × 1 0 - 4 1 / 8 1.3398 × 1 0 - 1 9.4894 × 1 0 - 4 1.4412 × 1 0 - 1 9.4891 × 1 0 - 4 1 / 16 4.4359 × 1 0 - 2 9.4947 × 1 0 - 4 4.5572 × 1 0 - 2 9.4947 × 1 0 - 4 1 / 32 1.7338 × 1 0 - 2 9.4974 × 1 0 - 4 1.7525 × 1 0 - 2 9.4974 × 1 0 - 4


## 3.1.3. The p53-mdm2 System

The steady statey∗=(PI∗,M∗,C∗,PA∗) is determined by the following equations:(41)sp-dpPI∗-kckudmjc+kuPI∗sm0+sm1PI∗+KmPI∗=0,M∗=sm0dm+sm1PI∗dmPI∗+Km,C∗=kcjc+kuPI∗M∗,PA∗=0.Putting in the form (7), system (27) has the rate matrix(42)Ω=∂f∂yy∗=-dp-kcM∗-kcPI∗jcjasm1Km+sm1-sm2PA∗PI∗+PA∗+Km2-kcM∗-dm+kcPI∗ku+jcsm2Km+sm2-sm1PI∗PI∗+PA∗+Km2kcM∗kcPI∗-ku+jc0000ja+dp.and the function(43)gy=-kcy1y2g2ykcy1y20,where y=(y1,y2,y3,y4)T=(PI-PI∗,M-M∗,C-C∗,PA-PA∗)T, and(44)g2y=-sm1-sm2PA∗+sm1Kmy1+sm2-sm1PI∗+sm2Kmy4y1+y4y1+PI∗+y4+PA∗+KmPI∗+PA∗+Km2-kcy1y2.We use the parameter values (see [33]) as follows:(45)sm0=2×10-3nMmin-1,ka=20min-1,ja=0.2min-1,sm1=0.15nMmin-1,kc=4min-1nM-1,jc=2×10-3min-1,sm2=0.2nMmin-1,ku=0.4min-1,dm=0.4min-1,sp=1.4nMmin-1,Km=100nM,dp=2×10-4min-1.For simplicity, we take the small function S(t)≡0. The system has a unique steady state (PI∗,M∗,C∗,PA∗)=(9.42094,0.0372868,3.49529,0). The rate matrix Ω has the eigenvalues (46)λ1=-38.4766,λ2,3=-0.0028±0.0220i,λ4=-0.2002.Since all the eigenvalues have negative real parts, the steady state is asymptotically stable.Figure7(a) shows three solution trajectories projected on the inactive p53-complex plane starting at (PI(0),M0(0),C0(0),PA(0))=(5,0.01,0,0),(30,0.01,0,0),(15,0.01,0,0), respectively. Figure 7(b) shows the time evolution of concentrations of p53 and mdm2.Figure 7
p53-mdm2 system: (a) phase curves projected on the inactive p53-complex plane; (b) time evolution of concentrations of p53 and mdm2.
(a) (b)With the values of parameters in (45) and initial data and initial values PI(0)=10, M(0)=0.1, C(0)=3, and PA(0)=0,8, we integrate system (25) on the time interval [0,100] by the methods expEuler, expHeun, expRK3/8, and expRK4 as well as their corresponding prototype methods with the step size h=1/16. In Figures 8 and 9, we plot the global error growth of the inactive p53 on the time interval [50,100].Figure 8
Accuracy comparison for the p53-mdm2 system: (a) Euler and expEuler; (b) Heun and expHeun.
(a) (b)Figure 9
Accuracy comparison for the p53-mdm2 system: (a) RK3/8 and expRK3/8; (b) RK4 and expRK4.
(a) (b)The problem is solved on the interval[0,50] for differential step sizes and the average errors are presented in Tables 5 and 6.Table 5
p53-mdm2 system: comparison of average errors for Euler, expEuler, Heun, and expHeun methods.
Step size Euler expEuler Heun expHeun 1 / 16 3.5455 9.7152 × 1 0 - 3 3.3784 2.47061 × 1 0 - 3 1 / 32 1.3141 × 1 0 - 3 5.3681 × 1 0 - 3 1.7446 × 1 0 - 4 1.2767 × 1 0 - 3 1 / 64 6.5825 × 1 0 - 4 2.7553 × 1 0 - 3 5.9564 × 1 0 - 5 4.383 × 1 0 - 4 1 / 128 3.0970 × 1 0 - 4 1.3769 × 1 0 - 3 4.7031 × 1 0 - 5 1.1482 × 1 0 - 4Table 6
p53-mdm2 system: comparison of average errors for RK3/8, expRK3/8, RK4, and expRK4 methods.
Step size RK3/8 expRK3/8 RK4 expRK4 1 / 16 3.1651 2.2126 × 1 0 - 4 3.2096 7.8709 × 1 0 - 4 1 / 32 1.2676 × 1 0 - 3 2.5953 × 1 0 - 5 1.2770 × 1 0 - 3 9.4086 × 1 0 - 5 1 / 64 6.4492 × 1 0 - 4 4.2016 × 1 0 - 5 6.4759 × 1 0 - 4 4.8393 × 1 0 - 5 1 / 128 3.0628 × 1 0 - 4 4.4244 × 1 0 - 5 3.0697 × 1 0 - 4 4.4697 × 1 0 - 5From Figures2, 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 [0,100] with the step sizes h=1/2j,j=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 h=1/2, expRK4 cannot produce smaller error, just increasing the computation effort.Figure 10
Efficiency curves with step sizesh=1/2j,j=1,2,3,4.Figure 11
RK4 simulation with step sizeh=1.32 compared with the exact solution.


## 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 h=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 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|r|≪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 (r=1), 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 p and p+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.



---

*Source: 2729683-2017-01-19.xml*