# Limit-Cycle-Preserving Simulation of Gene Regulatory Oscillators

Authors: Xiong You
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

In order to simulate gene regulatory oscillators more effectively, Runge-Kutta (RK) integrators are adapted to the limit-cycle structure of the system. Taking into account the oscillatory feature of the gene regulatory oscillators, phase-fitted and amplification-fitted Runge-Kutta (FRK) methods are designed. New FRK methods with phase-fitted and amplification-fitted updated are also considered. The error coefficients and the error constant for each of new FRK methods are obtained. In the numerical simulation of the two-gene regulatory system, the new methods are shown to be more accurate and more efficient than their prototype RK methods in the long-term integration. It is a new discovery that the best fitting frequency not only depends on the problem to be solved, but also depends on the method.


---

## Body

## 1. Introduction

Differential equations (DEs) have become one of the powerful tools for modeling the complex dynamics of gene regulatory systems, where the cellular concentrations of mRNAs, proteins, and other molecules are assumed to vary continuously in time (see, e.g., de Jong [1], Widder et al. [2], Polynikis et al. [3],  Altinok et al. [4],  Gérard and Goldbeter [5], and the references therein). An N-gene activation-inhibition network can be modeled by a system of ordinary differential equations of the form (Polynikis et al. [3])
(1.1)Transcription:r˙i=Ψi(g1R(p1),…,gNR(pN))-γiri,Translation:p˙i=giP(ri)-δipi,
where for i=1,…,N, ri is the concentration of mRNA Ri produced by gene Gi, pi is the concentration of protein Pi produced by mRNA Ri, γi is the degradation rate of Ri, and δi is the degradation rate of Pi. Function giP(ri) is the translation function. Function Ψi(g1R(p1),…,gNR(pN)) is the regulation function, typically defined sums and products of functions g1R(p1),…,gNR(pN). If protein Pj has no effect on gene Gi, gjR(pj) does not appear in Ψi. If ∂Ψi/∂pj>0, protein pj is an activator of gene Gi. If ∂Ψi/∂pj<0, protein pj is an inhibitor of gene Gi.Since the functionsgiP(ri), giR(pi), and Ψi(·) are nonlinear, the analytical solution of the system (1.1) is in general not acquirable. Therefore, in order to reveal the dynamics of such gene regulatory systems, one usually resorts to numerical simulation. To be theoretically clear, we first consider abstractly the initial value problem (IVP) of the autonomous system of ordinary differential equations
(1.2)y˙=f(y),t>0,
where y:[0,+∞)→ℝd, “y˙” represents the first derivative of y with respect to time, and f:ℝd→ℝd is a sufficiently smooth function. Based on experimental observation of biological facts about gene regulation systems, it is reasonable to make the following assumptions.(I)
The system (1.2) has a stable limit cycle Γ0.(II)
The functionf(y) satisfies f(y*)=0, that is, y* is an equilibrium point of the system (1.2).(III)
The equilibrium pointy* lies inside the limit cycle Γ0 and there is no other equilibrium point inside Γ0.With these assumptions, any solution near the limit cycle is oscillatory. Up till now, differential equations of gene regulatory systems are mostly simulated by Runge-Kutta methods, especially by the classical fourth-order Runge-Kutta method or by the Runge-Kutta-Felhberg adaptive method as in the MATLAB package (see [6–8]). Unfortunately, the general-purpose RK method has not taken into account the special structure of the system (1.1) and hence they cannot give satisfactory numerical results. There are mainly two deficiencies of the classical RK methods: (i) They cannot produce as accurate numerical results as required, even if they have a very high algebraic order; (ii) the true dynamical behavior of the system cannot be preserved as expected in long term integration.Recently, researchers have proposed to adapt traditional integrators to problems whose solutions are oscillatory or periodic (see [9–12]). Bettis [13] constructed a three-stage method and a four-stage method which can solve the equation y'=iωy(i2=-1) exactly. Paternoster [14] developed a family of implicit Runge-Kutta (RK) and Runge-Kutta-Nyström (RKN) methods by means of trigonometric fitting. For oscillatory problems which can be put in the form of a second-order equation y′′=f(x,y), Franco [15] improved the update of the classical RKN methods and proposes a family of explicit RKN methods adapted to perturbed oscillators (ARKN) and a class of explicit adapted RK methods in [16]. Anastassi and Simos [17] constructed a phase-fitted and amplification-fitted RK method of “almost” order five. Chen et al. [18] considered symmetric and symplectic extended Runge-Kutta-Nyström methods. You et al. [19] investigated trigonometrically fitted Scheifele two-step (TFSTS) methods, derived the necessary and sufficient conditions for TFSTS methods of up to order five based on the linear operator theory and constructed two practical methods of algebraic four and five, respectively. For other important work on frequency-dependent integrators for general second-order oscillatory equations, the reader is referred to [20–24]. Some authors aimed at obtaining effective integrators for specific categories of oscillatory problems. For instance, Vigo-Aguiar and Simos [24] constructed an exponentially fitted and trigonometrically fitted method for orbital problems. The papers [25–27] have designed highly efficient integrators for the Schrödinger equation.The objective of this paper is to develop effective integrators for simulating the gene regulation system (1.1) near its limit cycle. Motivated by the work of Van de Vyver [28] on the two-step hybrid methods (FTSH), this paper develops a novel type of phase-amplification-fitted methods of Runge-Kutta type. These new numerical integrators respect the limit cycle structure of the system and are expected to be more accurate than the traditional RK methods in the long-time integration of gene regulatory systems. In Section 2, we present the order conditions for the modified Runge-Kutta methods for solving initial value problems of autonomous ordinary differential equations with oscillatory solutions. Section 3 gives the conditions for a modified RK method and its update to be phase-fitted and amplification-fitted. In Section 4, we construct six phase-fitted and amplification-fitted RK (FRK) methods. For each of these methods, the error coefficients and the error constant are presented. In Section 5, the two-gene regulation system is solved by the new FRK methods as well as their prototype RK methods. Section 6 is devoted to the conclusions.


## 2. Modified Runge-Kutta Methods and Order Conditions

