Exponentially Fitted Two-Derivative Runge-Kutta Methods for Simulation of Oscillatory Genetic Regulatory Systems

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 second-order 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 two-gene system with cross-regulation 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.


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 Runge-Kutta (RK) methods, especially by the classical fourth-order Runge-Kutta method, or by the Runge-Kutta-Fehlberg 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 long-term 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 three-stage method and a four-stage method which can solve the equation = ( 2 = −1) exactly. Very recently You [16] developed a new family of phase-fitted and amplification fitted methods of RK type which have been proved very effective for genetic regulatory systems with a limit-cycle 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 one-gene network, a two-gene network, and a p53-mdm2 network showed that the new splitting methods constructed in this paper are remarkably more effective and more suitable for long-term computation with large steps than the traditional general purpose Runge-Kutta methods.
Motivated by the work of Chan and Tsai [18] and Fang et al. [19] on the two-derivative Runge-Kutta methods (TDRK), the objective of this paper is to develop a novel type of exponentially fitted two-derivative Runge-Kutta (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 long-term 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 sixth-order traditional RK and a sixth-order 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.
In particular, we will be concerned with the following two-gene system: where 1 and 2 are the concentrations of 1 and 2 , respectively, 1 and 2 are the concentrations of the corresponding 1 and 2 , respectively, 1 and 2 are the maximal transcription rates of gene 1 and gene 2, respectively, 1 and 2 are the degradation rates of 1 and 2 , respectively, 1 and 2 are the degradation rates of 1 and 2 , respectively, are the Hill functions for activation and inhibition, respectively, 1 and 2 are the Hill coefficients, and 1 and 2 are the thresholds.

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 V ; there it is degraded enzymatically, in a Michaelian manner, at a maximum rate V . The rate of synthesis of PER is characterized by an apparent first-order 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 five-variable system of equations as follows: Computational and Mathematical Methods in Medicine where the dependent variables , 0 , 1 , 2 , and are the concentrations of cytosolic per mRNA, unphosphorylated PER, monophosphorylated PER, bisphosphorylated PER, and nuclear PER, respectively, and ( = 1, . . . , 4) denote the maximum rate and Michaelis constant of the kinase(s) and phosphatase(s) involved in the reversible phosphorylation of 0 into 1 and of 1 into 2 , respectively.

Modified Two-Derivative Runge-Kutta Methods.
We begin by considering the general initial value problem (IVP) of ordinary differential equationṡ = ( ) , where : [0,+∞) → R , "" represents the first derivative of with respect to time, and : R → R 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, ( * ) = 0.
(ii) System (6) has oscillatory solution near * ; that is, the Jacobian = ( * ) has at least a pair of complex eigenvalues with nonzero imaginary part.
In applications, for some choice of the parameter values, system (2) has oscillatory solutions. This motivates us to consider the modified two-derivative Runge-Kutta (TDRK) methods where , , = 1, . . . , are constants and the coefficients (]), (]), (]), = 1, . . . , are the functions of ] = ℎ with ℎ being the step size and being the principal frequency of the problem. We assume that (0) = (0) = 1 so that as → 0 (] → 0) 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 as where = (1, 1, . . . , 1) , and 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 sixth-order conditions as follows: , and a dot "⋅" between two vectors indicates a pairwise multiplication. If we assume that then conditions (12) become

Exponentially Fitted Two-Derivative Runge-Kutta Methods.
With the modified TDRK method (8) we associate a linear operator L on 2 [0, ∞), the linear space of functions on [0, ∞) with continuous second derivatives, defined by Definition 1. The modified TDRK method (8) is said to be exponentially fitted if the linear operator L satisfies for the function ( ) in the reference set where , = 0, 1, . . . , are real or complex numbers and , are nonnegative integers.
In this subsection, we consider four-stage explicit modified TDRK methods given by the following tableau: The coefficients (]), (]), (]), = 1, 2, 3, 4 will be obtained by the exponential fitting conditions for some specific reference sets.

First EFTDRK Method.
We take the reference set and assume that the linear operator L in (15) vanishes for all functions in F . This leads to Computational and Mathematical Methods in Medicine The third and fourth conditions in (14) for = 4 with higher order terms neglected give We can solve system (20) 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.

Second EFTDRK Method.
We take the reference set and assume that the linear operator (15) 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.

Third EFTDRK Method.
We take the reference set and we assume that the linear operator (15) We denote this method as EFTDRK4s6c.
It is noted that as ] → 0, the new methods EFTDRK4s6a, EFTDRK4s6b, and EFTDRK4s6c reduce to a traditional TDRK method given by the following tableau:

