Computational and Mathematical Methods in Medicine

Volume 2017, Article ID 2729683, 16 pages

https://doi.org/10.1155/2017/2729683

## Steady-State-Preserving Simulation of Genetic Regulatory Systems

^{1}College of Sciences, Nanjing Agricultural University, Nanjing 210095, China^{2}College of Horticulture, Nanjing Agricultural University, Nanjing 210095, China^{3}Department of Mathematics, University of Lagos, Lagos 23401, Nigeria

Correspondence should be addressed to Xiong You; nc.ude.uajn@xuoy

Received 11 September 2016; Revised 11 December 2016; Accepted 18 December 2016; Published 19 January 2017

Academic Editor: Irini Doytchinova

Copyright © 2017 Ruqiang Zhang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A novel family of exponential Runge-Kutta (expRK) methods are designed incorporating the stable steady-state structure of genetic regulatory systems. A natural and convenient approach to constructing new expRK methods on the base of traditional RK methods is provided. In the numerical integration of the one-gene, two-gene, and p53-mdm2 regulatory systems, the new expRK methods are shown to be more accurate than their prototype RK methods. Moreover, for nonstiff genetic regulatory systems, the expRK methods are more efficient than some traditional exponential RK integrators in the scientific literature.

#### 1. Introduction

One of the challenges in systems biology is to understand how biochemical molecules, such as DNAs, mRNAs, and proteins, interact to form harmonic and uniform cellular systems which give rise to life (see [1, 2]). The synthetic genetic regulatory networks (GRNs) play an important role in the investigation of protein regulation processes in living organisms (see [3–6]). By introducing ordinary differential equations (ODEs) to describe the rates of change in the concentrations of biochemical molecules, such as mRNAs and proteins, more detailed understanding and insights of the dynamic behavior exhibited by biological systems can be achieved (see [7]). The first attempt to model the oscillation in genetic regulation in terms of ODEs was made by Goodwin [8]. A standard presentation of general regulatory dynamics can be found in the monographs by Thomas and D’Ari [9] and by Fall et al. [10]. Iwamoto et al. [11] presented a dynamical model of the DNA damage signaling pathway that includes p53 and whole cell cycle regulation and explored the relationship between p53 oscillation and cell fate selection. ODEs models admit mathematically qualitative and quantitative analysis to reveal the profound properties from steady states with stability, bistability, oscillation, and limit cycles to chaos (see [12–17] and the references therein).

A typical system of ODEs governing an -gene activation-inhibition system has the form (see Polynikis et al. [15]) where, for , is the concentration of mRNA produced by gene , is the concentration of protein translated from mRNA , is the degradation rate of , and is the degradation rate of . Function is the translation function. Function is the regulation function, typically taking the form of a sum of products of functions . If protein has no effect on gene , does not appear in . The partial derivative if protein is an* activator* of gene and if protein is an* inhibitor* of gene . The genetic regulatory system (1) can be written in matrix formwhere and are -dimensional vectors representing the concentrations of mRNAs and proteins at time , respectively, and , , , and are diagonal matrices.

The analytical solution of system (1) is in general not acquirable due to the nonlinearity of the functions and . Therefore, in order to explore the dynamics of the gene regulatory system (1), one usually resorts to numerical solution. For example, Shinto et al. [18] proposed a kinetic simulation model of metabolic pathways that describes the dynamic behaviors of metabolites in acetone-butanol-ethanol (ABE) production by* Clostridium saccharoperbutylacetonicum* N1-4 using a simulator WinBEST-KIT. So far, differential equations for genetic regulation are mostly solved by the classical four-stage Runge-Kutta (RK) method or by the Runge-Kutta-Fehlberg adaptive method (see Hairer et al. [19]). However, general-purpose RK methods have not taken into account the special structure of system (1) and fail to capture the dynamical features of the system effectively, especially in the long time simulation (see Figure 11). Thanks to new advances in the last two decades, new approaches have been developed aiming at preserving the intrinsic geometric or physical structures of the true solution. A comprehensive account of structure-preserving algorithms can be found in the monographs by Stuart and Humphries [20], Hairer et al. [21], and Wu et al. [22]. Recently Hochbruck and Ostermann [23] investigated exponential Runge-Kutta methods for initial value problems of parabolic differential equations. This type of methods simulates exactly the linear structure of the differential equations. Defterli et al. [24] and Weber et al. [25] considered discretizing and optimizing the so-called gene-environment networks based on usually finite data series.

