基于果蝇周期蛋白PER调控网络构建的两导数Runge-Kutta算法

Two-derivative Runge-Kutta algorithms for configuration of the regulatory network of period protein PER in Drosophila

为更准确地模拟基因调控网络的动力学行为,考虑采用新的两导数Runge-Kutta模拟算法。这种数值算法的特点是在格式中融入了真实解的二阶结构,使得模拟结果更为精确。与通用的经典四级四阶Runge-Kutta算法RK4相比,二级四阶两导数Runge-Kutta算法TDRK2s4具有3个优势:1)构造简单,易于编程实现;2)计算成本小,RK4每步需要计算4次非线性函数,而TDRK2s4只需要计算3次;3)存储空间小,RK4方法需要4个内级存储单元,而TDRK2s4只要2个。TDRK2s4用于模拟具有持续振荡特征的果蝇周期蛋白PER与per mRNA负反馈调控网络。与RK4相比较,对不同步长(h=1/2,1/4,1/8,1/16),TDRK2s4的误差远小于RK4的误差;而且TDRK2s4对大步长(h=1/2)的误差远小于RK4对小步长(h=1/16)的误差。对固定步长(h=1/4),随着模拟区间的延长,TDRK2s4的累积误差增长十分缓慢,而RK4的累积误差增长很快。对不同步长(h=1/2,1/4,1/8,1/16),TDRK2s4计算耗费的CPU时间也明显少于RK4。TDRK2s4的模拟效率远高于RK4。结论:当基因调控网络微分方程组右端函数的导数容易得到时,TDRK算法比传统的RK算法更适用于大步长、长时间的高效数值模拟。

The dynamical mechanism of genetic regulatory networks are usually simulated predominately by the traditional RungeKutta（RK） methods. The purpose of this paper is to pursue new simulation approaches which are more precise and more effective. For the numerical integration of the initial value problems z＇ = g（x,z）,z（x0） = z0,two-derivative Runge-Kutta（TDRK） algorithms of the form,z1= z0＋hg（x0,z0） ＋h2∑s i = 1bili,li= φ（x0＋cih,z0＋cihg（x0,z0） ＋h2∑s j = 1aijlj）（i=1,2,…,s） are considered. These numerical algorithms are characterized by incorporating the second-order structure of the true solution into the scheme which makes the scheme more accurate. Although the two-stage two-derivative Runge-Kutta algorithm（TDRK2s4） has the same algebraic order（four） and the same stability function R（w） = 1 ＋ w ＋ w2/2 ＋ w3/6 ＋ w4/24（that is the function in the iteration z1= R（w） z0,w=hλ,obtained by applying the algorithm to the linear test equation z＇ = λz with step size h） as the general-purposed classical four-stage Runge-Kutta algorithm（RK4）,it possesses three advantages： 1） It is simpler in construction and easier to implement. 2） It requires less computational cost and less storage space-TDRK2s4 involves three evaluations of nonlinear functions g（x0,z0）,l1= φ（x0,z0） and l2= φ（x0＋c2h,z0＋c2hg（x0,z0） ＋h2a21l1） at each step while RK4 involves four k1= g（z0）,k2= g（x0＋c2h,z0＋ha21k1）,k3= g（x0＋c3h,z0＋ha32k2） and k4= g（x0＋ c4h,z0＋ ha43k3）. 3） It requires less storage space-TDRK2s4 needs only two units for storing the internal stages l1and l2while RK4 involves three for k1,k2,k3and k4. TDRK2s4 is applied to the simulation of the negative-feedback regulatory network of period protein（PER） and per mRNA in Drosophila which is characterized by sustained oscillation. RK4 is also used for comparison. The numerical results show that the errors produced by TDRK2s4 are much less than those produced by RK4 for different step sizes h = 1/2,1/4,1/8,1/16,and the errors of TDRK2s4 for a large step size h = 1/2 are much less that those of RK4 for a small step size h = 1/16. For a fixed step size h = 1/4,the accumulated error of TDRK2s4 grows much slowly as the simulation time interval extends while that of RK4 grows rather rapidly. For different step sizes h = 1/2,1/4,1/8,1/16,TDRK2s4 consumes explicitly less CPU time than RK4 does. TDRK2s4 has a much higher efficiency of simulation than RK4. When the derivative of the right-hand side function of the differential equation of a genetic regulatory network is attainable,TDRK algorithms are more suitable for long-time highly effective numerical simulation with large step sizes than traditional RK algorithms are.