数学建模在昼夜节律生物钟中的应用*

李莹 郑明银 刘曾荣

（1.上海海洋大学信息学院，上海 201306）（2.泛亚汽车技术中心有限公司，上海 201201）（3.上海大学系统生物技术研究所，上海 200444）

引言

生物理论模型的发展具有很长的历史，早在遗传学和分子生物学出现之前，人们就利用概念上的数学模型来研究生物群体和生物系统中的一些现象特征.随着生命过程中复杂动力学分子调控机制的发现，人们基于实际的调控过程建立了详细的模型.从神经元活性到代谢、细胞周期、昼夜节律生物钟等细胞过程，人们建立了各种各样的动力学模型，这些模型仍然是理论研究的焦点.

生物钟的基因调控网络是由一些钟基因和蛋白以及它们之间的调控作用、转录、后转录（比如多磷酸化、蛋白复合物的形成等）组成.目前，已经发现了几种由正负反馈环组成的转录机制［1］，这些环路使得单个细胞能产生以24小时为周期的节律振子，细胞间的耦合机制使得细胞振子在组织水平上达到同步［2，3］.昼夜节律生物钟涉及到有机体的各个不同水平，遗传学家、分子生物学家、结构生物学家、生理学家和数学家等都对生物钟进行了广泛的研究.早在第一个时钟基因被发现之前，数理学者们就在探索生物钟建模的方法.在1960年前，有关昼夜节律生物钟的文献就有600多篇［4］.现在，生物钟的理论模型已经成为理解复杂昼夜节律生物钟的一个普遍而有效的方法.

数学模型可以帮助我们在单个细胞和细胞群体之间掌握和理解昼夜节律生物钟的动力学特征.通过数值模拟，数学模型可以分析关键参数在系统中的作用，可以预测在没有被实验验证的条件下的系统行为.数学模型也可以为实验中的直观发现提供合理的解释，或者揭示生物钟分子机制的设计原理.

尽管可以利用数学模型和计算机等工具来研究生命系统，但是没有现成的方法构建数学模型.实际上，当我们要写方程构建模型时，有很多问题需要考虑.建模时所做的选择不仅取决于所研究的生物问题，也与建模人员的专业背景和经验有关.本文将介绍用来研究昼夜节律生物钟动力学性质的各类模型以及分析模型的方法，然后介绍数学模型在研究生物钟中的优势、局限性和未来面临的挑战.

1 模型建立

1.1 模型框架

每年都有几百篇论文报到新的钟基因的发现、新的基因突变以及新的调控关系等，但是在昼夜节律生物钟研究中并不能包含所有的信息.建立数学模型需要收集和整合所研究系统涉及到的各种信息，区分出系统的关键成分，并从中找出系统各个成分之间的关键作用和连接.在昼夜节律生物钟领域，为了建立模型甚至在写方程之前，要画一个包含这些信息的框架图.

构建框架要求我们认清系统中各组成成分之间相互作用的本质并整理不同的分子过程.比如：假定一个蛋白可以被磷酸化并在细胞核中担任转录因子.这个信息或许来自不同的文献，但是我们并不清楚到底是这个蛋白先进入细胞核再磷酸化还是磷酸化之后再进入细胞核.建立模型就要将各种不同的分子过程进行连接和整理.另外一个问题是如何处置蛋白复合物.为了降低模型的复杂度，我们需要识别哪些蛋白复合物在生物钟中起关键作用［5］.对这类问题实验证据的缺失使得研究者们或者考虑简化后的模型（比如不区分细胞质和细胞核中蛋白的形式，并且只考虑有限几个蛋白复合物），或者需要验证或辨别各种模棱两可的生物学机制.

1.2 从框架到方程

一旦有了模型框架，我们要将其转化成数学方程.方程通常用来描述关键化合物（变量）的浓度或活性随时间的演化过程.模型可以是概念上的（也就是一般意义上的）也可以是根据生物有机体的精确生命过程建立起来的［6］.通常，模型是由描述变量动力学特性的常微分方程组成的方程组，模型类型见表 1［7］.

相振子方程是最抽象最直观的一类模型，其中变量是振子的振幅和相位，这类方程可以独立于分子机制来描述振子的动力学行为.相振子方程可以用来研究一般的动力学特性，比如相响应曲线［8］、日常光照的导引［9］、耦合振子的同步［10］、或为实验数据提供合理的解释［11］.另外一类是现象模型，其中的变量和参数与分子化合物没有直接联系，这类方程可以用来研究振子的一般特性.Van der Pol模型是一个物理模型［12］，研究者将其适当修改以产生节律周期以及响应光照的导引，用来描述昼夜节律生物钟［13］，进而研究哺乳动物生物钟振子的重置特性［14］、蓝藻生物钟振子与细胞分裂振子之间的相互作用［15］以及生物钟振子的同步［16］.