Assume that the principal frequency of the problem (1.2) is known or can be accurately estimated in advance. This estimated frequency matrix is denoted by ω, which is also called the fitting frequency. We consider the following s-stage modified Runge-Kutta method:
(2.1)ki=f(yn+h∑j=1s‍aijkj),i=1,…,s,yn+1=yn+h∑i=1s‍bi(ν)ki,
where h is the step size, aij, i, j=1,…,s are real numbers and bi(ν), i=1,…,s are even functions of ν=hΩ. The scheme (2.1) can be represented by the Butcher tableau
(2.2)cAbT(ν)=c1a11⋯a1s⋮⋮⋱⋮csas1⋯assb1(ν)⋯bs(ν)
or simply by (c,A,b(ν)). With the rooted tree theory of Butcher [6], the exact solution to the problem (1.2) and the numerical solution produced by the modified RK type method (2.1) have the B-series expressions as follows:
(2.3)y(x0+h)=∑j=0∞‍hjj!∑ρ(τ)=j‍α(τ)ℱ(τ)(y0),y1=∑j=0∞‍hjj!∑ρ(τ)=j‍α(τ)γ(τ)bTΦ(τ)ℱ(τ)(y0),
where the trees τ, the functions ρ(τ) (order), α(τ) (the number of monotonic labelings of tree τ), γ(τ) (density), Φ(τ)=(Φ1(τ),…,Φs(τ))T (the vector of elementary weights) and the elementary differentials ℱ(τ)(y0) are defined in [8]. Then the local truncation error has the following series expansion:
(2.4)LTE=y(x0+h)-y1=∑j=0∞‍hjj!∑ρ(τ)=j‍α(τ)(1-γ(τ)bT(ν)Φ(τ))ℱ(τ)(y0).
If for any (p+1)th differentiable function f(y), when the scheme (2.1) is applied to the problem (1.2), LTE=Cp+1hp+1+𝒪(hp+2),Cp+1≠0, then the method (2.1) is said to have (algebraic) orderp.Theorem 2.1 (Franco [16]).
The modified RK type method (2.1) has order p if and only if the following conditions are satisfied:
(2.5)bTΦ(τ)=1γ(τ)+𝒪(hp-ρ(τ)+1),ρ(τ)=1,…,p.If the method (2.1) is of order p, then we have
(2.6)LTE=hp+1(p+1)!∑ρ(τ)=p+1‍α(τ)(1-γ(τ)bT(ν)Φ(τ))ℱ(τ)(y0)+𝒪(hp+2):=hp+1(p+1)!∑ρ(τ)=p+1‍α(τ)ϵ(τ)ℱ(τ)(y0)+𝒪(hp+2),
where ϵ(τ)=1-γ(τ)bT(ν)Φ(τ) is called the error coefficient associated with the tree τ of order p+1 which is involved in the leading term of the local truncation error. Denote
(2.7)ECp+1(ν)=(∑ρ(τ)=p+1‍ϵ(τ)2)1/2.
The positive number
(2.8)Cp+1=limν→0ECp+1(ν)
is called the error constant of the method. According to Theorem 2.1, we list the conditions for modified RK type methods to be of up to order 5 as follows:
(2.9)(i)Order1:bTe=1+𝒪(h);(ii)Order2:bTe=1+𝒪(h2),bTc=12+𝒪(h);(iii)Order3:bTe=1+𝒪(h3),bTc=12+𝒪(h2),bTc2=13+𝒪(h),bTAc=16+𝒪(h);(iv)Order4:bTe=1+𝒪(h4),bTc=12+𝒪(h3),bTc2=13+𝒪(h2),bTAc=16+𝒪(h2),bTc3=14+𝒪(h),bT(c·Ac)=18+𝒪(h),bT(Ac2)=112+𝒪(h),bTA2c=124+𝒪(h);(v)Order5:bTe=1+𝒪(h5),bTc=12+𝒪(h4),bTc2=13+𝒪(h3),bTAc=16+𝒪(h3),bTc3=14+𝒪(h2),bT(c·Ac)=18+𝒪(h2),bT(Ac2)=112+𝒪(h2),bTA2c=124+𝒪(h2),bTc4=15+𝒪(h),bTA2c2=110+𝒪(h),bT(c·(Ac2))=115+𝒪(h),bT(c·A2c)=130+𝒪(h),bT(Ac·Ac)=120,bTAc3=120+𝒪(h),bTA(c·Ac)=140+𝒪(h),bTA2c2=160+𝒪(h),bTA3c=1120+𝒪(h),
where e=(1,1,…,1)T, c=(c1,…,cs)T with ci=∑j=1saij,i=1,…,s, the dot “·” between two vectors indicates componentwise product, c2 indicates the componentwise square of the vector c and so on. We will follow this convention in the sequel.It is obvious that if a method satisfies the conditions for some orderp, then it satisfies all the conditions for orders lower than p.In Section4, when we use the above order conditions as equations to solve for the coefficients of a method, the higher order terms will be neglected.


## 3. Phase-Fitting and Amplification-Fitting Conditions

For oscillatory problems, Paternoster [14] and Van der Houwen and Sommeijer [29] propose to analyze the phase properties (dispersion and dissipation) of numerical integrators. To this end, we consider the linear equation
(3.1)y′=iωy.
The exact solution of this equation with the initial value y(x0)=y0 satisfies
(3.2)y(x0+h)=R0(z)y0,
where R0(z)=exp(z),z=iν. This means that after a period of time h, the exact solution experiences a phase advance ν=hω and the amplification remains constant.An application of the modified RK method (2.1) to (3.1) yields
(3.3)y1=R(z)y0,
where
(3.4)R(z)=1+zbT(I-zA)-1e,z=iν,e=(1,…,1)T.
Thus, after a time step h, the numerical solution obtains a phase advance argR(z) and the amplification factor |R(iν)|. R(z) is called the stability function of the method (2.1). Denoting the real and imaginary part of R(z) by U(ν) and V(ν), respectively, we have
(3.5)U(ν)=1-ν2(bTAe)+ν4(bTA3e)-ν6(bTA5e)+ν8(bTA7e)-⋯,V(ν)=ν(bTe)-ν3(bTA2e)+ν5(bTA4e)-ν7(bTA6e)+⋯.
Thus, argR(z)=tan-1(V(ν)/U(ν)) and |R(z)|=U2(ν)+V2(ν). The above analysis leads to the following definition.Definition 3.1 (see [29]).
The quantities(3.6)P(ν)=ν-argR(iν),D(ν)=1-|R(iν)|
are called the phase lag (or dispersion) and the error of amplification factor (or dissipation) of the method, respectively. If
(3.7)P(ν)=𝒪(νq+1),D(ν)=𝒪(νp+1),
then the method is called  dispersive of order q and  dissipative of order p, respectively. If
(3.8)P(ν)=0,D(ν)=0,
the method is called phase-fitted (or zero-dispersive) and amplification-fitted (or zero-dissipative), respectively. By a phase-amplification-fitted RK method we mean a modified RK method that is both phase-fitted and amplification-fitted, which we refer to as an FRK method.It is interesting to consider the phase properties of the update of the scheme (2.1). Suppose that the internal stages have been exact for the linear equation (3.1), that is, Yi=exp(iciV)y0, then the update gives
(3.9)y1=Ru(z)y0,
where
(3.10)Ru(z)=1+z∑i=1s‍bi(ν)exp(ciz),z=iν.
Denote the real and imaginary part of Ru(z) by Uu(ν) and νu(ν), respectively. Then, for small h,
(3.11)Uu(ν)=1-ν∑i=1s‍bi(ν)sin(ciν),Vu(ν)=ν∑i=1s‍bi(ν)cos(ciν).Definition 3.2.
The quantities(3.12)Pu(ν)=ν-argRu(iν),Du(ν)=1-|Ru(iν)|
are called the dispersion (or phase lag) and the dissipation (or error of amplification factor) of the update of the method, respectively. If
(3.13)Pu(ν)=𝒪(νq+1),Du(ν)=𝒪(νp+1),
then the update is called dispersive of orderq and  dissipative of orderp, respectively. If
(3.14)Pu(ν)=0,Du(ν)=0,
the update is called phase-fitted (or zero-dispersive) and amplification-fitted (or zero-dissipative), respectively.In general, a classical RK method with constant coefficients carries a nonzero dispersion and a nonzero dissipation when applied to the linear oscillatory equation (3.1). Therefore, they are neither phase-fitted nor amplification-fitted. So does the update. For example, the classical RK method of order four with constant coefficients (denoted as RK4, see [8]) given by
(3.15)0131323-12111-1116262616
has a phase lag and an error of amplification factor
(3.16)P(ν)=1120ν5+𝒪(ν7),D(ν)=1144ν6+𝒪(ν8).
For the update,
(3.17)Pu(ν)=-12880ν5+𝒪(ν7),Du(ν)=-15760ν6+𝒪(ν8).
Therefore, both the method and its update are dispersive of order four and dissipative of order five.The following theorem gives the necessary and sufficient conditions for a modified RK method and its update to be phase-fitted and amplification-fitted, respectively. The proof of this theorem is immediate.Theorem 3.3.
(i) The method (2.1) is phase-amplification-fitted if and only if
(3.18)U(ν)=cos(ν),V(ν)=sin(ν).
(ii) The update of the method (2.1) is phase-amplification-fitted if and only if
(3.19)Uu(ν)=cos(ν),Vu(ν)=sin(ν).


