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 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.
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 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 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.
2. Models
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. [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: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 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 [18].
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).
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 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:
4. 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.
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. [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 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
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.
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.
5. 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 as a function , 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 is (much) larger than . That is, EFTDRK4s6a may be less effective than its prototype TDRK method (TDRK4s6, the case ) 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 and . 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.
Appendix
A. Bicolored Trees and Order Conditions for Modified TDRK Methods
A.1. Bicolored Trees
In this appendix we present the theory of bicolored rooted trees which forms the basis of the derivation of order conditions for modified TDRK methods.
Definition A.1. The set of bicolored rooted trees is defined recursively as follows:(i)The graph and belong to ; they are denoted by and , respectively. The black vertex in each of the two trees is called the root.(ii)If , then the tree obtained by grafting the roots of onto the white vertex of .
Definition A.2. For each tree , the elementary differential is a function defined recursively as follows:(i) and ().(ii)For ,
Definition A.3. The function is defined recursively as follows:(i) and () .(ii)For , For each , is called the order of the tree . Actually, is the number of vertices of the tree . The set of trees of order is denoted by .
Definition A.4. The function is defined recursively as follows:(i) and () .(ii)For , where is the multiplicity of
Definition A.5. The function is defined recursively as follows:(i) and () .(ii)For , For each , is called the density of the tree .
Definition A.6. For the modified TDRK method (8), the elementary weight vector () is defined as follows: (i)() .(ii)For , , ,
Theorem A.7. The exact solution and the numerical solution 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 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. [8]:where and the order , the density , the vector of elementary weights , the elementary differential at , and the function are defined in Appendix A.1. The modified TDRK method (8) has (algebraic) order if its local error satisfies .
The above analysis leads to the order conditions as stated in the following theorem.
Theorem A.8. The modified TDRK method (8) has (algebraic) order if and only if its coefficients satisfy the following conditions:
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 is the set of trees of order . Therefore for small values of step size the global error over the grids can be approximated as This shows that the global error depends on and the values of and the elementary differentials 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors are grateful to Professor Albert Goldbeter from Free University of Brussels for his invaluable comments and suggestions which help much to improve the paper. This work was supported by the Natural Science Foundation of China (NSFC) under Grant no. 11171155 and the Fundamental Research Funds for the Central Universities under Grant no. KYZ201424.