早在昼夜节律生物钟的分子机制被揭示之前，研究者们就为生物钟提出了一般的分子模型，其中最著名的就是由Goodwin于1965年提出的Goodwin振子模型［17］，该模型是一个含有一个负反馈环的最小的振子模型.二十世纪九十年代，果蝇和脉孢菌的生物钟转录调控的负反馈环被第一次证实［18，19］.随后，Ruoff及其合作者在 Goodwin振子的基础上建立了其微分方程模型［20］.这个简单的模型或其变形，至今仍被用来分析核心节律振子的根本性质或者耦合节律振子的同步问题［21-23］.

表1 昼夜节律生物钟数学模型的分类Table 1 Types ofmathematicalmodels proposed for circadian clocks

随着分子生物学的发展，生物钟的分子机制越来越清楚，反映生物细节的模型也越来越多.随着新的时钟基因和蛋白的发现，基于分子机制的模型也在不断被更新.模型中用5个到70个不等的变量描述参与生物钟的mRNA和蛋白的浓度，浓度的演变方程是基于合成（转录或翻译）、降解、修正（磷酸化等）、复合物的合成或分解以及细胞成分的传输等建立的.这些过程的动力学是基于标准的生化反应以及酶动力学来建立的，一般用Michaelis-Menten方程或者Hill函数来描述，模型中的参数分别代表合成率或降解率、传输率、结合亲和力等.目前，对哺乳动物、脉孢菌、果蝇、拟南芥和蓝藻等生物的生物钟已经有大量的计算模型［24-30］.详细的数学模型是依据详细的生物学机制建立的，模型变量个数较多，在此就不列出，可查阅相关文献.

建模方法的选择取决于所研究的问题.一般的模型用来研究昼夜节律生物钟的一般特性，比如大量振子之间的耦合，或者研究受生物钟调控的现象.在这些问题中，能产生昼夜节律振子的模型很重要，但是相关的分子机制却不是必需的.反之，如果研究的是生物钟的结构与相响应曲线的形状之间的关系，或者是某个基因突变对生物钟的影响，我们则需要一个详细的模型来明确地描述光照的作用或者与突变相关的参数.

2 模型的分析

2.1 数值模拟及参数的选取

大部分的模型中都包含若干个变量，其时间演化方程中含有Michaelis-Menten方程以及Hill函数等非线性项.当变量个数较多时，这类方程一般不能通过分析的方法求解.这就需要用数值模拟的方法求系统的数值解，进而研究系统的动力学性质.然而，数值模拟要求方程中的每个参数的值取定，这个要求往往是比较难办.因此，大家通常关注系统在不同参数取值条件下的动力学行为，分岔图是研究系统动力学与控制参数之间关系的一个有力工具.

目前有许多软件可以用来对非线性微分方程组进行数值模拟求数值解.其中被建模者用的较多的软件是XPPAUT［31］，该软件可以提供各种综合方法，并且可以用该软件画分岔图.其它软件还有Grind［32］、Berkeley Madonna和 COPASI，以及程序语言比如Matlab或 Mathematica.

系统的动力学行为依赖于参数取值，但是我们并不能通过改变参数取值得到系统任意的动力学行为，只能根据系统的一些特征把参数的取值缩小到一定范围内.为了得到极限环振子，系统应该包含一个负反馈环以及非线性项.在选取参数时除了要考虑周期近24小时外，其它的实验发现也需要考虑，比如振子的振幅、变量之间的相关系、突变的影响、被光照的导引以及相响应曲线等.满足条件的参数值往往是通过手工微调得到的，有时也用自动化方法来选取参数值，这种方法一般依赖于最优化方法.找到一组合适的参数值也不是永远有效，只是能使系统产生具有鲁棒性的振子，当参数在一定范围内变化时，其周期和振幅不变或者只有轻微的变化.为了量化参数变化的影响，研究者们提出几种敏感性分析并且应用在节律模型中［33］.

2.2 极限环振子、分岔图及相响应曲线

昼夜节律振子是一种内在的周期活动，即使在没有外界刺激下，节律仍然存在.昼夜节律生物钟表现出的振子行为是系统自身产生的自维持振子，而不是系统被动地响应外界周期环境信号产生的.

