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.

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 y′ = iωy (i 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.

2. Models

2.1. Cross-Regulation Systems

An N-gene regulatory system can be modeled by a system of ordinary differential equations,or in matrix form,where m(t) = (m 1(t),…, m (t)) and p(t) = (p 1(t),…, p (t)) are N-dimensional vectors which represent the concentrations of mRNAs and proteins at time t, respectively, a dot “·” over a variable represents time derivative, R(p(t)) = (R 1(p(t)),…, R (p(t))), T(m(t)) = (T 1(m 1(t)),…, T (m (t))), Γ = diag⁡(γ 1,…, γ ), and M = diag⁡(μ 1,…, μ ) are diagonal matrices. For i = 1,…, N, m is the concentration of mRNA R produced by gene G , p is the concentration of protein P produced by R , γ is the degradation rate of R , and μ is the degradation rate of P . Function T (m ) is the translation function. Function R (p(t)) is the regulation function, typically defined as sums of functions of p 1,…, p . If ∂R /∂p > 0, protein P is an activator of gene G . If ∂R /∂p < 0, protein P is an inhibitor of gene G . If protein P has no effect on gene G , p does not appear in R . Usually, T(m(t)) = (κ 1 m 1(t),…, κ m (t)) = diag⁡(κ 1,…, κ )m(t).

In particular, we will be concerned with the following two-gene system:where m 1 and m 2 are the concentrations of R 1 and R 2, respectively, p 1 and p 2 are the concentrations of the corresponding P 1 and P 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 R 1 and R 2, respectively, μ 1 and μ 2 are the degradation rates of P 1 and P 2, respectively,are the Hill functions for activation and inhibition, respectively, n 1 and n 2 are the Hill coefficients, and θ 1 and θ 2 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 n represents the degree of cooperativity, and K 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 k . 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 M, P 0, P 1, P 2, and P are the concentrations of cytosolic per mRNA, unphosphorylated PER, monophosphorylated PER, bisphosphorylated PER, and nuclear PER, respectively, V and K (i = 1,…, 4) denote the maximum rate and Michaelis constant of the kinase(s) and phosphatase(s) involved in the reversible phosphorylation of P 0 into P 1 and of P 1 into P 2, 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 y : [0, +∞) → ℝ , “” represents the first derivative of y with respect to time, and f : ℝ → ℝ 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):

System (6) has a steady state y ; that is, f(y ) = 0.

System (6) has oscillatory solution near y ; that is, the Jacobian J = f′(y ) has at least a pair of complex eigenvalues with nonzero imaginary part.

A special form of two-derivative Runge-Kutta (TDRK) method readswhere a , b , c , i, j = 1,…, s are real numbers and g(y)≔f′(y)f(y). 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) methods where a , i, j = 1,…, s are constants and the coefficients η(ν), β(ν), b (ν), i = 1,…, s are the functions of ν = hω with h 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 aswhereand I is d × d 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 c 2 = (0, c 2 2,…, c 2),  c 3 = (0, c 2 3,…, c 3), c 4 = (0, c 2 4,…, c 4), 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 C 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 ℒ satisfies for the function y(t) in the reference set where λ , k = 0,1,…, N are real or complex numbers and N, M are nonnegative integers.

In this subsection, we consider four-stage explicit modified TDRK methods given by the following tableau:The coefficients η(ν), β(ν), b (ν), i = 1,2, 3,4 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 s = 4 with higher order terms neglected giveWe can solve system (20)-(21) for η(ν), β(ν), and b (ν), i = 1,2, 3,4, 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 b (ν), i = 1,2, 3,4.

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 b (ν), i = 1,2, 3,4. For small values of |ν|, the following Taylor series should be used: 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:

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: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.

RK6: the classical RK method of order six presented in [8].

TDRK4s6: the classical TDRK method of order six presented in [18].

EFTDRK4s6a, EFTDRK4s6b, EFTDRK4s6c: the three four-stage EFTDRK methods of order six derived in Section 3 of this paper.

ETFRK4: the fourth-order exponentially and trigonometrically fitted RK method constructed by Simos [22].

EFRK4: the fourth-order exponentially fitted RK method constructed by Vanden Berghe et al. [23].

MRK4: the fourth-order modified RK method constructed by Van de Vyver [24].

4.1. The Two-Gene System

Denote y = (m 1, m 2, p 1, p 2). 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 ω = 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.

Figure 1

Time evolution of proteins in the two-gene system.

For the initial values (m 1(0), m 2(0), p 1(0), p 2(0)) = (0.6,0.8,0.4,0.6) near the equilibrium point, we solve system (3) on the interval [0,100] by the methods EFTDRK4s6a, EFTDRK4s6b, and EFTDRK4s6c with step sizes h = 1/2, i = 2,3, 4,5, with respect to best fitting frequencies ω. Tables 1 –4 display the global error (GE) of Protein 1 for comparison.

Table 1

Two-gene systems: global error of Protein 1 with step size h = 1/4.

	ETFRK4	EFRK4	MRK4	RK6
ω	0.659	0.100	0.412
GE	3.9971e − 08	2.1795e − 07	1.3827e − 04	1.3216e − 06
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		0	0	0
GE	2.0703e − 07	2.0703e − 07	2.0703e − 07	2.0703e − 07

Table 2

Two-gene systems: global error of Protein 1 with step size h = 1/8.

	ETFRK4	EFRK4	MRK4	RK6