From the dynamics point of view, there are two basic categories of genetic regulatory systems: Category 1 consists of systems having sustained oscillation, such as limit cycles; Category 2 consists of systems having steady states. For the genetic regulatory system (1) with a limit-cycle structure, You [26] proposed a new class of phase-fitted and amplification-fitted Runge-Kutta type methods which were shown to be more effective and more efficient than the traditional Runge-Kutta methods of the same order. Very recently, You et al. [27] developed a splitting approach for genetic regulatory systems with a stable steady state. In the numerical simulation, the new splitting methods constructed in that paper are shown to be remarkably more effective and more suitable for long-term computation with large steps than the general-purpose Runge-Kutta methods. In order to respect the oscillatory feature of the solution of some genetic regulatory systems, Chen et al. [28] developed a new type of exponentially fitted TDRK (EFTDRK) methods. Zhang et al. [29] constructed a family of phase-fitted symmetric splitting methods of order two and order four. The result of the numerical experiment on the Lotka–Volterra system shows that the new phase-fitted symmetric splitting methods are more effective than their prototype splitting methods and can preserve the invariant of the system in the long term compared with the classical Runge-Kutta method of order four.

The purpose of this paper is to develop a novel type of exponential RK methods for the simulation of genetic regulatory systems which have an asymptotically stable steady state. In Section 2 we present the general scheme of exponential Runge-Kutta (expRK) methods for solving initial value problems of ODEs based on a matrix form of the variation-of-constants formula. A convenient approach of transiting traditional RK methods into a special type of expRK methods is given. In Section 3 we integrate the above three regulatory systems by the new expRK methods as well as their prototype RK methods for comparison. Section 4 is devoted to conclusive remarks. In Appendix, we analyze the linear stability and phase properties of the expRK methods.

#### 2. Exponential Runge-Kutta Methods

##### 2.1. Formulation of Exponential RK Methods for Systems with a Stable Steady-State Structure

Prior to dealing with the genetic regulatory system (2) numerically, we first consider the general initial value problem (IVP) of the autonomous system of ODEs where and “” represents the derivative of with respect to time. We make the following assumptions on system (3):(**A1**)The origin is a steady state of the system; that is, .(**A2**)The Jacobian has eigenvalues of negative real parts in a neighborhood of the origin. Then, according to the theorem in Section 8.5 of Hirsch et al. [30], the origin is an asymptotically stable steady state; that is, for every solution of system (3) through a point in the neighborhood of the steady state, .(**A3**)The function is continuously differentiable and satisfies the Lipschitz condition; that is, there exists a constant (called the* Lipschitz constant*) such that for all .

Let us recall the general-purpose Runge-Kutta methods for IVP (3).

*Definition 1. *An -stage Runge-Kutta (RK) method for system (3) has the schemewhere , are real numbers.

Scheme (5) can be expressed briefly by the Butcher tableau of its coefficientsThe order conditions for the RK method (5) can be found in Hairer et al. [19]. Note that the general RK scheme (5) does not take into account the special structure of the equilibrium structure of the system so that the computational results are usually not satisfactory. In order to simulate more effectively system (3) with a steady state at the origin, we rewrite problem (3) in an equivalent formwhere the matrix , the function , is the Euclidean norm, and is called the* rate matrix* of system (3). The matrix form of the variation-of-constants formula for system (3) is given bywhere and are real numbers and . Approximating the integral on the right-hand side of (8) by some effective quadrature formulas leads to a numerical integrator. In general we have the following definition of the so-called exponential Runge-Kutta methods.

*Definition 2. *An -stage exponential Runge-Kutta method for system (3) has the schemewhere , are real numbers and and , are real matrix-valued functions of matrix .

It is convenient to express scheme (9) by the Butcher tableauwhere . For a comprehensive review of exponential integrators with the construction, analysis of convergence and error bounds, order conditions, example integrators, and applications, the reader is referred to Hochbruck and Ostermann [23].

It is noted that if we need to integrate a system near a stable steady state , the exponential RK scheme (9) should take the formwhere with .

##### 2.2. A Special Class of Exponential RK Methods

Based on an RK method (5), we can formulate, as a special case of the exponential RK method (9), the following scheme for system (3): where , , and , are the constant coefficients of the RK method (5). Note that as (), scheme (12) reduces to the RK method (5). The latter is called the* prototype RK method* of the former. In the sequel of this paper, scheme (12) will be referred to as an* expRK method*. Its Butcher tableau can be written as follows:where , , and .