大部分的昼夜节律生物钟模型都产生极限环振子，其特点是有特定的周期和振幅，且与初始条件无关.极限环振子使得系统在环境小扰动下具有鲁棒性.尽管在实验中发现的振子类型可能依赖于细胞类型、组织以及环境条件等，很难被刻画，但是理论学家们一般都用极限环模型来刻画它们［34，35］.

通常极限环的特征随着参数值变化而变化，这种变化可以通过分岔图展示出来.分岔图描述的是振子的振幅或者是变量的最大值、最小值或者是振子的周期与控制参数的函数关系.通过数值模拟可以画出每个变量关于控制参数的分岔图，如果系统收敛到稳态，则变量的最大最小值将在分岔图上重合成一个点.振子形状的变化以及周期随着控制参数的变化不仅依赖于系统的结构，也依赖于其他参数值，单纯的直观预测是很难得到的.

2.3 光照导引及相响应曲线

在昼夜节律研究中，光照始终是一个至关重要的因素.在自然条件下，生物有机体受到日夜交替的周期光照影响.为了适应这种周期环境，昼夜节律生物钟被光照周期导引（当振子周期等于光照周期时，称系统被导引）.

光照对模型的影响可以通过调整与光照相关的参数值进行模拟.比如，为了模拟日夜交替光照对哺乳动物生物钟振子的影响，per基因的转录速率在白天取较大的值，在晚上取较小的值.换句话说，通过对模型中依赖于光照的转录速率加一个周期驱动来模拟日夜交替的周期光照.当节律振子被周期驱动导引时，相对于外部驱动而言，振子的相位是固定的，也称振子被相锁.通常，当日夜交替周期光照的周期与节律振子的周期非常接近时，周期振子可以被导引并与周期光照信号相锁［36］.当周期光照信号的振幅比较小或者其周期与节律振子周期差别较大时，振子不能被导引而表现出非规则振动及非相锁等复杂的动力学行为.

光照脉冲通常诱导振子的相变，脉冲加入系统可能导致振子的相位超前或滞后.相变的符号和振幅关于加入脉冲时刻的函数图像称为相响应曲线（PRC），该曲线是光照对生物钟的影响独有的特征.光照脉冲可以通过在一个短时间周期内提高相关参数的值来实现.相变可以通过比较振子在光脉冲下的相位以及没有光脉冲下的相位得到.画出相变和脉冲加入时刻之间的函数关系就得到相响应曲线，并且可以将画出的相响应曲线与实验得到的相响应曲线对比［37］.

3 模型评价

3.1 模型的优势

理论模型可以提供新的方法用来研究复杂系统的动力学特征，并且有助于解释很难被实验解释的一些热点问题.数学模型已经被用于研究昼夜节律生物钟的基本特征以及生物钟与生理和行为之间的关系［38］.

模型的另一个优势是可以突出生物系统基本的设计原理.比如节律振子必须对温度、参数变化或者分子噪声鲁棒，这与其它遗传振子是类似的，于是这些不同的系统可以共享共同的调控结构［39，40］.研究发现，许多振子都是由相锁的正负反馈环组成.通过对这些系统的系统性研究，我们可以根据它们的正负反馈环回路的拓扑结构对它们进行分类［39］，并且将它们基本的动力学性质（比如：正弦振子或松弛振子、分岔结构、周期的可调性、鲁棒性等等）与拓扑类型相对应［39，40］.

3.2 模型的局限性

模型最普遍的局限性是方程的类型以及参数值的选取都是任意的.分子模型是基于确定的基因调控过程建立的，而分子机制的精确过程通常是不知道的.经常用Michaelis-Menten和Hill函数是因为饱和度和动力学阈值是酶动力学过程的现实描述，并且它们提供了产生极限环和其它非平凡行为所必需的非线性，但是这些描述背后的假设却很少被证明是满足的.分子水平上的理论模型可以研究昼夜节律生物钟量化的动力学性质，而量化数据（如参数值、绝对浓度或分子数量）很少可以被利用.因此，不可能为昼夜节律生物钟建立一个精确的模型.

昼夜节律生物钟的大部分模型是常微分方程模型，这些模型都忽略了细胞的空间结构.这些模型隐含的假设是mRNA和蛋白分子可以在细胞内自由移动.然而，细胞内是非常拥挤的，细胞过程在空间中被高度组合.空间和扩散似乎在细胞系统的动力学中起了关键作用.忽略扩散或许是一个粗略的近似，这种近似使得量化模型变得具有很大争议.然而，预测行为和实测数据之间的良好匹配可以提高模型及其动力学特征的可信赖度.

