基因调控网络的数值方法_刘雪萍.docx

《基因调控网络的数值方法_刘雪萍.docx》由会员分享，可在线阅读，更多相关《基因调控网络的数值方法_刘雪萍.docx（61页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 NUMERICAL METHODS FOR GENETIC REGULATORY NETWORKS By Xueping LIU Supervised by Professor Xiong YOU A Thesis Submitted To Nanjing Agriculture University In Partial Fulfillment of The Requirements For The Master Degree May, 2014 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独立进行研究工作 所取得的成果。除文中己经注明引用的内容外，本论文不包含任何其他个人或集体

2、已 经发表或撰写过的作品成果。对本文的研究作出重要贡献的个人和集体，均已在文中 以明确方式标明。本人完全意识到本声明的法律结果由本人承担。 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并 向国家有关部门或机构送交论文的复印件和电子版，允许论文被査阅和借阅。本人授 权南京农业大学可以将本学位论文的全部或部分内容编入有关数据库进行检索，可以 采用影印、缩印或扫描等复制手段保存和汇编学位论文。 保密 在 _ 年解密后适用本授权书。本学位论文属于不保密口。 (请在以上方框内打 “ V” ） 基因调控网络的数值方法 摘要 现代分子生物学与细胞生物学的一个核

3、心问题是揭示基因的表达与调控机制 .对 基因调控网络的探究有助于我们更清晰地了解细胞生物化学反应过程，同时对诸如生 物制药和临床医学实践有重要的指导意义 .基因调控网络微分方程的数值模拟，主要是 应用Rung-KTitta(RK)方法 .但是对于许多具有特别结构的基因调控网络，在大步长长 时间积分时，用 Rung-Kutta (RK)方法处理后所得到的结果并不理想，甚至完全失效 . 本文针对某些基因调控网络的特征结构，发展保结构的高效数值算法 . 本文共分四章： 第一章给出了本文的研究 背景 介绍了基因调控的基本原理、基因调控网络 （ ge* netic regulatory networks

4、, GRNs)的数学模型，特别是常微分方程模型以及基本研究方 法 . 第二章针对具有全局稳定平衡点的基因调控网络发展了基于 RungeKutta (RK) 方法的分裂法 .我们将基因调控微分方程的向量场在平衡点附近分解为线性主部和高 阶非线性部分，应用一类特殊的 LieTrotter 分裂法：线性部分精痛积分，非线性部分 用经典的 Runge-Kutta(RK)方法积分 .得到了四个实用的分裂法 Split(exact: Euler), Split(exact: Heun)， Split(exact: RK4)和 Split(exact: RK3/8).单基因自调控网络、两基 因交叉调控网络以

5、及抑癌蛋白 p53-Mdm2 调控网络数值模拟的结果证明，分裂法的计 算精度与计算效率明显高于 Euler, Heun, RK4 和 RK3/8 方法，并且更适用于基因调控 网络的大步长长时间积分 .经分析，分裂法比相应的 RK 方法有更大的稳定性区域 . 第三章考虑两导数 RungeKutta(TDRK)方法 这类方法能够精确积分真实解的 一阶与二阶结构 .通过对两类模拟蛋白质磷酸化过程的调控模型（果蝇周期蛋白模型 与哺乳动物生物钟模型）的数值解法，证明二级四阶 TDRK 方法计算精度和计算效率 明显高于经典的四级四阶 RK 方法 . 第四章考虑时滞基因调控网络 .通过对果绳周期蛋白调控模型

6、和哺乳动物生物钟 模型的数值模拟探讨时滞对基因调控网络动力学行为的影响 . 关键词：基因调控网络，分裂法，两导数 Runge-Kutta 方法，时滞 ABSTRACT NUMERICAL METHODS FOR GENETIC REGULATORY NETWORKS ABSTRACT The exploration on gene regulation is one of the key problems in molecular biology and cell biology. The investigation of genetic regulatory networks allows