## 4. Construction of FRK Methods

In this section, we construct modified RK type methods that are phase-amplification-fitted based on the internal coefficients of familiar classical RK methods. For convenience we restrict ourselves to explicit methods.


### 4.1. FRK Methods of Type I

The coefficients of Type I methods are obtained by solving the phase-amplification-fitting conditions and some low-order conditions.


#### 4.1.1. Three-Stage Methods of Order Three

We begin by considering the following three-stage method(4.1)012121-12b1(ν)b2(ν)b3(ν)
The phase-fitting and amplification-fitting conditions (3.18) for the scheme (4.1) have the form
(4.2)1-12(b2+2b3)ν2=cos(ν),(b1+b2+b3)-b3ν2=sin(ν)ν.
On the other hand, the first-order condition is given by
(4.3)bTe=b1(ν)+b2(ν)+b3(ν)=1.
Solving (4.2) and (4.3), we obtain
(4.4)b1(ν)=-ν-ν3-2νcos(ν)+sin(ν)ν3,b2(ν)=-2νcos(ν)+2sin(ν)ν3,b3(ν)=ν-sin(ν)ν3.
As ν→0, bi(ν) have the following Taylor expansions
(4.5)b1(ν)=16+340ν2-1350400ν4+17362880ν6-11900800ν8+⋯,b2(ν)=23-115ν2+1420ν4-122680ν6+11995840ν8+⋯,b3(ν)=16-1120ν2+15040ν4-1362880ν6+139916800ν8+⋯.It is easy to verify that the method defined by (4.1) and (4.4) is of order three and we denote it by FRK3s3aI.The error coefficients for FRK3s3aI are given as follows:(4.6)ϵ(τ41)=1-8ν-6sin(ν)-2νcos(ν)2ν3,ϵ(τ42)=1-8(ν-sin(ν))ν3,ϵ(τ43)=1-6(ν-sin(ν))ν3,ϵ(τ44)=1.
The error constant is
(4.7)C4=limν→0ϵ(τ41)2+ϵ(τ42)2+ϵ(τ43)2+ϵ(τ44)2=103.Next we consider the following three-stage method based on Heun's method of order three(4.8)0131323023b1(ν)b2(ν)b3(ν)
Solving the phase-fitting and amplification-fitting conditions and the first-order condition, we obtain
(4.9)b1(ν)=3ν+2ν3+6νcos(ν)-9sin(ν)2ν3,b2(ν)=-3(2ν+νcos(ν)-3sin(ν))ν3,b3(ν)=-9(ν-sin(ν))2ν3.
As ν→0, bi(ν),i=1,2,3, have the following Taylor expansions:
(4.10)b1(ν)=14+780ν2-113360ν4+116128ν6-1926611200ν8+⋯,b2(ν)=-120ν2-1420ν4-120160ν6+11663200ν8+⋯,b3(ν)=34-380ν2+11120ν4-180640ν6+18870400ν8+⋯.
The method given by (4.8) and (4.9) is also of order three and we denote it by FRK3s3bI.The error coefficients for FRK3s3bI are given as follows:(4.11)ϵ(τ41)=1+-40ν+4νcosν+36sin(ν)9ν3,ϵ(τ42)=1-16(ν-sin(ν))3ν3,ϵ(τ43)=1-4(ν-sin(ν))ν3,ϵ(τ44)=1.
The error constant is
(4.12)C4=2239.


#### 4.1.2. Four-Stage Methods of Order Four

Next we consider the following four-stage method based on the so-called 3/8 rule (see the right method in Table  1.2 of Hairer et al. [8])
(4.13)0131323-13111-11b1(ν)b2(ν)b3(ν)b4(ν)Fors=4, the phase-fitting and amplification-fitting conditions (3.18) have the form
(4.14)1-13(b2+2b3+3b4)ν2=cos(ν),b1+b2+(1-13ν2)(b3+b4)=sin(ν)ν.
On the other hand, the first-order and second-order conditions become
(4.15)bTe=b1+b2+b3+b4=1,bTc=13b2+23b3+b4=12.
Solving (4.29) and (4.15), we obtain
(4.16)b1(ν)=-6-9ν2+ν4-6cos(ν)+6νsin(ν)2ν4,b2(ν)=3(2-5ν2+ν4-2cos(ν)+4νsin(ν))2ν4,b3(ν)=3(2+ν2-2cos(ν)-2νsin(ν))2ν4,b4(ν)=3(-2+ν2+2cos(ν))2ν4.
As ν→0, bi(ν) have the following Taylor expansions:
(4.17)b1(ν)=18-7240ν2+34480ν4-11120960ν6+1159667200ν8-11937295360ν10+⋯,b2(ν)=38+13240ν2-1713440ν4+15760ν6-16386688ν8+2929059430400ν10+⋯,b3(ν)=38-148ν2+11920ν4-1134400ν6+114515200ν8-12235340800ν10+⋯,b4(ν)=18-1240ν2+113440ν4-11209600ν6+1159667200ν8-129059430400ν10+⋯.
The method given by (4.13) and (4.16) is also of order four and we denote it by FRK4s4aI.Corresponding to the nine fifth-order rooted treesτ5j, j=1,…,9, the error coefficients are given by
(4.18)ϵ(τ51)=420-180ν2+17ν4-420cos(ν)-30νsin(ν)32ν4,ϵ(τ52)=90-75ν2+8ν4-90cos(ν)+30νsin(ν)8ν4,ϵ(τ53)=90-135ν2+16ν4-90cos(ν)+90νsin(ν)16ν4,ϵ(τ54)=1-45(-2+ν2+2cos(ν))2ν4,ϵ(τ55)=-15ν+ν3+15sin(ν)ν3,ϵ(τ56)=1-15(ν-sin(ν))2ν3,ϵ(τ57)=1-15(-2+ν2+2cos(ν))ν4,ϵ(τ58)=1-45(-2+ν2+2cos(ν))2ν4,ϵ(τ59)=1.
The error constant is
(4.19)C5=2033464.Now we consider the following four-stage method based on the classical Runge-Kutta method (see Hairer et al. [8, Page 138], the left method in Table 1.2)
(4.20)01212120121001b1(ν)b2(ν)b3(ν)b4(ν)The phase-fitting and amplification-fitting conditions (3.18) and the first-order and second-order conditions give
(4.21)b1(ν)=2(-2+ν2-2cos(ν))2ν4,b2(ν)=-4ν+ν3+4sin(ν)ν3,b3(ν)=8-8cos(ν)-4νsin(ν)ν4,b4(ν)=b1(ν).
As ν→0, bi(ν) have the following Taylor expansions:
(4.22)b1(ν)=16-1180ν2+110080ν4-1907200ν6+1119750400ν8-121794572800ν10+⋯,b2(ν)=13+130ν2-11260ν4+190720ν6-19979200ν8+11556755200ν10+⋯,b3(ν)=13-145ν2+11680ν4-1113400ν6+111975040ν8-11816214400ν10+⋯,b4(ν)=16-1180ν2+110080ν4-1907200ν6+1119750400ν8-121794572800ν10+⋯.It is easily verified that the method defined by (4.20) and (4.21) is of order four and we denote it by FRK4s4bI. The method was originally obtained by Simos [30].Corresponding to the nine fifth-order rooted treest5j, j=1,…,9, the error coefficients of FRK4s4bI are given by
(4.23)ϵ(τ51)=280-140ν2+11ν4-280cos(ν)16ν4,ϵ(τ52)=30-20ν2+2ν4-30cos(ν)+5νsin(ν)2ν4,ϵ(τ53)=30-30ν2+4ν4-30cos(ν)+15νsin(ν)4ν4,ϵ(τ54)=1-15(-2+ν2+2cos(ν))ν4,ϵ(τ55)=10-10ν2+ν4-10cos(ν)+5νsin(ν)ν4,ϵ(τ56)=-5ν+ν3+5sin(ν)ν3,ϵ(τ57)=1-10(-2+ν2+2cos(ν))ν4,ϵ(τ58)=ϵ(τ54),ϵ(τ59)=1.
Then the error constant is
(4.24)C5=288148.