最后，我们不得不认识到，再具体的模型都不可能对问题给出非常明确的回答.然而，它可以提供新的思路供我们思考.尽管正反馈环在昼夜节律振子中的作用没有明确解释，但是不同的模型都可以对这些环路的作用提供可能的提示，比如提高系统对参数变化的鲁棒性、增加周期的可调性等等.所以，模型使得我们可以从一个新的角度看待和思考问题.

4 结论

本文提到的不同模型说明生物钟的复杂结构和它的动力学特征之间存在着难以理解的关系.目前，仍有很多公开的问题有待解决.针对昼夜节律生物钟分子机制的建模越来越多，随着模型的发展，新的挑战也逐渐浮出.分子模型不断地被加入最新的实验数据，变得越来越具体.利用这些模型可以研究节律中非常明确的问题，比如具体的突变或者某个特定磷酸化的作用.最近的细胞模型不仅包含了生物钟的基因调控网络，而且还包含了相关的成分，比如神经递质、信号通路、钙和其它离子的传输过程［41］.未来的模型将把详细的分子机制与实际的细胞网络结构相结合，这种模型毫无疑问地将为昼夜节律生物钟的功能提供新的观点.

近年来，越来越多的昼夜节律研究的实验团队与建模协力进行实验研究.基于新的实验数据模型得以进一步改进，并且基于模型的实验也将带来新的数据，这样的相互促进迭代过程在昼夜节律研究中越来越普遍.建模也成了合成生物学中的关键步骤，研究者们已经从实验上设计和构造出了人造生物钟［42］.设计一个鲁棒可调的生物钟显然依赖于在生物钟建模中所获取的所有知识.

1 Shearman L P，Sriram S，Weaver D R，et al.Interacting molecular loops in themammalian circadian clock.Science，2002，288：1013～1019

2 Hastings M H，Herzog E D.Clock genes，oscillators，and cellular networks in the suprachiasmatic nuclei.Journal of Biological Rhythms，2004，19：400～413

3 Reppert SM，Weaver D R.Coordination of circadian timing inmammals.Nature，2002，418：935～941

4 Roenneberg T，Chua E J，Bernardo R，etal.Modeling biological rhythms.Current Biology，2008，18：826～835

5 Yamada Y，Forger D.Multiscale complexity in the mammalian circadian clock.Current Opinion in Genetics＆Development，2010，20：626～633

6 Roenneberg T，Chua E J，Bernardo R，etal.Modelling biological rhythms.Current Biology，2008，18：826～835

7 Gonze D.Modeling circadian clocks：From equations to oscillations.Central European Journal of Biology，2011，6（5）：699～711

8 Granada A E，Herzel H.How to achieve fast entrainment.The timescale to synchronization.PLoS One，2009，4：e7057

9 Roenneberg T，Dragovic Z，Merrow M.Demasking biological oscillators：properties and principles of entrainment exemplified by the Neurospora circadian clock.Proceedingsof the National Academy of Sciences of the United States of A-merica，2005，102：7742～7747

10 Rougemont J，Naef F.Collective synchronization in populations ofglobally coupled phase oscillatorswith drifting frequencies.Physical Review E，2006，73：011104

11 Yang Q，Pando B F，Dong G，et al.Circadian gating of the cell cyclerevealed in single cyanobacterial cells.Science，2010，327：1522～1526

12 Van der Pol B.A theory of the amplitude of free and forced triode vibrations.Radio Review，1，1920，701～710，754～762

13 Pittendrigh，C S，Bruce V G.An oscillatormodel for biological clocks.In：Rudnick D editor.Rhythmic and synthetic processes in growth.Princeton NJ：Princeton University Press，1957：75～109

14 Forger D B，JewettM E，Kronauer R E.A simplermodel of the human circadian pacemaker.Journal of Biological Rhythms，1999，14，532～537

15 Gonze D，RousselM R，Goldbeter A.Amodel for the enhancement of fitness in cyanobacteria based on resonance of a circadian oscillator with the external light-dark cycle.Journal of Theoretical Biology，2002，214：577～597

16 Kunz H，Achermann P.Simulation of circadian rhythm generation in the suprachiasmatic nucleus with locally coupled self-sustained oscillators.Journal of Theoretical Biology，2003，224：63～78

17 Goodwin B C.Oscillatory behavior in enzymatic control processes.Advances in Enzyme Regulation，1965，3：425～438

18 Hardin P E，Hall JC，Rosbash M.Feedback of the Drosophila period gene producton circadian cycling of itsmessenger RNA levels.Nature，1990，343：536～540