7、us to understand the cell mechanisms more clearly, and also has important theoretical and practical significance in biopharmaceuticals and clinical medicine. The principal methods for the simulation of genetic regulatory networks (GRNs) have been Rung-Kutta (RK) methods, which often fail to produce

8、satisfactory results or even lose effect in long-term integration with large step sizes. The purpose of this paper is to develop structure-preserving and highly effective algorithms for the simulation of genetic regulatory networks based on some of their characteristic structures. This paper is divi

9、ded into four chapters: Chapter 1 gives the background of this paper and presents basic principles of genetic regulations, mathematical models of GRNs, especially ordinary differential equation models. We also discussed fundamental research methods. In chapter 2, splitting methods based on Runge-Kut

10、ta (RK) methods are developed to solve the gene regulation systems with a structure of globally stable equlibrium. The vector field of the ODEs of GRNs into two parts, the linear principal part and the higher order nonlinear part. A special type Lie-Trotter splitting method is obtained by using the

11、exact flow for the linear principal part and the numerical flow of a classical Runge-Kutta method for the nonlinear part. Four practical splitting methods, Split (exact: Euler), Split(exact: Heun), Split(exact: RK4) and Split(exact: RK3/8) are derived. When applied to the one- gene self-regulatory n

12、etwork, the two-gene cross-regulatory network and tumor-repressing proteins p53-Mdm2 regulatory network, splitting methods are shown to be more accurate and more efficient than Euler, Heun, RK4 and RK3/8 methods. Moreover, splitting methods are more suitable for the long-time integration with large

13、step sizes. By analysis， a splitting method has a larger stability region than the corresponding RK method. In chapter 3, we consider two-derivative Runge-Kutta (TDRK) methods. These methods have the property that they integrate exactly the first order and second order structures of the true solutio

14、n. When applied to solve the two models with phosphorylation, Drosophila iii 基因调控网络的数值方法 period protein and mammalian circadian clock, the two-stage fourth-order method TDRK2s4 is shown to be more accurate and more efficient than the classical four-stage fourth-order method RK4. In Chapter 4, we con

15、sider genetic regulatory networks with delays. By means of simulating the period protein PER in Drosophila and the circadian clock in mammals, the effect of delays on the dynamic behavior of genetic regulatory networks is explored. Keywords: gene regulatory network; splitting method; two-derivative

16、Runge-Kutta method; delay IV 目录 目录 摘要 i ABSTRACT iii 第一章基因调控网络模型及其动力学方程 1 1.1 基因调控的基本原理 . 1 1.2 基因调控网络与数学模型 . 1 1.3 基因调控网络的常微分方程模型 . 3 1.4 本文的结构 . 5 第二章基因调控网络的分裂法 7 2.1 基于 Rung-Kutta (RK)方法的分裂法 . 7 2.2 基因调控的 LieTVotter 分裂法 . 9 2.3 分裂法用于单基因自调控网络的计算 . 10 2.4 分裂法用于两基因交叉调控网络的计算 . 13 2.5 分裂法用于 p53-Mdm2 调

17、控网络 . 17 2.6 RK 方法和分裂法的稳定性分析 . 22 第三章含有磷酸化过程的基因调控网络的两导数 RimgeKutta 方法 25 3_1 模型 . 25 3.1.1 果媚周期蛋白模型 . 25 3.1.2 哺乳动物生物钟模型 . 28 3_2 两导数 RungeKutta 方法 . 31 3.3 数值模拟 . 32 3.3.1 果蝇周期蛋白模型的计算 . 32 3.3.2 哺乳动物生物钟模型的计算 . 36 第四章时滞基因调控网络 39 V _ 基因调控网络的数值方法 _ 4.1 时滞对果蝇周期基因动力学行为的影响 . 39 4.2 哺乳动物生物钟模型动力学行为的影响 . 42

