Research Article  Open Access
Exponentially Fitted TwoDerivative RungeKutta Methods for Simulation of Oscillatory Genetic Regulatory Systems
Abstract
Oscillation is one of the most important phenomena in the chemical reaction systems in living cells. The general purpose simulation algorithms fail to take into account this special character and produce unsatisfying results. In order to enhance the accuracy of the integrator, the secondorder derivative is incorporated in the scheme. The oscillatory feature of the solution is captured by the integrators with an exponential fitting property. Three practical exponentially fitted TDRK (EFTDRK) methods are derived. To test the effectiveness of the new EFTDRK methods, the twogene system with crossregulation and the circadian oscillation of the period protein in Drosophila are simulated. Each EFTDRK method has the best fitting frequency which minimizes the global error. The numerical results show that the new EFTDRK methods are more accurate and more efficient than their prototype TDRK methods or RK methods of the same order and the traditional exponentially fitted RK method in the literature.
1. Introduction
The qualitative analysis and quantitative simulation of gene expression and regulation play an important role in understanding the dynamics of complex processes in cells. Ordinary differential equations (ODEs) have proved to be one of the powerful tools for modeling the complex dynamics of genetic regulation in cells, 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], and Gérard and Goldbeter [5] and the references therein).
Due to the nonlinearity of the ODE models, the closed form of solution is usually not acquirable. Therefore, in order to reveal the dynamics of such gene regulatory systems, one usually resorts to numerical simulation. Up till now, differential equations of gene regulatory systems are mostly simulated by RungeKutta (RK) methods, especially by the classical fourthorder RungeKutta method, or by the RungeKuttaFehlberg adaptive method (see Butcher [6, 7] and Hairer et al. [8]).
As is often observed in experiments, in a variety of cell processes, genes exhibit an oscillatory behavior. Among examples are sustained oscillations associated with circadian clocks, enzyme synthesis, or the cell cycle (see Goldbeter [9] and Jolley et al. [10]). Unfortunately, when applied to these oscillatory systems, the general purpose RK method often fails to produce satisfactorily efficient numerical results since it did not take into account the special structure of the true solution. 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 and (ii) the true dynamical behavior of the system cannot be preserved as expected in longterm integration.
Recently, some authors have proposed to adapt traditional integrators to problems whose solutions are oscillatory or periodic (see Bettis [11], Gautschi [12], Martín and Ferándiz [13], and Raptis and Simos [14]). Bettis [15] constructed a threestage method and a fourstage method which can solve the equation exactly. Very recently You [16] developed a new family of phasefitted and amplification fitted methods of RK type which have been proved very effective for genetic regulatory systems with a limitcycle structure. You et al. [17] considered a splitting approach for the numerical simulation of genetic regulatory networks with a stable steady state structure. The numerical results of the simulation of a onegene network, a twogene network, and a p53mdm2 network showed that the new splitting methods constructed in this paper are remarkably more effective and more suitable for longterm computation with large steps than the traditional general purpose RungeKutta methods.
Motivated by the work of Chan and Tsai [18] and Fang et al. [19] on the twoderivative RungeKutta methods (TDRK), the objective of this paper is to develop a novel type of exponentially fitted twoderivative RungeKutta (EFTDRK) methods for simulating genetic regulatory systems with an oscillatory structure. 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 longterm integration of gene regulatory systems. In Section 2 we present the main models: one is for gene systems with cross regulations and the other is for the circadian oscillation of the period protein in Drosophila. In Section 3, three EFTDRK methods of algebraic order six are constructed. In Section 4 the new EFTDRK methods are applied to the simulation of the two genetic regulatory systems given in Section 2 and their efficiency is compared with that of a sixthorder traditional RK and a sixthorder TDRK method and three exponentially fitted RK methods. Section 5 is devoted to some conclusive remarks and discussions. The mathematical theory of order conditions and the evaluation of best fitting frequencies for EFTDRK methods are presented in the Appendix.
2. Models
2.1. CrossRegulation Systems
An gene regulatory system can be modeled by a system of ordinary differential equations,or in matrix form,where and are dimensional vectors which represent the concentrations of mRNAs and proteins at time , respectively, a dot “” over a variable represents time derivative, , , , and are diagonal matrices. For , is the concentration of mRNA produced by gene , is the concentration of protein produced by , is the degradation rate of , and is the degradation rate of . Function is the translation function. Function is the regulation function, typically defined as sums of functions of . If , protein is an activator of gene . If , protein is an inhibitor of gene . If protein has no effect on gene , does not appear in . Usually,
In particular, we will be concerned with the following twogene system:where and are the concentrations of and , respectively, and are the concentrations of the corresponding and , respectively, and are the maximal transcription rates of gene and gene , respectively, and are the degradation rates of and , respectively, and are the degradation rates of and , respectively,are the Hill functions for activation and inhibition, respectively, and are the Hill coefficients, and and are the thresholds.
2.2. Circadian Rhythms
Another model we are interested in is for circadian oscillations in the period protein (PER) in Drosophila. A crucial mechanism for oscillations in the model is the negative feedback exerted by nuclear PER on the production of per mRNA. This negative feedback will be described by a Hill type equation, where the Hill coefficient represents the degree of cooperativity, and represents the threshold inhibition constant. It is assumed that per mRNA is synthesized in the nucleus and immediately transfers to the cytosol, where it accumulates at a maximum rate ; there it is degraded enzymatically, in a Michaelian manner, at a maximum rate . The rate of synthesis of PER is characterized by an apparent firstorder rate constant . PER experiences a series of phosphorylations (Edery et al. [20]). For simplicity, it is assumed that there are three states of the protein: unphosphorylated, monophosphorylated, and bisphosphorylated. Goldbeter [21] formulated the fivevariable system of equations as follows:where the dependent variables , , , , and are the concentrations of cytosolic per mRNA, unphosphorylated PER, monophosphorylated PER, bisphosphorylated PER, and nuclear PER, respectively, and denote the maximum rate and Michaelis constant of the kinase(s) and phosphatase(s) involved in the reversible phosphorylation of into and of into , respectively.
3. Methods
3.1. Modified TwoDerivative RungeKutta Methods
We begin by considering the general initial value problem (IVP) of ordinary differential equationswhere , “” represents the first derivative of with respect to time, and is a sufficiently smooth function. Based on experimental observation of oscillatory behavior in genetic regulatory systems, it is reasonable to make the following assumptions on system (6):(i)System (6) has a steady state ; that is, .(ii)System (6) has oscillatory solution near ; that is, the Jacobian has at least a pair of complex eigenvalues with nonzero imaginary part.
A special form of twoderivative RungeKutta (TDRK) method readswhere , , , are real numbers and The order conditions for TDRK methods can be found in Chan and Tsai [18].
In applications, for some choice of the parameter values, system (2) has oscillatory solutions. This motivates us to consider the modified twoderivative RungeKutta (TDRK) methodswhere , are constants and the coefficients , , , are the functions of with being the step size and being the principal frequency of the problem. We assume that so that as the modified TDRK method (8) reduces to a traditional special TDRK method called the limit method or the prototype method of (8).
In Kronecker’s notation, scheme (8) can be written compactly aswhereand is identity matrix. The TDRK method (7) can be briefly expressed by the following Butcher tableau of coefficients:
The order conditions for modified TDRK methods will be derivative via the theory of biordered trees in Appendix. For purpose of construction of practical methods, we list the sixthorder conditions as follows:where , , , and a dot “” between two vectors indicates a pairwise multiplication. If we assume that then conditions (12) become
3.2. Exponentially Fitted TwoDerivative RungeKutta Methods
With the modified TDRK method (8) we associate a linear operator on , the linear space of functions on with continuous second derivatives, defined by
Definition 1. The modified TDRK method (8) is said to be exponentially fitted if the linear operator satisfies for the function in the reference set where are real or complex numbers and are nonnegative integers.
In this subsection, we consider fourstage explicit modified TDRK methods given by the following tableau:The coefficients , , , will be obtained by the exponential fitting conditions for some specific reference sets.
3.2.1. First EFTDRK Method
We take the reference setand assume that the linear operator in (15) vanishes for all functions in This leads toThe third and fourth conditions in (14) for with higher order terms neglected giveWe can solve system (20)(21) for , , and , whose expressions are extremely complicated. For small values of , we have their Taylor series used as follows: It is easily verified that these coefficients satisfy all conditions in (14). Therefore the method has algebraic order six and we denote this method as EFTDRK4s6a.
3.2.2. Second EFTDRK Method
We take the reference set and assume that the linear operator (15) vanishes for all functions in . Then we obtain the expression of , , and
For small values of , the following Taylor series should be used:
It is easily verified that these coefficients satisfy all conditions in (14). Therefore the method has algebraic order six and we denote this method as EFTDRK4s6b.
3.2.3. Third EFTDRK Method
We take the reference setand we assume that the linear operator (15) vanishes for all functions in . Then we obtain the expression of , , and For small values of , the following Taylor series should be used: We denote this method as EFTDRK4s6c.
It is noted that as , the new methods EFTDRK4s6a, EFTDRK4s6b, and EFTDRK4s6c reduce to a traditional TDRK method given by the following tableau:
4. Results
In order to examine the effectiveness of the EFTDRK methods proposed in this paper, we apply these methods as well as a sixthorder RK method and a sixthorder TDRK method to the two genetic regulatory systems presented in Section 2. The numerical methods we will use are listed as follows:(i)RK6: the classical RK method of order six presented in [8].(ii)TDRK4s6: the classical TDRK method of order six presented in [18].(iii)EFTDRK4s6a, EFTDRK4s6b, EFTDRK4s6c: the three fourstage EFTDRK methods of order six derived in Section 3 of this paper.(iv)ETFRK4: the fourthorder exponentially and trigonometrically fitted RK method constructed by Simos [22].(v)EFRK4: the fourthorder exponentially fitted RK method constructed by Vanden Berghe et al. [23].(vi)MRK4: the fourthorder modified RK method constructed by Van de Vyver [24].We will compare the efficiency of these methods by plotting the decimal logarithm of the maximal global error against the computational effort measured by the number of function evaluations.
4.1. The TwoGene System
Denote . The Jacobian of system (3) is given by and the functionWe take the values of parameters as follows (Polynikis et al. [3]): We solve system by Newton iteration for the unique positive steady state of system (3) which is given byThe Jacobian matrix at the steady state has the eigenvalues Since these eigenvalues have nonzero imaginary parts, the solution near the steady state is oscillatory with frequency . This oscillatory behavior of the two proteins is shown in Figure 1 which is plotted straightly by the classical RungeKutta method of order four.
For the initial values near the equilibrium point, we solve system (3) on the interval by the methods EFTDRK4s6a, EFTDRK4s6b, and EFTDRK4s6c with step sizes , , with respect to best fitting frequencies . Tables 1–4 display the global error (GE) of Protein 1 for comparison.




In Figure 2 we compare the efficiency of the eight methods by plotting the global error against the number of evaluations of nonlinear functions and .
4.2. PER Oscillations in Drosophila
Denote . The Jacobian of system (5) is given byand the function withwhere The parameter values in model (5) are given by (Goldbeter [21])
System is solved for the unique positive steady state Thus the five eigenvalues of the Jacobian matrix at the steady state areSince have nonzero imaginary parts, the solution near the steady state is oscillatory with frequency . Figure 3 displays time evolution of the concentration of the nuclear protein .
For the initial values near the equilibrium point, we integrate system (3) on the interval with step sizes The numerical results are presented in Tables 5–9.