In Kronecker’s notation, scheme (12) can be written as where is the unit matrix and .

Among simple examples are the following:(a)The* exponential Euler method* (explicit), denoted by expEuler:(b)The* exponential backward Euler method* (implicit), denoted by expImEuler:(c)The* exponential trapezoidal rule* (implicit), denoted by expTrap:(d)The* exponential Heun rule* (explicit), denoted by expHeun:(e)The* exponential midpoint rule* (implicit), denoted by expMid:

Two typical expRK methods, denoted by expRK3/8 and expRK4, have the prototype RK methods of order four, denoted by RK3/8 and RK4, respectively, whose respective Butcher tableaux are given in Page 138 of [19]

The (algebraic) order is a measure of the accuracy of numerical method. A method is said to have* order * if its local error .

Theorem 3. *The expRK method (12) has the same algebraic order as its prototype RK method.*

*Proof. *The conclusion follows by expanding the exponential functions , , and in scheme (12) in series of , hence in series of , and comparing this series with that of the true solution.

The next theorem asserts that exponential RK method (12) preserves the steady state of system (3).

Theorem 4. *Suppose that in system (3) satisfies the Lipschitz condition and the origin is a steady state of the system; that is, . Then the origin is also a fixed point of the expRK method (12) for small step size .*

The proof is given in Appendix A.

In order to apply the expRK method (9) or (12) to system (2), we first use a coordinate transform to translate the steady state of the system to the origin and yieldswhere is the Jacobian matrix of at point and . Then system (2) can be written in the form (7) with where the rate matrix and the function .

#### 3. Numerical Illustrations

From Theorem 4, contrast to traditional RK methods, exponential RK methods, especially expRK methods, retain the rate of growth, phase, and amplification of the exact solution of the test equation (B.1) without error. Then it is reasonable to expect expRK (12) to be more effective than their prototype RK methods. In this section, in order to compare their effectiveness, we apply them to three test systems—one-gene self-regulation system, two-gene cross-regulation system, and the p53-mdm2 system.

*(I) A One-Gene System of Self-Regulation*. The first model we consider is the one-gene system with self-regulation given bywhere represents the action of an inhibitory protein that acts as a dimer and are positive constants. For a similar model with delays see Xiao and Cao [31].

*(II) A Two-Gene System with Cross-Regulation*. The second model is a two-gene activation-inhibition system with cross-regulation (studied by Polynikis et al. [15], Widder et al. [17], Chen et al. [32], You [26], and You et al. [27])where, for , is the concentration of mRNA produced by gene , is the concentration of protein , is the maximal transcription rate of gene , is the translation rate of mRNA , and and are the degradation rates of mRNA and protein , respectively. The functionsare the* Hill functions* of activation and repression, respectively. The parameters are the* expression thresholds*. The integer value of , called the* Hill coefficient*, stands for the number of protein monomers required for saturation of binding to DNA. It is easy to see that the activation function is increasing in and the repression function is decreasing in .

*(III) The p53-mdm2 System*. The third model is for the damped oscillation of the p53-mdm2 regulatory pathway. Strictly speaking, this system is not of the form (1). We adopt this model since its solutions also have a stable steady-state structure of interest. The system, given by van Leeuwen et al. [33] with the small transient stress stimulus , has the formwhere represents the concentration of the p53 tumor suppressor, (mdm2) is the concentration of the p53’s main negative regulator, is the concentration of the p53-mdm2 complex, is the concentration of an active form of p53 that is resistant against mdm2-mediated degradation, are de novo synthesis rates, are production rates, are reverse reactions (e.g., dephosphorylation), is the degradation rate of active p53, and is the saturation coefficient.

##### 3.1. Accuracy Test

###### 3.1.1. The One-Gene System

Steady states of system (24) can be determined by the cubic equation and the relation . If the system is written in the form (7), the rate matrixand the function where . With the parameter values (provided by [31])this system has a unique positive steady state where the rate matrix has eigenvalues , where is the imaginary unit satisfying . Since the two eigenvalues both have negative real parts, the steady state is asymptotically stable. Figure 1(a) shows three solution trajectories on the phase plane starting at , respectively. Figure 1(b) shows the time evolution of concentrations of mRNA and protein starting at .