18、 结论与展望 45 参考文献 47 致谢 53 VI 第一章基因调控网络模型及其动力学方程 第一章基因调控网络模型及其动力学方程 1.1 基因调控的基本原理 英国科学家 Crick 和美国科学家 Watson 1953 年共同提出了 DNA 的反向平行双 螺旋结构，使得遗传学的研究深入到分子水平 .1954 年 Crick 又提出了遗传信息传递 的 “ 中心法则 ” ，即遗传信息只能由核酸向蛋白质传递，而不能相反 .1961 年法国科学家 Mcmod和 Jacob 提出的操纵子学说打开了人类认识基因调控的窗口 .随后，人们开始 系统研究原核生物以及真核生物的基因表达的过程来揭示基因调控的规律

19、. 基因表达是指 DNA 经过转录生成信使 RNA (message RNA， mRNA)， mRNA 再翻 译生成蛋白质的过程 .基因表达的过程始于 RNA 聚合酶和启动子 （ promoter, DNA 中 控制基因表达开始的 DNA 片段）的结合 ,产生启动子一聚合酶复合物，在 RNA 聚合酶 的作用下，开始转录过程 .当将整条基因链转录完毕，阻抑子 （ repressor)将被激活，阻 止转录的进行 .转录结束后就进入更为复杂的翻译过程，此过程以 mRNA 为模板，通 过转移 RNA (transfer RNA， tRNA)实现 .图 1.1 给出了基因的三种状态： “ 裸 ” DNA

20、 通常只有低水平的转录活动，称基因处于基本状态 (basal state); (b)只有当 “ 激活子 ” 与 基因的调控位区结合时，基因的表达才能达到正常水平 ,称基因处于激发状态 （ active state); (c)只要 “ 阻抑子，与其调控位区结合，无论是否有激活子存在，基因都处于非激 活状态 （ inactive state). 1.2 基因调控网络与数学模型 所谓基因调控网络 （ genetic regulatory network)实质上是由一组相互作用的基因、 蛋白质、 RNA 和小分子构成的系统 .一 般说来，基因转录翻译过程的调控可以在三个 水平上发生： DNA 水平、

21、RNA7K 平和蛋白质水平 .其中在 RNA 水平上的调控最为重 要，所以通常所说的基因调控是指转录水平的调控 . 例如，图 1.2给出了一个两基因交叉调控网络 (见 Widder 等 55).假设 DNA 上只 有两个基因 A 和 G2，它们的产物蛋白质分别是 Pi 和 P2. P2对 Gi 的表达起促进作 用， 对G2的表达起抑制作用 .该调控机制可以下式表示： Fi(p2,n2) Gl + H2P2 - Gl ( F*2)n2， F2(Pi； ni) 11 G2 + niPX 一 - - - G2 ( Pl)ni) 基因调控网络的数值方法 图 1.1 基因的三种状态 2 第一章基因调控网

22、络模型及其动力学方程 图 1.2 两基因交叉调控模型 图 1.3 Bool 模型图示 其中巧，朽为 Hill 函数，负和內分别是 Pi 和 P2的平衡浓度， rn， n2为 Hill 系数 . 数学模型已经成为基因调控网络研究的一个十分有效的工具 .常见的基因调控网 络模型大致分为两类，即离散模型 (如 Bool 模型）和连续性模型 (如微分方程模型 ).其 他模型如加权矩阵模型等 .Bool 模型可表述为：在一个有 ri 个节点，每个节点有来自 m 个其他节点的输入边 ,每个节点的状态只有 “ 开 ” 和 “ 关 ”（ on and off),状态 “ 幵 ” 表 示基因转录表达，形成基因产