### 4.2. FRK Methods of Type II

The coefficients of Type II methods are determined by solving the phase-amplification-fitting conditions and update phase-amplification-fitting conditions. Due to the limitation of the number of parameters, we can only consider methods of at least four stages.For the method given by (4.13), the conditions in Theorem 3.3 assume
(4.25)1-13(b2+2b3+3b4)ν2+13b4ν4=cos(ν),13ν(3b1+3b2+(b3+b4)(3-ν2))=sin(ν),b1(ν)+b2(ν)cos(13ν)+b3(ν)cos(23ν)+b4(ν)cos(ν)=sin(ν)ν,b2(ν)sin(13ν)+b3(ν)sin(23ν)+b4(ν)sin(ν)=1-cos(ν)ν.
Solving this system, we can obtain bi(v), i=1,…,4. Their analytic expressions are complicated. Here we present their Taylor expansions as follows:
(4.26)b1(ν)=18-4997200v2+11210718144000v4-1169505729393280000v6+2827383107116397388800000ν8+𝒪(v10),b2(ν)=38+4433600v2-34702318144000v4+54752513674160000v6-11085686363116397388800000ν8+𝒪(v10),b3(ν)=38-11288v2+227103680v4-6946695878656000v6+145605312116316160000ν8+𝒪(v10),b4(ν)=18-7450v2+2031118144000v4-102308314696640000v6+492479531116397388800000ν8+𝒪(v10).
It can be verified that the method given by (4.13) and (4.26) is of order four. We denote the method by FRK4s4bII.The error coefficients of FRK4s4aII are given by(4.27)ϵ(t51)=964+118123040v2-883929030400v4+𝒪(v6),ϵ(t52)=-332+2411920v2-134411612800v4+𝒪(v6),ϵ(t53)=-1164+4993840v2-803479676800v4+𝒪(v6),ϵ(t54)=-78+730v2-203111209600v4+𝒪(v6),ϵ(t55)=-32+43160v2-5003302400v4+𝒪(v6),ϵ(t56)=-14+43320v2-5003604800v4+𝒪(v6),ϵ(t57)=-14+745v2-203111814400v4+𝒪(v6),ϵ(t58)=ϵ(t54),ϵ(t59)=1.
Then the error constant is
(4.28)C5(ν)=2033464.For the method given by (4.13), the conditions in Theorem 3.3 assume
(4.29)1-13(b2+2b3+3b4)ν2=cos(ν),(b1+b2+b3+b4)ν-(13b3+13b4)ν3+13b4ν4=sin(ν),b1(ν)+b2(ν)cos(12ν)+b3(ν)cos(12ν)+b4(ν)cos(ν)=sin(ν)ν,b2(ν)sin(12ν)+b3(ν)sin(12ν)+b4(ν)sin(ν)=1-cos(ν)ν.
Solving this system, we obtain
(4.30)b1(ν)=4(ν-2sin(ν/2))sin(ν/2)ν2(-4+ν2+4cos(ν/2)),b2(ν)=Lν4(-4+ν2+4cos(ν/2)),b3(ν)=-8(νcos(ν/2)-2sin(ν/2))sin(ν/2)ν4,b4(ν)=b1(ν),
where
(4.31)L=2sin(ν2)(8ν-4ν3+ν5+4ν(-4+ν2)cos(ν2)+8νcos(ν)2sin(ν2)+32sin(ν2)-8ν2sin(ν2)-16sin(ν)+4ν2sin(ν)(ν2)).
It can be verified that the method given by (4.20) and (4.30) is of order four. We denote the method by FRK4s4bII.The error coefficients of FRK4s4bII are given by(4.32)ϵ(t51)=W8ν2Q,ϵ(t52)=1+5sin(ν/2)ν4Q×(ν(-4+ν2)cos(ν2)+2(ν-2ν3+νcos(ν)+(4+3ν2)sin(ν2)-2sin(ν))),ϵ(t53)=1-15sin(ν/2)(-νcos(ν/2)Q+2(ν3-(4+ν2)sin(ν/2)+2sin(ν)))2ν4Q,ϵ(t54)=1-30(ν-2sin(ν/2))sin(ν/2)ν2Q,ϵ(t55)=1+10sin(ν/2)ν4Q×(ν(-4+ν2)cos(ν2)+2(ν-ν3+νcos(ν)+(4+ν2)sin(ν2)-2sin(ν))),ϵ(t56)=1-40cos(ν/4)sin3(ν/4)(ν3+4νcos(ν/2)-8sin(ν/2))ν4Q,ϵ(t57)=1-20(ν-2sin(ν/2))sin(ν/2)ν2Q,ϵ(t58)=1-ϵ(t54)ϵ(t59)=1,
where
(4.33)W=160-32ν2+8ν4+2(-5+16ν2)cos(ν2)-160cos(ν)+10cos(3ν2)-140νsin(ν2)-5ν3sin(ν2),Q=-4+ν2+4cos(ν2).
Then the error constant is
(4.34)C5(ν)=288148.


## 4.1. FRK Methods of Type I

The coefficients of Type I methods are obtained by solving the phase-amplification-fitting conditions and some low-order conditions.


### 4.1.1. Three-Stage Methods of Order Three