19 Aronson B D，Johnson K A，Loros J J，et al.Negative feedback defining a circadian clock：autoregulation of the clock gene frequency.Science，1994，263：1578～1584

20 Ruoff P，Mohsenzadeh S，Rensing L.Circadian rhythms and protein turnover：the effectof temperature on the period lengths of clock mutants simulated by the Goodwin oscillator.Naturwissenschaften，1996，83：514～517

21 Li Y，Liu Z R，Zhang JB.Phase synchronization of circadian oscillators induced by a light-dark cycle and extracellular noise.Lecture Notes in Computer Science，2007，4689：1～7

22 Li Y，Liu ZR，Zhang JB，etal.Synchronizationmechanisms of circadian rhythms in the suprachiasmatic nucleus.IET Systems Biology，2009，3（2）：100～112

23 Li Y，Liu ZR.Synchronization of clocks coupled by neurotransmitter in the SCN.Journal of Biological Systems，2013，21（1）：1350006（12）

24 Leloup JC，Goldbeter A.Towards a detailed computationalmodel for the mammalian circadian clock.Proceedings of the National Academy of Sciencesof the United States of America，2003，100：7051～7056

25 Francois P.A model for the Neurospora circadian clock.Biophysical Journal，2005，88：2369～2383

26 Hong C I，Jolma IW，Loros J J，et al.Simulating dark expressions and interactions of frq and wc-1 in the Neurospora circadian clock.Biophysical Journal，2008，94：1221～1232

27 Leloup JC，Goldbeter A.A model for circadian rhythms in Drosophila incorporating the formation of a complex between the PER and TIM proteins.Journal of Biological Rhythms，1998，13：70～87

28 Locke JC，Millar A J，Turner M S.Modelling genetic networks with noisy and varied experimental data：the circadian clock in Arabidopsis thaliana.Journal of Theoretical Biology，2005，234：383～393

29 Gen K，Kazuyuki A，Yoh I.A model for the circadian rhythm of cyanobacteria that maintains oscillation without gene expression.Biophysical Journal，2006，91：2015～2023

30 Li Y，Wu H，Luo JH.Comparing two models based on the transcriptional regulation by KaiC of cyanobacteria rhythm.Xi’an：In IEEE 6th International Conference on Systems Biology，2012：216～219

31 Ermentrout B.Simulating，analyzing and animating dynamical systems：a guide to XPPAUT for researchers and students.Philadelphia：SIAM，2002：1～286

32 Hoops S，Sahle S，Gauges R，et al.COPASI-a complex pathway simulator.Bioinformatics，2006，22：3067～3074

33 Leloup JC，Goldbeter A.Modeling themammalian circadian clock：sensitivity analysis and multiplicity of oscillatory mechanisms.Journal of Theoretical Biology，2004，230：541～562

34 Leloup JC，Gonze D，Goldbeter A.Limit cycle models for circadian rhythms based on transcriptional regulation in Drosophila and Neurospora.Journal of Biological Rhythms，1999，14：433～448

35 Li Y，Liu Z R，Zhang J B.Circadian oscillators and phase synchronization under a light-dark cycle.International Journal of Nonlinear Science，2006，1（3）：131～138

36 Gonze D，Goldbeter A.Entrainment versus chaos in a model for a circadian oscillator driven by light-dark cycles.Journal of Statistical Physics，2000，101：649～663

37 Leloup JC，Goldbeter A.Modeling themolecular regulatorymechanism of circadian rhythms in Drosophila.BioEssays，2000，22：84～93

38 李莹，刘曾荣.哺乳动物生物钟的数学建模及研究进展.动力学与控制学报，2012，10（3）：263～267（Li Y，Liu ZR.Advances andmathematicmodeling ofmammalian circadian clock.Journal of Dynamics and Control，2012，10（3）：263～267（in Chinese））

39 Novak B，Tyson J J.Design principles of biochemical oscillators.Nature ReviewsMolecular Cell Biology，2008，9：981～991

40 Tsai T Y，Choi Y S，Ma W，et al.Robust，tunable biological oscillations from interlinked positive and negative feedback loops.Science，2008，321：126～129

41 Vasalou C，Henson M A.A multiscale model to investigate circadian rhythmicity of pacemaker neurons in the suprachiasmatic nucleus. PLOS Computational Biology，2010，6：e1000706

42 Tigges M，Marquez-Lago T T，Stelling J，etal.A tunable syntheticmammalian oscillator.Nature，2009，457：309～312