23、物，状态 “ 关 ” 表示基因未转录 .图 1.3 给出了一个简单的 二基因 Bool 模型的网络拓扑，其中 Gi, G2, G3 分别代表基因， Gi 对 G2 起促进作用， Gi 和 G3相互促进， G2对 G3起阻抑作用 .表 1.1 给出了各基因的真值 . 1.3 基因调控网络的常微分方程模型 常微分方程模型表述 RNA、 蛋白质及其它分子的浓度随时间变化的规律 .1965 年 Goodwin 丨 22首先建立了一个基因调控振荡的微分方程模型 .随后，关于基因调控的数 学分析引起了广泛而持续的兴趣 (见 Tiwari 等 50, Tyson 和 Othmer 51， Smith 46)

24、. 基因调控网络的数值方法 4 第一章基因调控网络模型及其动力学方程 Widder 等 55研究了直到 Hill 系数 n = 4 的两基因交叉调控网络，给出了一维 鞍-结点分叉与二维 Hopf 分叉值的解析表达式 .Polynikis 等 42对两基因网络的不同 ODE模型做了比较，用分段线性函数 (piecewise4inear)近似代替 Hill 函数，并建立了完 全非线性模型 (complete nonlinear model (CNM)，简化非线性模型 (simplified nonlinear model (SNM)，完全分段线性模型 （ complete piecewise li

25、near model (CPWLM),简化 分段模型 （ simplified piecewise linear model (SPWLM)和离散时间模型（ discretetime model),并分析了各类模型的稳定性和分支 Middleton 等 38建立了 Aux/IAA 的负 反馈环数学模型，研究了 Aux/IAA 蛋白质和 mRNA 转换率的比值，得到如下结论 :当 比率充分小时，系统有一个稳定的极限环 ,对应于 Aux/IAA 表达水平的振荡 ;否则，系 统有一个稳定的结点，或者螺旋稳定的 . 1.4 本文的结构 本文考虑从数值方法的角度研究基因调控网络 .具体安排如下 : 5

26、基因调控网络的数值方法 第二章讨论具有全局稳定平衡点的基因调控网络的高效数值方法 .我们将基因调 控微分方程的向量场在平衡点附近分解为线性主部和髙阶非线性部分 ,应用一类特殊 的分裂法 :线性部分精确积分 ,非线性部分用经典的 Runge*Kutta (RK)方法积分 .新方 法用于单基因自调控网络、两基因交叉调控网络以及抑癌蛋白 P53-Mdm2 调控网络的 计算 ,数值结果证明，新方法的计算精度与计算效率明显高于 Euler, Heun, HK4, RK3/8 方法 .经分析，分裂法比相应的 RK 方法有更大的稳定性区域 . 第三章考虑两导数 RungeKutta (TDRK)方法 .这类

27、方法能够精确积分真实解的 一阶与二阶结构 .通过对两类模拟蛋白质磷酸化过程的调控模型 (果蝇周期蛋白模型 与哺乳动物生物钟模型）的数值解法，证明二级四阶 TDRK 方法计算精度和计算效率 明显高于经典的四级四阶 RK 方法 . 第四章考虑时滞基因调控网络 .通过对果蝇周期蛋白调控模型和哺乳动物生物钟 模型的数值模拟探讨时滞对基因调控网络动力学行为可能的影响 . 6 第二章基因调控网络的分裂法 第二章基因调控网络的分裂法 本章研究具有全局稳定平衡点的基因调控网络的数值方法 .我们将微分方程的向 量场在平衡点附近分解为线性主部和非线性高阶项，应用一类特殊的分裂法：线性部 分精确积分，非线性部分用经

28、典的 RungeKutta (RK)方法积分 .经分析 ,分裂法比相应 的RK 方法有更大的稳定性区域 . 2.1 基于 Rung-Kutta (RK)方法的分裂法 或简单的表示为 （ c,A6).如果内级系数矩阵 4 为严格下三角形，即， = 0, 1 S i Trotter 分裂法计 算误差比较 步长 Euler Heun Split(Exact: Euler)Split(Exact: Heun) 0.5 9.3227 x 10-1 3.4257 x 10-2 3.8025 x 1(T4 2.4890 x 102 1.0 1.4748 x 10 3.2048 x 10-1 4.4358 x

29、 1 4 2.4150 x 1 -2 1.5 2.8044 x 10 7.5948 x 10 1.7117x 103 2.5143 x 102 5.0 1.6373 x 1013 2.8356 x 1010 1.1526 x 1 -3 2.3476 x 10_2 表 2.6 两基因交叉调控网络：取不同步长， RK4 法、 RK3/8 法与相应的 Lie-Trotter 分裂法计 算误差比较 步长 RK4 RK3/8 Split (Exact: RK4)Split (Exact: RK3/8) 0.1 9.8188 x 1-2 9.3019 x 102 9.9338 x 1 -4 9.9338

30、x 10_4 1.2 2.5157 x 10-1 5.5798 x 1-2 7.6803 x 103 7.6803 x 103 L5 4.3953 x 10 2.9958 x 10 9.5832 x 10-3 9.5832 x 1 -3 2.0 4.1821 x 1024 7.5238 x 1017 1.6686 x 1 -2 1.6686 x 1 -2 5.0 1.9692 x 1069 3.2225 x 1068 3.1227 x 10 一 2 3.1227 x 1 -2 14 第二章基因调控网络的分裂法 表 2.7 两基因交叉调控网络：取步长 /i = 1.5, Euler 法、 Heu

31、n 法与相应的 Lie-TVotter 分裂法 在不同区间上积分的误差比较 积分区间 Euler Heun SpIit(Exact: Euler)Split(Exact: Heun) 0,100 2.8044 x 10 7.5948 x 10 1.7117 x 103 2.5143 x 10-2 0,500 5.9438 x 10 3.2406 x 10 5.6993 x 107 1.2988 x l -6 0,1000 1.5087 x 10 2.7533 x 10 5.6993 x 107 5.6993 x l -7 表 2.8 两基因交叉调控网络：取步长 ft = 2, RK4 法、 R

32、K3/8 法与相应的 LieTtotter 分裂法 在不同区间上积分的误差比较 积分区间 RK4 RK3/8 Split (Exact: RK4)Split (Exact: RK3/8) 0,500 4.4414 x 10 1.5342 x 10 1.9352 x 10_3 1.9352 x 10_3 0,1000 4.4396 x 10 1.0172 x 10 9.6936 x l -4 9.6936 x 104 0,1500 NaN 2.3844 x 1095 1.1409 x 1CT3 1.1409 x 103 0,2000 NaN NaN 8.5613 x l -4 8.5613 x

33、l -4 0,2500 NaN NaN 6.8520 x l -4 6.8520 x 10_4 15 基因调控网络的数值方法 图 2.4 两基因交叉调控网络的不同方法计算效率比较 (1) 图 2.5 单基因自调控网络的误差比较 (2) 16 第二章基因调控网络的分裂法 图 2.6 单基因自调控网络的计算效率比较 (2) 我们依次用 Euler 法、 Heun 法、 RK4 法、 RK3/8 法以及相应的基于上述向量场分 解的 LieTrotter 分裂法对方程组 （ 2.13)积分 mRNA 和蛋白质初值取为 mO) = 0.6, m2(0) = 0.8, Pl(0) = 0.6, p2(0)