We begin by considering the following three-stage method(4.1)012121-12b1(ν)b2(ν)b3(ν)
The phase-fitting and amplification-fitting conditions (3.18) for the scheme (4.1) have the form
(4.2)1-12(b2+2b3)ν2=cos(ν),(b1+b2+b3)-b3ν2=sin(ν)ν.
On the other hand, the first-order condition is given by
(4.3)bTe=b1(ν)+b2(ν)+b3(ν)=1.
Solving (4.2) and (4.3), we obtain
(4.4)b1(ν)=-ν-ν3-2νcos(ν)+sin(ν)ν3,b2(ν)=-2νcos(ν)+2sin(ν)ν3,b3(ν)=ν-sin(ν)ν3.
As ν→0, bi(ν) have the following Taylor expansions
(4.5)b1(ν)=16+340ν2-1350400ν4+17362880ν6-11900800ν8+⋯,b2(ν)=23-115ν2+1420ν4-122680ν6+11995840ν8+⋯,b3(ν)=16-1120ν2+15040ν4-1362880ν6+139916800ν8+⋯.It is easy to verify that the method defined by (4.1) and (4.4) is of order three and we denote it by FRK3s3aI.The error coefficients for FRK3s3aI are given as follows:(4.6)ϵ(τ41)=1-8ν-6sin(ν)-2νcos(ν)2ν3,ϵ(τ42)=1-8(ν-sin(ν))ν3,ϵ(τ43)=1-6(ν-sin(ν))ν3,ϵ(τ44)=1.
The error constant is
(4.7)C4=limν→0ϵ(τ41)2+ϵ(τ42)2+ϵ(τ43)2+ϵ(τ44)2=103.Next we consider the following three-stage method based on Heun's method of order three(4.8)0131323023b1(ν)b2(ν)b3(ν)
Solving the phase-fitting and amplification-fitting conditions and the first-order condition, we obtain
(4.9)b1(ν)=3ν+2ν3+6νcos(ν)-9sin(ν)2ν3,b2(ν)=-3(2ν+νcos(ν)-3sin(ν))ν3,b3(ν)=-9(ν-sin(ν))2ν3.
As ν→0, bi(ν),i=1,2,3, have the following Taylor expansions:
(4.10)b1(ν)=14+780ν2-113360ν4+116128ν6-1926611200ν8+⋯,b2(ν)=-120ν2-1420ν4-120160ν6+11663200ν8+⋯,b3(ν)=34-380ν2+11120ν4-180640ν6+18870400ν8+⋯.
The method given by (4.8) and (4.9) is also of order three and we denote it by FRK3s3bI.The error coefficients for FRK3s3bI are given as follows:(4.11)ϵ(τ41)=1+-40ν+4νcosν+36sin(ν)9ν3,ϵ(τ42)=1-16(ν-sin(ν))3ν3,ϵ(τ43)=1-4(ν-sin(ν))ν3,ϵ(τ44)=1.
The error constant is
(4.12)C4=2239.


### 4.1.2. Four-Stage Methods of Order Four

Next we consider the following four-stage method based on the so-called 3/8 rule (see the right method in Table  1.2 of Hairer et al. [8])
(4.13)0131323-13111-11b1(ν)b2(ν)b3(ν)b4(ν)Fors=4, the phase-fitting and amplification-fitting conditions (3.18) have the form
(4.14)1-13(b2+2b3+3b4)ν2=cos(ν),b1+b2+(1-13ν2)(b3+b4)=sin(ν)ν.
On the other hand, the first-order and second-order conditions become
(4.15)bTe=b1+b2+b3+b4=1,bTc=13b2+23b3+b4=12.
Solving (4.29) and (4.15), we obtain
(4.16)b1(ν)=-6-9ν2+ν4-6cos(ν)+6νsin(ν)2ν4,b2(ν)=3(2-5ν2+ν4-2cos(ν)+4νsin(ν))2ν4,b3(ν)=3(2+ν2-2cos(ν)-2νsin(ν))2ν4,b4(ν)=3(-2+ν2+2cos(ν))2ν4.
As ν→0, bi(ν) have the following Taylor expansions:
(4.17)b1(ν)=18-7240ν2+34480ν4-11120960ν6+1159667200ν8-11937295360ν10+⋯,b2(ν)=38+13240ν2-1713440ν4+15760ν6-16386688ν8+2929059430400ν10+⋯,b3(ν)=38-148ν2+11920ν4-1134400ν6+114515200ν8-12235340800ν10+⋯,b4(ν)=18-1240ν2+113440ν4-11209600ν6+1159667200ν8-129059430400ν10+⋯.
The method given by (4.13) and (4.16) is also of order four and we denote it by FRK4s4aI.Corresponding to the nine fifth-order rooted treesτ5j, j=1,…,9, the error coefficients are given by
(4.18)ϵ(τ51)=420-180ν2+17ν4-420cos(ν)-30νsin(ν)32ν4,ϵ(τ52)=90-75ν2+8ν4-90cos(ν)+30νsin(ν)8ν4,ϵ(τ53)=90-135ν2+16ν4-90cos(ν)+90νsin(ν)16ν4,ϵ(τ54)=1-45(-2+ν2+2cos(ν))2ν4,ϵ(τ55)=-15ν+ν3+15sin(ν)ν3,ϵ(τ56)=1-15(ν-sin(ν))2ν3,ϵ(τ57)=1-15(-2+ν2+2cos(ν))ν4,ϵ(τ58)=1-45(-2+ν2+2cos(ν))2ν4,ϵ(τ59)=1.
The error constant is
(4.19)C5=2033464.Now we consider the following four-stage method based on the classical Runge-Kutta method (see Hairer et al. [8, Page 138], the left method in Table 1.2)
(4.20)01212120121001b1(ν)b2(ν)b3(ν)b4(ν)The phase-fitting and amplification-fitting conditions (3.18) and the first-order and second-order conditions give
(4.21)b1(ν)=2(-2+ν2-2cos(ν))2ν4,b2(ν)=-4ν+ν3+4sin(ν)ν3,b3(ν)=8-8cos(ν)-4νsin(ν)ν4,b4(ν)=b1(ν).
As ν→0, bi(ν) have the following Taylor expansions:
(4.22)b1(ν)=16-1180ν2+110080ν4-1907200ν6+1119750400ν8-121794572800ν10+⋯,b2(ν)=13+130ν2-11260ν4+190720ν6-19979200ν8+11556755200ν10+⋯,b3(ν)=13-145ν2+11680ν4-1113400ν6+111975040ν8-11816214400ν10+⋯,b4(ν)=16-1180ν2+110080ν4-1907200ν6+1119750400ν8-121794572800ν10+⋯.It is easily verified that the method defined by (4.20) and (4.21) is of order four and we denote it by FRK4s4bI. The method was originally obtained by Simos [30].Corresponding to the nine fifth-order rooted treest5j, j=1,…,9, the error coefficients of FRK4s4bI are given by
(4.23)ϵ(τ51)=280-140ν2+11ν4-280cos(ν)16ν4,ϵ(τ52)=30-20ν2+2ν4-30cos(ν)+5νsin(ν)2ν4,ϵ(τ53)=30-30ν2+4ν4-30cos(ν)+15νsin(ν)4ν4,ϵ(τ54)=1-15(-2+ν2+2cos(ν))ν4,ϵ(τ55)=10-10ν2+ν4-10cos(ν)+5νsin(ν)ν4,ϵ(τ56)=-5ν+ν3+5sin(ν)ν3,ϵ(τ57)=1-10(-2+ν2+2cos(ν))ν4,ϵ(τ58)=ϵ(τ54),ϵ(τ59)=1.
Then the error constant is
(4.24)C5=288148.


