Research Article | Open Access
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.
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 , Widder et al. , Polynikis et al. , Altinok et al. , and Gérard and Goldbeter  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. ).
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  and Jolley et al. ). 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 , Gautschi , Martín and Ferándiz , and Raptis and Simos ). Bettis  constructed a three-stage method and a four-stage method which can solve the equation exactly. Very recently You  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.  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  and Fang et al.  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.
2.1. Cross-Regulation 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 two-gene 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 first-order rate constant . PER experiences a series of phosphorylations (Edery et al. ). For simplicity, it is assumed that there are three states of the protein: unphosphorylated, monophosphorylated, and bisphosphorylated. Goldbeter  formulated the five-variable 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.1. Modified Two-Derivative Runge-Kutta 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 two-derivative Runge-Kutta (TDRK) method readswhere , , , are real numbers and The order conditions for TDRK methods can be found in Chan and Tsai .
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) 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).
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:where , , , and a dot “” between two vectors indicates a pairwise multiplication. If we assume that then conditions (12) become
3.2. Exponentially Fitted Two-Derivative Runge-Kutta 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 four-stage 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:
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 .(ii)TDRK4s6: the classical TDRK method of order six presented in .(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 .(v)EFRK4: the fourth-order exponentially fitted RK method constructed by Vanden Berghe et al. .(vi)MRK4: the fourth-order modified RK method constructed by Van de Vyver .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 Two-Gene System
Denote . The Jacobian of system (3) is given by and the functionWe take the values of parameters as follows (Polynikis et al. ): 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 Runge-Kutta 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
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 .