34、 = 0.8 首先取不同步长 ,在区间 0,100上积分，数值结果 如表 2.5和表 2.6 所示 . 然后，对同样的初值，取固定步长 /I = 1.5,在不同区间上用 Euler 法、 Herni 法和相 应的 Lie*Ttotter 分裂法积分，数值结果如表 2.7 所示 .再取固定步长 /I = 2,在不同区间 上用RK4 法、 RK3/8 法和相应的 LieTrotter 分裂法积分 ,数值结果如表 2.8 所示 . 我们在区间 0,500上积分，比较 Euler 法、 Heun 法、 RK4 法、 RK3/8 法以及相应 的基于上述向量场分解的 LieIVotter 分裂法的计算效率

35、 ,数值结果如图 2.4 -图 2.6 所 不 .与传统 Runge-Kutta 方法相比，新的分裂法更适合于大步长长时间积分 . 2.5 分裂法用于 p53-Mdm2 调控网络 P53 是一种肿瘤抑制蛋白，当细胞受到损伤时， P53 就处于激活状态，它能对细胞产 生可变的影响，即能使细胞停歇或修复这些损伤，或细胞凋亡 （ apoptosis),或者细胞能 正常存在但不能分裂 .P53 可以被 Mdm2 ( 种泛素蛋白连接酶 )负调控 ,Mdm2 通过泛 17 基因调控网络的数值方法 图 2.7 p53-Mdm2 调控网络 素化 （ ubiquitination)能让 p53 降解或消失 .图