## 4.1.1. Three-Stage Methods of Order Three

We begin by considering the following three-stage method(4.1)012121-12b1(ν)b2(ν)b3(ν)
The phase-fitting and amplification-fitting conditions (3.18) for the scheme (4.1) have the form
(4.2)1-12(b2+2b3)ν2=cos(ν),(b1+b2+b3)-b3ν2=sin(ν)ν.
On the other hand, the first-order condition is given by
(4.3)bTe=b1(ν)+b2(ν)+b3(ν)=1.
Solving (4.2) and (4.3), we obtain
(4.4)b1(ν)=-ν-ν3-2νcos(ν)+sin(ν)ν3,b2(ν)=-2νcos(ν)+2sin(ν)ν3,b3(ν)=ν-sin(ν)ν3.
As ν→0, bi(ν) have the following Taylor expansions
(4.5)b1(ν)=16+340ν2-1350400ν4+17362880ν6-11900800ν8+⋯,b2(ν)=23-115ν2+1420ν4-122680ν6+11995840ν8+⋯,b3(ν)=16-1120ν2+15040ν4-1362880ν6+139916800ν8+⋯.It is easy to verify that the method defined by (4.1) and (4.4) is of order three and we denote it by FRK3s3aI.The error coefficients for FRK3s3aI are given as follows:(4.6)ϵ(τ41)=1-8ν-6sin(ν)-2νcos(ν)2ν3,ϵ(τ42)=1-8(ν-sin(ν))ν3,ϵ(τ43)=1-6(ν-sin(ν))ν3,ϵ(τ44)=1.
The error constant is
(4.7)C4=limν→0ϵ(τ41)2+ϵ(τ42)2+ϵ(τ43)2+ϵ(τ44)2=103.Next we consider the following three-stage method based on Heun's method of order three(4.8)0131323023b1(ν)b2(ν)b3(ν)
Solving the phase-fitting and amplification-fitting conditions and the first-order condition, we obtain
(4.9)b1(ν)=3ν+2ν3+6νcos(ν)-9sin(ν)2ν3,b2(ν)=-3(2ν+νcos(ν)-3sin(ν))ν3,b3(ν)=-9(ν-sin(ν))2ν3.
As ν→0, bi(ν),i=1,2,3, have the following Taylor expansions:
(4.10)b1(ν)=14+780ν2-113360ν4+116128ν6-1926611200ν8+⋯,b2(ν)=-120ν2-1420ν4-120160ν6+11663200ν8+⋯,b3(ν)=34-380ν2+11120ν4-180640ν6+18870400ν8+⋯.
The method given by (4.8) and (4.9) is also of order three and we denote it by FRK3s3bI.The error coefficients for FRK3s3bI are given as follows:(4.11)ϵ(τ41)=1+-40ν+4νcosν+36sin(ν)9ν3,ϵ(τ42)=1-16(ν-sin(ν))3ν3,ϵ(τ43)=1-4(ν-sin(ν))ν3,ϵ(τ44)=1.
The error constant is
(4.12)C4=2239.


## 4.1.2. Four-Stage Methods of Order Four

Next we consider the following four-stage method based on the so-called 3/8 rule (see the right method in Table  1.2 of Hairer et al. [8])
(4.13)0131323-13111-11b1(ν)b2(ν)b3(ν)b4(ν)Fors=4, the phase-fitting and amplification-fitting conditions (3.18) have the form
(4.14)1-13(b2+2b3+3b4)ν2=cos(ν),b1+b2+(1-13ν2)(b3+b4)=sin(ν)ν.
On the other hand, the first-order and second-order conditions become
(4.15)bTe=b1+b2+b3+b4=1,bTc=13b2+23b3+b4=12.
Solving (4.29) and (4.15), we obtain
(4.16)b1(ν)=-6-9ν2+ν4-6cos(ν)+6νsin(ν)2ν4,b2(ν)=3(2-5ν2+ν4-2cos(ν)+4νsin(ν))2ν4,b3(ν)=3(2+ν2-2cos(ν)-2νsin(ν))2ν4,b4(ν)=3(-2+ν2+2cos(ν))2ν4.
As ν→0, bi(ν) have the following Taylor expansions:
(4.17)b1(ν)=18-7240ν2+34480ν4-11120960ν6+1159667200ν8-11937295360ν10+⋯,b2(ν)=38+13240ν2-1713440ν4+15760ν6-16386688ν8+2929059430400ν10+⋯,b3(ν)=38-148ν2+11920ν4-1134400ν6+114515200ν8-12235340800ν10+⋯,b4(ν)=18-1240ν2+113440ν4-11209600ν6+1159667200ν8-129059430400ν10+⋯.
The method given by (4.13) and (4.16) is also of order four and we denote it by FRK4s4aI.Corresponding to the nine fifth-order rooted treesτ5j, j=1,…,9, the error coefficients are given by
(4.18)ϵ(τ51)=420-180ν2+17ν4-420cos(ν)-30νsin(ν)32ν4,ϵ(τ52)=90-75ν2+8ν4-90cos(ν)+30νsin(ν)8ν4,ϵ(τ53)=90-135ν2+16ν4-90cos(ν)+90νsin(ν)16ν4,ϵ(τ54)=1-45(-2+ν2+2cos(ν))2ν4,ϵ(τ55)=-15ν+ν3+15sin(ν)ν3,ϵ(τ56)=1-15(ν-sin(ν))2ν3,ϵ(τ57)=1-15(-2+ν2+2cos(ν))ν4,ϵ(τ58)=1-45(-2+ν2+2cos(ν))2ν4,ϵ(τ59)=1.
The error constant is
(4.19)C5=2033464.Now we consider the following four-stage method based on the classical Runge-Kutta method (see Hairer et al. [8, Page 138], the left method in Table 1.2)
(4.20)01212120121001b1(ν)b2(ν)b3(ν)b4(ν)The phase-fitting and amplification-fitting conditions (3.18) and the first-order and second-order conditions give
(4.21)b1(ν)=2(-2+ν2-2cos(ν))2ν4,b2(ν)=-4ν+ν3+4sin(ν)ν3,b3(ν)=8-8cos(ν)-4νsin(ν)ν4,b4(ν)=b1(ν).
As ν→0, bi(ν) have the following Taylor expansions:
(4.22)b1(ν)=16-1180ν2+110080ν4-1907200ν6+1119750400ν8-121794572800ν10+⋯,b2(ν)=13+130ν2-11260ν4+190720ν6-19979200ν8+11556755200ν10+⋯,b3(ν)=13-145ν2+11680ν4-1113400ν6+111975040ν8-11816214400ν10+⋯,b4(ν)=16-1180ν2+110080ν4-1907200ν6+1119750400ν8-121794572800ν10+⋯.It is easily verified that the method defined by (4.20) and (4.21) is of order four and we denote it by FRK4s4bI. The method was originally obtained by Simos [30].Corresponding to the nine fifth-order rooted treest5j, j=1,…,9, the error coefficients of FRK4s4bI are given by
(4.23)ϵ(τ51)=280-140ν2+11ν4-280cos(ν)16ν4,ϵ(τ52)=30-20ν2+2ν4-30cos(ν)+5νsin(ν)2ν4,ϵ(τ53)=30-30ν2+4ν4-30cos(ν)+15νsin(ν)4ν4,ϵ(τ54)=1-15(-2+ν2+2cos(ν))ν4,ϵ(τ55)=10-10ν2+ν4-10cos(ν)+5νsin(ν)ν4,ϵ(τ56)=-5ν+ν3+5sin(ν)ν3,ϵ(τ57)=1-10(-2+ν2+2cos(ν))ν4,ϵ(τ58)=ϵ(τ54),ϵ(τ59)=1.
Then the error constant is
(4.24)C5=288148.