Results
In order to examine the effectiveness of the EFTDRK methods proposed in this paper, we apply these methods as well as a sixth-order RK method and a sixth-order 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 four-stage EFTDRK methods of order six derived in Section 3 of this paper. (iv) ETFRK4: the fourth-order exponentially and trigonometrically fitted RK method constructed by Simos [22]. (v) EFRK4: the fourth-order exponentially fitted RK method constructed by Vanden Berghe et al. [23]. (vi) MRK4: the fourth-order 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.
Since these eigenvalues have nonzero imaginary parts, the solution near the steady state is oscillatory with frequency = 0.989478. This oscillatory behavior of the two proteins is shown in Figure 1 which is plotted straightly by the classical Runge-Kutta method of order four.
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 .
In Figure 4 we plot the efficiency curves for the eight methods.    It can be seen from Tables 1-9 and Figures 2 and 4 that the EFTDRK methods are more efficient than the other methods used for comparison.

Conclusions and Discussions
Oscillation is frequently observed in all living cells. Whether or not this feature is accurately preserved in simulation has a critical effect on the comprehension of genetic regulation systems. Now that the traditional simulation algorithms perform poorly in simulating the oscillatory genetic regulatory networks, new and more effective simulation technology is called for. In this paper, classical two-derivative Runge-Kutta methods are adapted to the oscillatory character of genetic regulatory systems. The newly developed exponentially fitted two-derivative Runge-Kutta methods (EFTDRK) adopt functions of ] = ℎ, the product of the fitting frequency and the step size ℎ, as weight coefficients in the update. As the fitting frequency tends to zero, EFTDRK methods reduce to their prototype TDRK methods of the same algebraic order. It should be noted that, in applying an exponentially/trigonometrically fitted method to oscillatory problems, a fitting frequency , an accurate estimate of the principal frequency, must be obtained in advance. However, for a given oscillatory system, the true frequency is in general not available. In the existing literature, all the methods of fitted type share a common value of the fitting frequency once it is well estimated. According to the argument in Appendix A.3, for a given differential equation (given function ), the global error (GE) of an EFTDRK method (e.g., EFTDRK4s6a) depends on the coefficients (and thus ] = ℎ). If we consider GE as a function , GE( ) then we take the fitting frequency as the value of that minimizes the global error. We call it the "best fitting frequency." It is also noted that it may happen that GE( ) is (much) larger than GE(0). That is, EFTDRK4s6a may be less effective than its prototype TDRK method (TDRK4s6, the case = 0) if an inappropriate value of the fitting frequency is employed.
EFTDRK methods have two advantages: they respect the second-order structure of the true solution and they can simulate exactly some standard oscillatory functions such as exp( ) and exp (2 ). These characteristic properties contribute to their high accuracy and high efficiency. The EFTDRK methods developed in this paper, a category of structure-preserving algorithms, open a new approach to simulating genetic regulatory systems with an oscillation structure.
Definition A.5. The function : BT → N is defined recursively as follows: Theorem A.7. The exact solution ( + ℎ) and the numerical solution +1 produced by the modified TDRK method (8) have the following expansions: A.2. Order Conditions. The algebraic order conditions for the TDRK method with constant coefficients presented in [18] are not applicable for the modified TDRK method (8) whose coefficients depend on ] = ℎ. We expand the exact solution ( +ℎ) and the numerical solution +1 in powers of ℎ under the local assumption = ( ): where and and their partial derivatives are evaluated at . The elementary differentials in the above expressions will be represented geometrically by a set of bicolored (rooted) trees, which are a variation of the rooted trees presented in Butcher [6] and Hairer et al. [8]. For the autonomous system (6), the bicolored trees are made up of two types of vertices: black and white. The first and the simplest tree, which represents , consists of a single black vertex. The time derivative of is expressed by a black vertex connecting upward by a branch to a white vertex . A partial derivative of with respect to is expressed by a white vertex connecting by a branch to a black vertex. Tables 10 and 11 give all the trees of order up to five with the values of the related functions defined on them.
The local truncation error of the modified TDRK method (8) can be expressed in terms of rooted trees which are similar to those defined in Hairer et al.    Proof. The result follows from comparing the expansions of (A.7) and independence of elementary differentials.
A. 3. The Choice of the Fitting Frequency. It has been known that the choice of the fitting frequency is crucial to the effectiveness of exponentially or trigonometrically fitted Runge-Kutta methods when applied to initial value problems with oscillatory solutions. In the existing literature, there have been several approaches for estimating the frequency. See, for example, Ixaru and Vanden Berghe [25], Vanden Berghe et al. [26], Ixaru et al. [27], Van de Vyver [28], and Ramos and Vigo-Aguiar [29].
In this paper, we evaluate the "best fitting frequency" for an EFTDRK method by minimizing the global error. According to (A.8) in Appendix A.2, we can write the local truncation error of th order modified TDRK method (8) as follows: where +1 is the set of trees of order +1. Therefore for small values of step size ℎ the global error over the grids , = 0, 1, . . . , can be approximated as This shows that the global error depends on ] = ℎ and the values of ( ) and the elementary differentials F( ) at each and thus depends on the step size ℎ, the fitting frequency of the method, the length of the integration interval, and the function (i.e., the problem to be solved). The "best fitting frequency" can be taken as the minimizer of the previous approximation of GE.