36、 2.7 是 p53-Mdni2 调控网络图 .图中 M 是Mdm2 蛋白風巧是非激活状态下的 p53 蛋白風 C 是 P53Mdm2 复合物， PA是激 活状态下的 p53 蛋白质 ， ri (i = 1,3,5,6,7,8,9)表示对应的反应速率 ,rm, no 表示降 解率， d表示 P53 的降解过程， IR 表示电离辐射 .图 2.7 中模型的微分方程为 ： 其中 r 糾 ,i = 1,2, ， 10是在时刻 *时第 i 个反应的速率 . 对图 2.7 中的模型作了如下合理的简化假设： (1) n是 P53 的基础合成率 ,假设为常数； 降解率 r2， r4,ri0正比于对应底物的浓

37、度； Mdm2 与 p53 的连接可以假设成是可逆的； 电离辐射可以诱导 P53 的磷酸化，阻止它与 Mdm2 的结合 ; (5) 是一个关于辐射程度的增函数； (6) r3是一个关于 P53 总水平的饱和函数 . (2.19) 18 第二章基因调控网络的分裂法 (2.20) (2.21) 其中 Pi 表示肿瘤抑制器， M (Mdm2)是 p53 主要的负调控子， C1是 P53Mdm2 复合 物，PA是 p53 的激活形式， S 短暂压力刺激物，形式为邓 ） =-ecs = 7fcu， s*(* = p,m0, ml)是 de mwo 的合成率，心卜 =a, c,u)是生成率 ,j*(* =

38、 a, c)是消退反 应 (比如降解 )， rfp是激活状态下的降解率，是饱和系数 此模型参见文献 52. 为了简化模型，我们假设 S=0 方程（ 2.21)的平衡点为 （ P/,M*， C,i)= (9.42094,0.0372868,3.49529,0).因为特征值为的 = 38.4766, w2 = 0.0 28 + 0. 22 i, = -0.0028 - 0.0220i,叫 =-0.2002 的实部全为负数，因此，平衡点是渐近稳定的 .对 此系统作如下分解 (2.22) 19 基因调控网络的数值方法 我们取初值乃 ( ) = 9.42 nM, C7(0) = 0.037 nM, M(

39、0) -： 3.49 nM, PA(0) = 0 nM, 用 RK4法、 RK3/8法以及相应的 LieTrotter分裂法在区间 0,100上积分 .数值结果如 表2.10 和表 2.11 所示 .在积分区间 0,500上积分的数值结果如图 2.S 和图 2.9 所示 . 表 2.10 p53-Mdm2 调控网络：以不同步长用 RK4 法、 RK3/8 法以及相应的 LieTrotter 分 裂法计算误差比较 stepsize RK4 RK3/8 Split (Exact: RK4) Split(Exact: RK3/8) 0.05 3.7686 x l -5 3.7046 x 105 1.

40、2091 x 106 1.2065 x 10 一 6 0.08 4.6118 x 10 4.5057 x 10 2.1383 x 106 2.1290 x 10-6 0.10 NaN 5.3550 x 10 2.8098 x 10_6 2.7945 x l -6 0.12 NaN NaN 3.4986 x l -6 3.4762 x l -6 5.00 NaN NaN 7.3001 x l -4 3.8208 x 10 一 4 从表 2.10 和表 2.11,以及图 2.8和图 2.9， 我们看出对较小的步长 Split (Exact: RK4) 和 Split (Exact: RK3/8)比