## 4.2. FRK Methods of Type II

The coefficients of Type II methods are determined by solving the phase-amplification-fitting conditions and update phase-amplification-fitting conditions. Due to the limitation of the number of parameters, we can only consider methods of at least four stages.For the method given by (4.13), the conditions in Theorem 3.3 assume
(4.25)1-13(b2+2b3+3b4)ν2+13b4ν4=cos(ν),13ν(3b1+3b2+(b3+b4)(3-ν2))=sin(ν),b1(ν)+b2(ν)cos(13ν)+b3(ν)cos(23ν)+b4(ν)cos(ν)=sin(ν)ν,b2(ν)sin(13ν)+b3(ν)sin(23ν)+b4(ν)sin(ν)=1-cos(ν)ν.
Solving this system, we can obtain bi(v), i=1,…,4. Their analytic expressions are complicated. Here we present their Taylor expansions as follows:
(4.26)b1(ν)=18-4997200v2+11210718144000v4-1169505729393280000v6+2827383107116397388800000ν8+𝒪(v10),b2(ν)=38+4433600v2-34702318144000v4+54752513674160000v6-11085686363116397388800000ν8+𝒪(v10),b3(ν)=38-11288v2+227103680v4-6946695878656000v6+145605312116316160000ν8+𝒪(v10),b4(ν)=18-7450v2+2031118144000v4-102308314696640000v6+492479531116397388800000ν8+𝒪(v10).
It can be verified that the method given by (4.13) and (4.26) is of order four. We denote the method by FRK4s4bII.The error coefficients of FRK4s4aII are given by(4.27)ϵ(t51)=964+118123040v2-883929030400v4+𝒪(v6),ϵ(t52)=-332+2411920v2-134411612800v4+𝒪(v6),ϵ(t53)=-1164+4993840v2-803479676800v4+𝒪(v6),ϵ(t54)=-78+730v2-203111209600v4+𝒪(v6),ϵ(t55)=-32+43160v2-5003302400v4+𝒪(v6),ϵ(t56)=-14+43320v2-5003604800v4+𝒪(v6),ϵ(t57)=-14+745v2-203111814400v4+𝒪(v6),ϵ(t58)=ϵ(t54),ϵ(t59)=1.
Then the error constant is
(4.28)C5(ν)=2033464.For the method given by (4.13), the conditions in Theorem 3.3 assume
(4.29)1-13(b2+2b3+3b4)ν2=cos(ν),(b1+b2+b3+b4)ν-(13b3+13b4)ν3+13b4ν4=sin(ν),b1(ν)+b2(ν)cos(12ν)+b3(ν)cos(12ν)+b4(ν)cos(ν)=sin(ν)ν,b2(ν)sin(12ν)+b3(ν)sin(12ν)+b4(ν)sin(ν)=1-cos(ν)ν.
Solving this system, we obtain
(4.30)b1(ν)=4(ν-2sin(ν/2))sin(ν/2)ν2(-4+ν2+4cos(ν/2)),b2(ν)=Lν4(-4+ν2+4cos(ν/2)),b3(ν)=-8(νcos(ν/2)-2sin(ν/2))sin(ν/2)ν4,b4(ν)=b1(ν),
where
(4.31)L=2sin(ν2)(8ν-4ν3+ν5+4ν(-4+ν2)cos(ν2)+8νcos(ν)2sin(ν2)+32sin(ν2)-8ν2sin(ν2)-16sin(ν)+4ν2sin(ν)(ν2)).
It can be verified that the method given by (4.20) and (4.30) is of order four. We denote the method by FRK4s4bII.The error coefficients of FRK4s4bII are given by(4.32)ϵ(t51)=W8ν2Q,ϵ(t52)=1+5sin(ν/2)ν4Q×(ν(-4+ν2)cos(ν2)+2(ν-2ν3+νcos(ν)+(4+3ν2)sin(ν2)-2sin(ν))),ϵ(t53)=1-15sin(ν/2)(-νcos(ν/2)Q+2(ν3-(4+ν2)sin(ν/2)+2sin(ν)))2ν4Q,ϵ(t54)=1-30(ν-2sin(ν/2))sin(ν/2)ν2Q,ϵ(t55)=1+10sin(ν/2)ν4Q×(ν(-4+ν2)cos(ν2)+2(ν-ν3+νcos(ν)+(4+ν2)sin(ν2)-2sin(ν))),ϵ(t56)=1-40cos(ν/4)sin3(ν/4)(ν3+4νcos(ν/2)-8sin(ν/2))ν4Q,ϵ(t57)=1-20(ν-2sin(ν/2))sin(ν/2)ν2Q,ϵ(t58)=1-ϵ(t54)ϵ(t59)=1,
where
(4.33)W=160-32ν2+8ν4+2(-5+16ν2)cos(ν2)-160cos(ν)+10cos(3ν2)-140νsin(ν2)-5ν3sin(ν2),Q=-4+ν2+4cos(ν2).
Then the error constant is
(4.34)C5(ν)=288148.


## 5. Numerical Simulation

In this section, we proceed to solve numerically a two-gene regulatory system without self-regulation (see Widder et al. [2] and Polynikis et al. [3])
(5.1)r˙1=m1H+(p2;θ2,n2)-γ1r1,r˙2=m2H-(p1;θ1,n1)-γ2r2,p˙1=k1r1-δ1p1,p˙2=k2r2-δ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 Gi, ki is the translation rate of Ri, γi is the degradation rate of Ri and δi is the degradation rate of Pi. The two functions
(5.2)H+(p2;θ2,n2)=p2n2p2n2+θ2n2,H-(p1;θ1,n1)=θ1n1p1n1+θ1n1
are the Hill functions of activation and repression, respectively. The parameters n1, n2 are the Hill coefficients, θ1, θ2 are the expression thresholds.The solution(r1*,r2*,p1*,p2*) of the system r˙1=r˙2=p˙1=p˙2=0 is an equilibrium point of the system (5.1). Now we take the values of parameters as follows:
(5.3)n1=n2=3,m1=1.15,m2=2.35,k1=k2=1,γ1=γ2=1,δ1=δ2=1,θ1=θ2=0.21.
For any initial point (r1(0),r2(0),p1(0),p2(0)) near the equilibrium point, the system (5.1) has a stable limit cycle. Figure 1 presents three phase-orbits projected in the protein plane starting at (0.6,0.8,0,0), (0.6,0.8,0.1,1), and (0.8,0.2,0.5,0.2), respectively.Figure 1
Efficiency curves.In order to test the effectiveness of the new FRK methods proposed in this paper, we apply them to the above two-gene regulatory system. For comparison, we also employ several highly efficient integrators from the literature. The numerical methods we choose for experiments are as follows:(i)
FRK3s3aI: the three-stage RK method of order three given by (4.1) and (4.4) in Section 4 of this paper;(ii)
FRK3s3bI: the three-stage RK method of order three given by (4.8) and (4.9) in Section 4 of this paper;(iii)
FRK4s4aI: the four-stage RK method of order four given by (4.13) and (4.16) in Section 4 of this paper;(iv)
FRK4s4bI: the four-stage RK method of order four given by (4.20) and (4.21) in Section 4 of this paper;(v)
FRK4s4aII: the four-stage RK method of order four given by (4.13) and (4.26) in Section 4 of this paper;(vi)
FRK4s4bII: the four-stage RK method of order four given by (4.20) and (4.30) in Section 4 of this paper;(vii)
ARK4: the second four-stage adapted RK method of order four given in Subsection 3.2 of Franco [16];(viii)
EFRK4: the four-stage exponentially fitted RK method of order four given in [31];(ix)
RK4: the classical RK method of order four presented in [8] (the prototype method of FRK4s4bI and FRK4s4bII);(x)
RK38:  the so-called 3/8 rule,  the right method in Table  1.2 of Hairer et al. [8] (the prototype method of FRK4s4aI and FRK4s4aII).