ω	0.653	1.815	0.640
GE	1.0441e − 09	1.3170e − 06	6.7673e − 06	2.9632e − 08
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		3.145	2.395	1.573
GE	4.0028e − 09	3.6251e − 12	3.1940e − 09	4.4587e − 12

Table 3

Two-gene systems: global error of Protein 1 with step size h = 1/16.

	ETFRK4	EFRK4	MRK4	RK6
ω	0.751	0.022	0.813
GE	1.6283e − 10	1.6144e − 08	2.7404e − 07	5.3117e − 10
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		2.964	2.292	1.482
GE	6.6183e − 11	5.6954e − 14	7.6938e − 14	5.6954e − 14

Table 4

Two-gene systems: global error of Protein 1 with step size h = 1/32.

	ETFRK4	EFRK4	MRK4	RK6
ω	1.356	0.011	0.931
GE	3.7377e − 12	2.9958e − 09	9.2907e − 09	8.7911e − 12
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		2.891	2.185	1.445
GE	1.0523e − 12	1.2212e − 15	1.6653e − 15	1.8874e − 15

In Figure 2 we compare the efficiency of the eight methods by plotting the global error against the number of evaluations of nonlinear functions f and g.

Figure 2

Two-gene system: global error versus function evaluations.

4.2. PER Oscillations in Drosophila

Denote y = (M, P 0, P 1, P 2, P ). The Jacobian of system (5) is given byand the function f′(y)f(y) = g(y) = (g 1, g 2, g 3, g 4, g 5) 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 ξ 3,4 have nonzero imaginary parts, the solution near the steady state (M , P 0 , P 1 , P 2 , P ) is oscillatory with frequency ω = 0.265113. Figure 3 displays time evolution of the concentration of the nuclear protein P .

Figure 3

Time evolution of nuclear PER in Drosophila.

For the initial values (M(0), P 0(0), P 1(0), P 2(0), P (0)) = (0.1,0.25,0.25,0.25,0.25) near the equilibrium point, we integrate system (3) on the interval [0,100] with step sizes h = 1/2, i = 1,2, 3,4, 5. The numerical results are presented in Tables 5 –9.

Table 5

PER oscillations in Drosophila: global error of P with step size h = 1/2.

	ETFRK4	EFRK4	MRK4	RK6
ω	0.811	1.311	0.290
GE	2.2113e − 04	0.3821	0.0011	0.0011
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		1.971	1.369	0.991
GE	1.2580e − 06	8.4100e − 10	1.1408e − 09	1.3104e − 09

Table 6

PER oscillations in Drosophila: global error of P with step size h = 1/4.

	ETFRK4	EFRK4	MRK4	RK6
ω	0.010	0.002	0
GE	4.1924e − 10	2.4214e − 09	3.0138e − 08	7.4781e − 08
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		1.934	1.363	0.968
GE	1.3468e − 08	1.0491e − 11	2.7490e − 11	2.0292e − 11

Table 7

PER oscillations in Drosophila: global error of P with step size h = 1/8.

	ETFRK4	EFRK4	MRK4	RK6
ω	0	0	0.047
GE	3.7713e − 09	3.7713e − 09	5.5061e − 11	3.5174e − 10
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		1.895	1.356	0.947
GE	1.5941e − 10	1.5565e − 13	1.9101e − 13	4.9882e − 13

Table 8

PER oscillations in Drosophila: global error of P with step size h = 1/16.

	ETFRK4	EFRK4	MRK4	RK6
ω	0	0	0.057
GE	3.0163e − 10	3.0163e − 10	4.2553e − 12	2.8773e − 12
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		1.871	1.354	0.935
GE	2.1183e − 12	3.8858e − 16	3.1641e − 15	6.3283e − 15

Table 9

PER oscillations in Drosophila: global error of P with step size h = 1/32.

	ETFRK4	EFRK4	MRK4	RK6
ω	0	0	0.061
GE	2.0380e − 11	2.0380e − 11	2.3392e − 13	3.1641e − 14
	TDRK4s6	EFTDRK4s6a	EFTDRK4s6b	EFTDRK4s6c
ω		1.880	1.370	0.931
GE	3.1086e − 14	3.3307e − 16	0	1.1102e − 16

In Figure 4 we plot the efficiency curves for the eight methods.

Figure 4

Two-gene system: global error versus function evaluations.

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 ν = ωh, the product of the fitting frequency ω and the step size h, 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 f), the global error (GE) of an EFTDRK method (e.g., EFTDRK4s6a) depends on the coefficients (and thus ν = ωh). 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⁡(iωt) and exp⁡(2iωt). 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.

Table 10

Trees of order up to order five and the values of corresponding functions.

Tree	ρ	α	γ	Φi	ℱ
	1	1	1		f
	2	1	2	1	g
	3	1	6	c i	g′f
	4	1	12	c i 2	g′′(f, f)
	4	1	24	∑ja_ij	g y g
	5	1	20	c i 3	g (3)(f, f, f)
	5	3	40	c_i∑ja_ij	g′′(f, g)
	5	1	120	∑ja_ijc_j	g′g

Table 11

Trees of order six and the values of corresponding functions.

Tree	ρ	α	γ	Φi	ℱ
	6	1	30	c i 4	g (4)(f, f, f, f)
	6	6	60	ci2∑ja_ij	g (3)(f, f, g)
	6	3	120	∑ja_ij^2	g′′(g, g)
	6	4	180	c_i∑ja_ijc_j	g′′(f, g′f)
	6	1	360	∑ja_ijcj2	g′g′′(f, f)
	6	1	720	∑j,ka_ija_jk	g′g′′g