41、 RK4 和 RK3/8 更为精确，更为高效；当步长增大到 RK4 和 RK3/8 失效时， Split (Exact: RK4)和 Split (Exact: RK3/8)仍然能高精度地计算 20 第二章基因调控网络的分裂法 表 2.11 P53-Mdm2 调控网络：以步长 A = 10 用 RK4 法、 RK3/8 法以及相应的 Libltotter 分裂法计算误差比较 time interval RK4 RK3/8 Split(Exact: RK4)Split (Exact: RK3/8) 0,100 NaN NaN 6.8810 x 1(T2 6.6940 x 102 0,500 Na

42、N NaN 2.1634 x l -2 2.0911 x 102 0,1000 NaN NaN 1.1625 x l -2 1.1214 x l -2 0,1500 NaN NaN 7.8314 x 103 7.5529 x l -3 0,2000 NaN NaN 5.8881 x 103 5.6786 x 103 图 2.8 p53-Mdm2 调控网络的不同方法计算效率比较 21 基因调控网络的数值方法 图 2.9 p53-Mdm2 调控网络的计算误差比较 2.6 RK 方法和分裂法的稳定性分析 本小节分析 RK4、 RK3/8 方法和相应的 Liyitotter 分裂法的稳定性，用以解释分

43、 裂法相对于传统 Runge-Kutta 方法的优势 . 首先 ,考虑下面的测试方程 y = Xy + ey, (2.23) 其中 A 是测试率， e 是估计误差 .当数值方法应用到测试问题 （ 2.23)中时，我们可以得 到 J/n+l (2.24) 称及 (14， 1； )为稳定性函数 .w = eh. 定义 2.1. u-u 平面区域 咒 S = (乜，以 )丨丨雕， U)| S 1 (2.25) 是方法 的稳定性区域，由 |i?(u,u)| $ 1 所定义的曲线 称为稳定性域的边界 . 可以得到 RK 方法的稳定性函数如下 i?(w, v) = 1 + (w + v)bT(I (u +

44、 v)A)le 22 第二章基因调控网络的分裂法 图 2.11 (a) Spit (Exact:RK4)方法的稳定性区域 (左） （ b) Spit (Exact:RK3/8)方法的稳定性区 域 (右） 其中 e = (1， 1， .， If , J 是 s x s单位矩阵，分裂法的稳定性函数为 R(u, v) = exp(u)(l + ubT(I Ji4)1e). 我们得到 RK4 和 RK3/8 的稳定性区域如图 2.10, Split (Exact:RK4)和 Split (Exact: RK3/8) 稳定域如图 2.11. 容易看出 RK4 和 RK3/8 要比 Split(Exact

45、: RK4) 和 Split(Exact: RK3/8)的稳定域小得多 .这说明从稳定性的目的看，前者比后者有更多 的步长限制，即当步长超过某个值时， RK4和 RK3/8将失效，而 Split(Exact: RK4)和 Split(Exact: RK3/8)仍然能够有效地计算 . 23 基因调控网络的数值方法 24 第三章含有磷酸化过程的基因调控网络的两导数 Riing+Kutta 方法 第三章含有磷酸化过程的基因调控网络的两导数 Runge-Kutta 方法 本章主要研究两类模拟蛋白质磷酸化过程的调控模型（果蝇周期蛋白模型与哺乳 动物生物钟模型 )的数值解法 .我们考虑两导数 Runge_Kutta (TDRK)方法 .这类方法精 确积分真实解的一阶与二阶结构，计算因而计算精度明显髙于传统的 Runge-Kutta (RK) 方法 .通过对两种含有憐酸化过程的基因调控网络的数值计算 ,证明二级四阶 TDRK 方 法（