### 5.1. Experiment 1: Protein Error Growth

We first integrate the problem (5.1)–(5.3) on the interval [0,400] with the fixed step size h=1 and the initial values
(5.4)(r10,r20,p10,p20)=(0.6,0.8,0.4,0.6).
In Figure 2, we plot the absolute global errors of Protein 1 produced by the methods FRK3s3aI and FRK3s3bI compared with their prototype RK methods. In Figure 3, we plot the absolute global errors of Protein 1 produced by the methods FRK4s4aI, FRK4s4aII, FRK4s4bI and FRK4s4bII compared with their prototype RK methods.Comparison of (a) RK3s3a and FRK3s3aI (left), (b) RK3s3b and FRK3s3bI (right).
(a)
(b)Comparison of (a) RK3/8, FRK4s4aI and FRK4s4aII (left), (b) RK4, FRK4s4bI and FRK4s4bII (right).
(a)
(b)It can be seen from Figures2 and 3 that the FRK methods are more accurate than their prototype RK methods. Moreover, for appropriate fitting frequency (here ω=1.006), the FRK methods of type II can be more accurate than the corresponding FRK methods of type I and their prototype RK methods.


### 5.2. Experiment 2: Efficiency

Next we compare the computational efficiency of the methods under consideration. In this experiment, we also take the same initial data as Experiment I. The problem is integrated on the interval[0,30] with step sizes h=1/2i,i=1,2,3,4 for the four-stage methods and with 3h/4 for the three-stage methods. The numerical results are presented in Figure 4, where the horizontal axis stands for  the computational effort measured by the number of function evaluations required and the vertical axis stands for the maximal global error (MGE) of Protein 1, both in the digital logarithm scale.Efficiency curves.
(a)
(b)It can be seen from Figure4 that the FRK methods are more efficient than their prototype RK methods and are more efficient than the other frequency depending methods of the same algebraic order.


## 5.1. Experiment 1: Protein Error Growth

We first integrate the problem (5.1)–(5.3) on the interval [0,400] with the fixed step size h=1 and the initial values
(5.4)(r10,r20,p10,p20)=(0.6,0.8,0.4,0.6).
In Figure 2, we plot the absolute global errors of Protein 1 produced by the methods FRK3s3aI and FRK3s3bI compared with their prototype RK methods. In Figure 3, we plot the absolute global errors of Protein 1 produced by the methods FRK4s4aI, FRK4s4aII, FRK4s4bI and FRK4s4bII compared with their prototype RK methods.Comparison of (a) RK3s3a and FRK3s3aI (left), (b) RK3s3b and FRK3s3bI (right).
(a)
(b)Comparison of (a) RK3/8, FRK4s4aI and FRK4s4aII (left), (b) RK4, FRK4s4bI and FRK4s4bII (right).
(a)
(b)It can be seen from Figures2 and 3 that the FRK methods are more accurate than their prototype RK methods. Moreover, for appropriate fitting frequency (here ω=1.006), the FRK methods of type II can be more accurate than the corresponding FRK methods of type I and their prototype RK methods.


## 5.2. Experiment 2: Efficiency

Next we compare the computational efficiency of the methods under consideration. In this experiment, we also take the same initial data as Experiment I. The problem is integrated on the interval[0,30] with step sizes h=1/2i,i=1,2,3,4 for the four-stage methods and with 3h/4 for the three-stage methods. The numerical results are presented in Figure 4, where the horizontal axis stands for  the computational effort measured by the number of function evaluations required and the vertical axis stands for the maximal global error (MGE) of Protein 1, both in the digital logarithm scale.Efficiency curves.
(a)
(b)It can be seen from Figure4 that the FRK methods are more efficient than their prototype RK methods and are more efficient than the other frequency depending methods of the same algebraic order.


## 6. Conclusions and Discussions

In this paper, classical Runge-Kutta methods are adapted to stable limit cycles of first-order dynamical systems of genetic regulation. The newly developed phase-fitted and amplification-fitted Runge-Kutta methods (FRK) adopt functions of the productν=ωh of the fitting frequency ω and the step size h as weight coefficients in the update. FRK methods have zero dispersion and zero dissipation when applied to the standard linear oscillator y′=iωy. That is to say, they can preserve initial phase and amplification as time extends. Therefore, FRK integrators are a kind of integrators which preserve the oscillation structure of the problem.We note that as the fitting frequency tends to zero, FRK methods reduce to their classical prototype methods. Furthermore, an FRK method has the same algebraic order and the same error constant with its prototypes method. Numerical experiments illustrate the high efficiency of FRK methods compared with their prototype methods and some other frequency depending methods like exponentially fitted RK methods.Theorem3.3 gives a pair of sufficient conditions for modified RK type methods to be phase-fitted and amplification-fitted. The coefficients of the FRK methods of type I (FRK3s3aI, FRK3s3bI, FRK4s4aI, and FRK4s4b) are obtained with these conditions combined with an appropriate number of low order conditions. Since the FRK methods of type II (FRK4s4aII and FRK4s4bII) have updates that are also phase-amplification-fitted, they can be more accurate than those methods whose updates do not have this property.In practical computations of oscillatory problems, the true frequency is, in general, not available. The fitting frequency contained in an FRK method is just an estimate of the true frequency. For techniques of estimating principal frequencies we refer to the papers [32–37]. In this experiment, we have a new discovery that for a specific FRK method, the choice of the value of the fitting frequency ω affects its efficiency to some extent. For instance, in Section 5, we have discovered that ω=1.1966 is the best choice for FRK4s4aII, while ω=1.006 is the best for FRK4s4bII. This interesting discovery constitutes a challenge to the traditional viewpoint on the choice of fitting frequency for frequency depending methods for which a “best” estimate of frequency can be obtain from the problem under consideration. The results of Experiment 2 of Section 5 show that this best frequency also depends on the method itself.



---

*Source: 673296-2012-12-18.xml*