系统生物学数学基础 84页

系统生物学数学基础(初稿)雷锦誌清华大学周培源应用数学研究中心2007年9月\n系统生物学数学基础前言什么是系统生物学?“Systemsbiologyisthescienceofdiscovering,modeling,understandingandultimatelyengineer-ingatthemolecularlevelthedynamicrelationshipsbetweenthebiologicalmoleculesthatdefinelivingorganisms.”LeroyHood,Ph.D.,M.D.,PresidentInstituteforSystemsBiologyAllbiologicalphenomena,whetherit’sdigestionofasugarmolecule,beatingofthehumanheart,orneutralizinganinvadingvirus,aretheresultofcomplexsystems.Thusourapproachistofocusresearchonbiologicalsystemsasawhole,ratherthanpursuethetraditionalapproachoffocusingonindividualgenes,proteins,orpartsofanorganism.系统生物学的研究内容?系统生物学研究方法?Scientistsfrommultipledisciplines(biology,chemistry,mathematics,physics,etc.)workcloselytogethertofullyunderstandallaspectsoftheinherentlycomplexsystemsintrinsictolivingorganisms.Suchin-depthunderstandingisultimatelyessentialtorealizingourgoalofpredictive,preventive,personalizedmedicine.系统生物学与数学?•“它(物理学)的范畴可定义为我们全部知识中能够用数学语言表发的那个部分”－－爱因斯坦•统计，模型，分析，模拟......总结定性、半定建立数学模收集经验数据--性规律型@I@@@?用经验资料验证求解、发展数学模型理论“Mostreadersofthispublicationwillknowthat‘post-genomics’and‘proteomics’arephrasesthatmeanlittlethatisspecificbutheraldanencyclopaediceraofinformationaboutthewaybiologicalcellsandtheirgenesandproteinsbehave.Buthowbesttomakesenseofitall?Itis,atlast,possibletoanticipatemathematicsbecomingusefulinthemodellingofthesystems.”–Nature4072000,819.内容简介?•基因表达•基因调控网络–Toggleswitches–生物振荡2\n系统生物学数学基础–生命节律–胚胎发育•细胞分裂与分化•系统生物学前沿介绍•随机过程(Masterequation,Langevinequation,Fokker-Plankequation)•微分方程(建模，定性理论，数值求解)•随机微分方程(建模，数值求解，稳定性分析)•反应扩散方程(建模，数值求解)补充阅读材料:1.Mackey,M.C.,Santillan,M.,Mathematics,Biology,andPhysicss:Interactionsandinterepen-dence,NoticesAMS,52(2005)(8).2.Sontga,E.D.,Molecularsystemsbiologyanddynamics:anintroductionfornon-biologists.3.Alon,U.,Anintroductiontosystemsbiology,Chapman&Hall/CRC,London,2007.4.Fall,C.P.,Marland,E.S.,Wagner,J.M.,Tyson,J.J.,(Eds.)Computationalcellbiology,Springer-Verlag,NewYork,2001.5.Alberghina,L.,Westerhoff,H.V.(Eds.)Systemsbiology:Definitionsandperspectives.Springer,Berlin,2005.3\n目录第一章化学反应的数学描述...................................................1§1.1引言..............................................1§1.2常微分方程..........................................2§1.3化学主方程..........................................2§1.4化学朗之万方程........................................5§1.5计算模拟...........................................7§1.5.1常微分方程的数值模拟................................7§1.5.2求解化学主方程...................................7§1.5.3求解Fokker-Plank方程...............................8§1.5.4求解化学朗之万方程.................................8§1.6Michaelis-MentenandHillEquations............................9§1.7补充阅读材料.........................................12第二章基因表达...................................................................13§2.1引言..............................................13§2.2实验事实...........................................13§2.3数学模型...........................................13§2.4基因表达的随机性......................................16§2.5反馈控制...........................................16§2.6补充阅读材料.........................................16第三章基因调控...................................................................18§3.1ToggleSwitches........................................18§3.1.1Bistability.......................................18§3.1.2Amodelforrepressorexpression[21]........................20§3.1.3Noiseinduceswitches–extrinsicnoise........................22§3.1.4Noiseinduceswitches–intrinsicnoise........................23§3.2生物振荡...........................................26§3.2.1AtkinsonOscillator..................................26§3.2.2Asyntheticgene-metabolicoscillator........................28§3.2.3Mechanismsofnoise-resistanceingeneticoscillators...............29§3.2.4振荡的同步......................................34§3.2.5常微分方程定性分析.................................39§3.3生命节律...........................................42§3.3.1DimerizationandproteolysisofPERandTIM..................42§3.3.2Circadianrhythmgenerator.............................45§3.4胚胎发育...........................................49§3.5Morphogengradient.....................................50§3.6模型..............................................52§3.7二阶常微分方程边值问题的数学基础............................57i\n§3.7.1解的存在唯一性...................................57§3.7.2问题的求解......................................60第四章神经科学...................................................................63§4.1离子通道与Nernst方程...................................63§4.2细胞膜模型..........................................64§4.3离子通道的激活与失活....................................64§4.4Morris-Lecar模型......................................64§4.5Hodgkin-Huxley模型....................................64第五章细胞增生与分化..........................................................65§5.1一些数据...........................................65§5.2造血干细胞的数学模型与参数估计.............................66§5.2.1数学模型.......................................66§5.2.2参数估计.......................................68§5.3造血干细胞模型的动力学分析................................71§5.4整体模型...........................................72第六章细胞调亡...................................................................74第七章蛋白质折叠与随机动力学.............................................75参考文献...................................................................................76名词索引...................................................................................78插图索引...................................................................................79\n第一章化学反应的数学描述{ch3}§1.1引言考虑N≥1种分子{S1,···,SN}的化学反应.假设所有分子充分混合在一个体积(体积为Ω)固定的容器中.并且假设温度不变。共有M≥1个反应通道{R1,···,RM}.以X(t)=(X1(t),···,XN(t))记该系统在时刻t的状态,其中Xi(t)表示系统中分子Si在时刻t的个数,(i=1,···,N).系统状态X(t)是一个随机过程,这是因为每个化学反应发生的时间是随机的,受热力学涨落的映像.对每个反应通道,可以定义相应的PropensityFunctionaj,使得aj(x)dt表示给定X(t)=x时,反应Rj在时间区间[t,t+dt)内,在容器中某处发生一次的概率(j=1,...,M).每个反应都会引起分子个数的改变.定义反应通道Rj的状态改变向量(state-changevector)vj如下:vji表示分子Si因为反应Rj所引起的改变量(j=1,···,M;i=1,···,N).这里vji>0表示反应Rj产生分子Si,vji<0表示反应Rj消耗分子Si.函数aj和向量vji一起给出了反应通道Rj的完整描述.根据这些量,我们可以建立数学模型,用于描述所研究的化学反应.对于两个分子的反应,Propensityfunctionaj(x)可以按以下形式定义:aj(x)=cjhj(x).(1.1.1)其中cj是反应通道Rj的specificprobabilityrateconstant,定义为cjdt表示随机选取的反应Rj的一对反应物分子,在下一个无穷小时间dt内发生反应的概率.这个概率等于这一对分子在下dt时间内发生碰撞的概率,乘以这对已经发生碰撞的分子确实发生了反应Rj的概率.常数cj也可以用反应数率常数(reactionrateconstant)kj来表示.函数hj(x)表示在状态X(t)=x时,可以发生反应Rj的所有不同的反应物的组合数.(相应地,我们有v1=(+1,−1,0,···,0)和v2=−v1.)如果Rj是单分子反应,上面的讨论仍然适用,但是cj与相应分子的量子性质,例如降解率有关.而hj(X)等于相应分子的个数.例:对以化学反应R1:X1+X2→2X1,我们有a1(x)=c1x1x2.而逆过程R2:2X1→X1+X2对应的propensityfunction为a2(x)=c2x1(x1−1)/2.1\n系统生物学数学基础§1.2常微分方程假设在t时刻,系统的状态为X(t).则在时间区间[t,t+dt)内,反应Rj发生的概率为aj(X)dt.如果发生了反应,分子Si的个数的改变量为Xi(t)+vji.因此,在t+dt时刻,分子Si的个数平均改变量为XMXi(t+dt)−Xi(t)=aj(X)vjidt(i=1,···,N).j=1两边除以dt,并且令dt→0,我们得到相应的常微分方程XMdXi=aj(X)vji(i=1,···,N).(1.2.2){eq:1.2:1}dtj=1在这里,分子的个数是整数,因此Xi(t)是整数值的函数.因此上面对时间t的微分数学上并不严格成立.然而,如果分子的个数充分大,上面公式可以给出很好的近似.通常,在描述化学反应时,适用浓度描述系统的状态,即Z=(Z1,···,ZN),其中Zi=Xi/Ω.此时上面方程可以表示为化学速率方程(chemicalrateequation)XMdZi=a˜j(Z)vji(i=1,···,N).(1.2.3){eq:cre}dtj=1其中aj(ΩZ)a˜j(Z)=.Ω例:对前面的例子,我们有常微分方程模型dZ1=k1Z1Z2−k2Z1(Z1−1/Ω)dt(1.2.4){eq:1.2:3}dZ2=−k1Z1Z2+k2Z1(Z1−1/Ω)dt这里k1=c1Ω,k2=c2Ω/2分别为反应R1和R2的反应速率常数.当Ω→+∞时,我们得到了孰知的方程(chemicalrateequation)dZ12=k1Z1Z2−k2Z1dt(1.2.5){eq:1.2:4}dZ22=−k1Z1Z2+k2Z1dt由上面的例子可以看到,方程(1.2.5)仅仅是当体积充分大(当然分子数也充分大)时的近似形式.而对于大量的生物学问题,相应的反应都是在细胞内完成的,体积并不是很大,而且分子数也不大,所以上面方程只能得到近似的描述,而且还可能得到错误的描述.因此,我们需要好的数学模型.§1.3化学主方程以为系统的状态随时间的变化是随机过程,为了得到更加精确的描述,我们使用随机描述来建立数学模型,既化学主方程(chemicalmasterequation).定义条件概率函数P(x,t|x0,t0)如下:P(x,t|x0,t0)=Prob{X(t)=x,giventhatX(t0)=x0}.(1.3.6){eq:1.3:1}2\n系统生物学数学基础这样,化学反应的动态过程可以通过函数P(x,t|x0,t0)随时间的变化描述出来.为此,取dt充分小,使得在时间dt内发生两次或者更多次化学反应的概率可以忽略.这样,我们可以根据t时刻的条件概率写出t+dt时刻的概率:XMXMP(x,t+dt|x0,t0)=P(x,t|x0,t0)×1−aj(x)dt+[P(x−vj),t|x0,t0)aj(x−vj)dt].j=1j=1令dt→0,我们就可以得到化学主方程XM∂P(x,t|x0,t0)=[aj(x−vj)P(x−vj,t|x0,t0)−aj(x)P(x,t|x0,t0)].(1.3.7){eq:cme}∂tj=1方程(1.3.7)从本质上反应了我们所研究的系统.如果可以求解出P,则可以完整地刻划随机过程X(t).然而,除了很特殊的情况,方程(1.3.7)的精确解一般是得不到的,当分子数或者是反应的数量很大时,即使是数值解也不容易得到.为了进一步研究前面的常微分方程模型的含义,我们来看X(t)作为随机过程的平均量的动力学.为此,定义条件期望XE(t|x0,t0)=xP(x,t|x0,t0).x表示当X(t0)=x0时,在时刻t≥t0的平均状态.这里求和的范围可以取所有允许的状态.如果把函数P拓展定义到全空间RN×[t,+∞),其中当X(t)=x是不允许出现的状态时,定义P(x,t|x,t)=0,000则上面的求和可以拓展到全空间RN:XE(t|x0,t0)=xP(x,t|x0,t0).(1.3.8){eq:1.3:2}x∈RN为简单,我们省略初始条件,而记E(t)=(hX1(t)i,···,hXN(t)i).由(1.3.7),可以得到dhXiXMXi=vjiaj(x)P(x,t)(1.3.9){eq:1.3:3}dtj=1x∈RN方程(1.3.9)给出了平均动力学方程.记Xhaj(X)i=aj(x)P(x,t)x∈RN则XMdhXii=vjihaj(X)i(1.3.10){eq:1.3:4}dtj=1如果系统的随机性可以忽略,则haj(X)i=aj(hXi),我们可以得到化学反应速率方程(1.2.2).但是,一般地,(1.3.10)并不等价于方程(1.2.2).如果是单分子反应,则aj(X)=cjXjk,则haj(X)i=aj(hXi).方程(1.3.10)等价于(1.2.2).因此,我们看到只对单分子反应,常微分方程(1.2.2)才反应平均动力学.因此,对于多分子反应,常微分方程模型(1.2.2)的结果的解释需要特别的小心.如果分子的数目很大,即x≫vj,我们可以把方程(1.3.7)的右边展开成x的泰勒级数.由此,可以得到(我们在这里省略初始条件)∂XMXN∂1X∂2P(x,t)=aj(x)P(x,t)−aj(x)P(x,t)vji+aj(x)P(x,t)∂t∂xi2∂xi∂xkj=1i=11≤i,k≤NXM−aj(x)P(x,t)+···j=13\n系统生物学数学基础其中忽略了vji的高阶项.令XMXMAi(x)=vjiaj(x),Bik=vjivjkaj(x)(1.3.11){eq:1.3:5}j=1j=1则有方程∂XN∂1X∂2P(x,t)=−Ai(x)P(x,t)+Bik(x)P(x,t).(1.3.12){eq:fk}∂t∂xi2∂xi∂xki=11≤i,k≤N这个就是Fokker-Plank方程.其中Ai(x)和Bik(x)的含义我们将在后面介绍.从上面的推导可以看到,Fokker-Plank方程是当分子数很大时,对化学主方程的近似.为描述化学反应的随机性,我们可以分析协方差σik=h(Xi−hXii)(Xk−hXki)i(1≤i,k≤N).(1.3.13){eq:var}协方差σik时时间t的函数.通过概率转移函数P,上面协方差可以表示为Xσik(t)=(xi−hXi(t)i)(xk−hXk(t)i)P(x,t).x∈RN我们可以推导出σik满足的方程:dσikXdhXiiXdhXki=(−)(xk−hXki)P(x,t)+(−)(xi−hXii)P(x,t)dtdtdtx∈RNx∈RNX∂+(xi−hXi(t)i)(xk−hXk(t)i)P(x,t)∂tx∈RNXXM=(xi−hXii)(xk−hXki)(aj(x−vj)P(x−vj,t)−aj(x)P(x,t))x∈RNj=1XMX=(xi−hXii)(xk−hXki)aj(x−vj)P(x−vj,t)j=1x∈RNXMX−(xi−hXii)(xk−hXki)aj(x)P(x,t)j=1x∈RNXMX=(xi+vji−hXii)(xk+vjk−hXki)aj(x)P(x,t)j=1x∈RNXMX−(xi−hXii)(xk−hXki)aj(x)P(x,t)j=1x∈RNXX=[Ai(x)(xk−hXki)+Ak(x)(xi−hXii)]P(x,t)+Bik(x)P(x,t).x∈RNx∈RN这里Ai(x),Bik(x)如前面所定义.由此,我们得到σik所满足的方程.dσikXX=[Ai(x)(xk−hXki)+Ak(x)(xi−hXii)]P(x,t)+Bik(x)P(x,t)(1.3.14){eq:var:1}dtx∈RNx∈RN4\n系统生物学数学基础当随机性很小时,即xi−hXii很小,我们可以把Ai(x)和Bik(x)展开成泰勒级数:XN∂Ai(hXi)Ai(x)=Ai(hXi)+(xl−hXli)+···∂xll=1XN∂Ak(hXi)Ak(x)=Ak(hXi)+(xl−hXli)+···∂xll=1XN∂Bik(hXi)Bik(x)=Bik(hXi)+(xl−hXli)∂xll=1X∂2B(hXi)ik+(xp−hXpi)(xq−hXqi)+···∂xp∂xq1≤p,q≤N带入方程(1.3.14),并且注意到关系XXx∈RN(x−hXi)P(x,t)=0,x∈RNP(x,t)=1ii我们可以得到方程NdσX∂A(hXi)∂A(hXi)X∂2B(hXi)ikikik=σil+σlk+Bik(hXi)+σpq.(1.3.15){eq:var:2}dt∂xl∂xl∂xp∂xql=11≤p,q≤N定义矩阵∂2B(hXi)ikσ=(σik),A=(∂Ai(hXi)/∂xl),B=(Bik(hXi),C=∂xp∂xq上面方程可以简写为dσT=(Aσ+Aσ+Cσ)+B.(1.3.16){eq:df}dt这个就是所谓的Fluctuation-DissipationTheorem(通常的形式是忽略了高阶导数项C).§1.4化学朗之万方程现在我们建立数学模型,用于直接描述随机过程X(t)本身.假设在时刻t,系统的状态为X=xt.令Kj(xt,τ)(τ>0)表示反应Rj在下个时间区间[t,t+τ]内发生的次数.因为每次这样的反应都把分子Si的个数增加vji,系统中分子Si在时刻t+τ的个数为XMXi(t+τ)=xt,i+Kj(xt,τ)vji,(i=1,···,N).(1.4.17){eq:1.4:1}j=1这里,Kj(xt,τ)是随机变量.要得到对所有τ的精确描述,我们需要求解化学主方程.然而,我们可以在下面的条件下给出很好的近似.条件一:首先,取τ充分小,使得在时间区间[t,t+τ]内,系统的状态只有微小的改变,因此,所有的propensityfunction几乎保持不变:′′aj(X(t))≈aj(xt),∀t∈[t,t+τ],∀j∈[1,M].(1.4.18){eq:1.4:2}通常地,每次反应都只使某种分子地个数增加或减少1,所以,当系统地反应物地数量远大于1时,只要取τ充分小,上面的条件一是很容易满足的.根据条件一,在时间区间[t,t+τ]内发生的所有反应都不改变系统的propensityfunction.因此,所有反应在时间区间[t,t+τ]发生的概率可以认为是相互独立的.因此,Kj(xt,τ)等于当propensityfunc-tion等于aj(xt)时,反应通道Rj在时间τ内的发生次数.这个次数满足独立Possion分布Pj(aj(xt),τ).5\n系统生物学数学基础这里P(a,t)表示当某个事件在任意无穷小事件区间dt内出现的概率为adt是,在长度为t的时间区间内出现的次数.令Q(n;a,t)表示P(a,t)等于n(整数)的概率,则由关系Q(0;a,t+dt)=Q(0;a,t)×(1−adt)可以得到∂Q(0;a,t)=−at,Q(0;a,0)=1.∂t由此容易得到Q(0;a,t)=e−at.对任意n≥1,根据概率的乘法定律,由关系ZtQ(n;a,t)=Q(n−1;a,t′)×adt′×Q(0;a,t−t′).t′=0通过数学归纳法,可以得到一般的公式e−at(at)nQ(n;a,t)=,(n=0,1,2,···).n!现在,我们可以计算随机变量P(a,t)的均值和方差hP(a,t)i=var{P(a,t)}=at.当at≫1时,可以证明e−at(at)n(n−at)2≈(2πat)−1/2exp−.n!2at因此,当at≫1,随机变量P(a,t)可以由具有相同的均值和方差的正态分布来近似:P(a,t)≈N(at,at),ifat≫1.(1.4.19){eq:1.A3}因此由条件一,方程(1.4.17)可以近似为XMxi(t+τ)=xt,i+vjiPj(aj(xt),τ),(i=1,···,N).(1.4.20){eq:1.4:3}j=1条件二:时间区间τ充分大,使得在时间区间[t,t+τ]内发生的化学反应的次数的期望值大于1,即hPj(aj(xt),τ)i=aj(xt)τ≫1,∀j∈[1,M].(1.4.21){eq:1.4:4}很显然,这个条件和条件一是矛盾的,可能会出现这样的情况:两个条件无法同时满足.在这种情况下,我们的模型不能满足.但是,在很多情况下,这这两个条件是可以同时满足的,例如,当发生系统中的反应每种分子的个数都足够大时.这时aj(xt)是大数,即使τ很小,上面的条件也是可以满足的.当条件二满足时,我们可以把Possion分布Pj(aj(xt),τ)近似为具有相同的均值和方差的正则分布.因此,我们由下面的关系XMxi(t+τ)=xt,j+vjiNj(aj(xt),τ),(i=1,···,N).(1.4.22){eq:1.4:5}j=1这里N(m,σ2)表示均值为m,方差为σ2的正则分布.注意到在这里,我们把整数的Possion分布变成为连续实数的正则分布.这样,分子数Xi也相应的变成为是实数的.另外,M个正则分布是相互独立的.这是因为我们假定所有的Possion分布Pj都是相互独立的.利用正则分布的简单关系2N(m,σ)=m+σN(0,1),我们可以把(1.4.22)改写为另外的形式:XMXMx(t+τ)=x+va(x)τ+v[a(x)τ]1/2N(0,1),(j=1,···,N).(1.4.23){eq:1.4:6}jt,jjijtjijtjj=1j=16\n系统生物学数学基础这里的正则分布Nj(0,1)都是独立的.下面,假设τ同时满足条件一和条件二,并记τ为dt.另外,我我们用白噪声ξj(t)记满足独立正则分布Nj(0,1)的随机变量.这里,白噪声满足关系hξ(t)i=0,hξ(t)ξ(t′)i=δδ(t−t′),∀i,j∈[1,M],∀t.jijij并记Xj(t)=xt,j,则方程(1.4.23)可以表述为XMXM1/2(x)ξ(t)(dt)1/2,(j=1,···,N).(1.4.24)xi(t+dt)=xi(t)+vjiaj(xt)dt+vjiajtj{eq:1.4:7}j=1j=1引进维纳过程(Winerprocess)Wj,使得dW=W(t+dt)−W(t)=ξ(t)(dt)1/2jjjj可以改写上面的方程为XMXM1/2dxi(t)=vjiaj(xt)dt+vjiaj(xt)dWj,(j=1,···,N).(1.4.25){eq:1.4:8}j=1j=1这里dxj(t)=xj(t+dt)−xj(t).这个就是化学朗之万方程(ChemicalLangevinEquation).§1.5计算模拟前面我们介绍了描述化学反应的几种数学模型,分别涉及到常微分方程,差分方程(化学主方程),偏微分方程(Fokker-Plank方程),随机微分方程.这些方程的计算模拟分别涉及不同的数学领域,可以参考相应的数学专业教材.这里简单介绍如下.§1.5.1常常常微微微分分分方方方程程程的的的数数数值值值模模模拟拟拟差分法,软件:xppaut.§1.5.2求求求解解解化化化学学学主主主方方方程程程Gilliespie算法:1.初始化Xi，并令初始时间t=0.PM2.计算aν(ν=1,···,M)，并令a0=ν=1aν.Pµ−13.产生[0,P1]上的平均分布随机数r1和r2,并令τ=(1/a0)ln(1/r1)，取µ为满足条件ν=1<µr2a0≤ν=1aν的整数，则(τ,µ)是满足概率密度为ae−a0τ,if0≤τ<∞andµ=1,···,MµP(τ,µ)=0otherwise得随机数。4.令t=t+τ,并根据反应反应通道Rµ更新分子个数，即Xi→Xi+vµi.HereP(τ,µ)isthe“reactionprobabilitydensityfunction”thatdefinedasP(τ,µ)dτ=probabilitythat,giventhestate(X1,···,XN)attimet,thenextreactioninVwilloccurinthein-finitesimaltimeinterval(t+τ,t+τ+dτ),andwillbeanRµraction.7\n系统生物学数学基础TheprobabilityPistheproductofP0(τ),theprobabilitythat,giventhestate(X1,···,XN)attimet,noreactionwilloccurinthetimeinterval(t,t+τ);timesaµdτ,thesubsequenceprobabilitythatanRµreactionwilloccurinthetimeinterval(t+τ,t+τ+dτ):P(τ,µ)dτ=P0(τ)·aµdτ.PTofindandexpressionforP(τ),wefirstnotethat[1−adτ′]istheprobabilitythatno0ννreactionwilloccurintimedτ′fromthestate(X,···,X).Therefore,1NXP(τ′+dτ′)=P(τ′)·[1−adτ′]00ννfromwhichitisreadilydeducedthatXMP0(τ)=exp[−aντ],ν=1Thus,weobtainthereactionprobabilitydensityfunctionae−a0τ,if0≤τ<∞andµ=1,···,MµP(τ,µ)=(1.5.26){eq:1.5:1}0otherwise§1.5.3求求求解解解Fokker-Plank方方方程程程差分法.§1.5.4求求求解解解化化化学学学朗朗朗之之之万万万方方方程程程随机微分方程数值方法:以Wt表示随机过程,如果满足下面条件:1.连续;2.独立增量过程:如果t11sec.onoff(1.6.34)给出了在一段时间内(例如,大于1sec),自由的基因位点占总数的百分比与抑制子X的浓度的关系.当X=Kd时,%50的位点是自由的.假设当位点是自由的时,相应基因的转录率为β.则mRNA的产生率(成为promoteractivity)和抑制子X的关系为βpromoteractivity=.(1.6.35){eq:mm4}1+X/Kd这里Kd称为repressioncoefficient.现在,我们考虑另一种情况,X要和诱导物SX结合称复合体[SXX]以后才有活性.如果每个X只能结合一个SX,则有关系X+[SXX]=XT其中XT表示X的总数.假定X和SX碰撞并结合成复合体的速率常数为jon,复合物的分解速率常数为joff.则复合物的动力学方程(质量作用定理)为d[SXX]=konXSX−joff[SXX].(1.6.36){eq:m21}dt9\n系统生物学数学基础并且假定诱导物SX的数量充分大,其数量的变化可以忽略.在平衡态时,有关系KX[SXX]=XSX这里KX是复合物[SXX]的解离常数.则复合物[SXX]的数量和诱导物的数量SX的关系可以由Michaelis-Menten方程(也称为米氏方程)表示出来:XTSX[XSX]=.(1.6.37){eq:mm11}SX+KX现在如果X上有n个结合位点,可以同时和n个SX结合成复合体[nSXX]并且被激活.则有关系[nSXX]+X0=XT.(1.6.38){eq:mm16}这里X0表示自由的X,并且中间态(与X结合的诱导物的个数少于n个)都忽略.复合物[nSXX]的形成是通过X和n个SX分子的碰撞而形成.设反应速率常数为jon,则collisionrate=jXSn.(1.6.39){eq:mm12}on0X令离解常数为joff:dissociationrate=joff[nSXX].(1.6.40){eq:mm13}参数joff通常对应与X和SX之间连接的化学键的强度.复合物[nSXX]的动力学方程为d[nSXX]n=jonX0SX−joff[nSXX](1.6.41){eq:mm14}dt这里假设细胞内S的数量很大,其数目的变化可以忽略.在平衡态的时候,有关系njoff[nSXX]=jonX0SX.(1.6.42){eq:mm15}由关系(1.6.38),我们有(j/j)[nSX]=(X−[nSX])Sn.offonXTXX由此可以得到结合的X所占的比例[nSX]SnXX=(1.6.43){eq:mm17}XTKn+SnXX其中Kn=j/j.这个就是Hillequation,系数n通常称为是Hill系数(Hillcoefficient).当n>1时,Xoffon通常称为是合作的.没有结合的抑制子X的浓度是X01=.(1.6.44){eq:mm18}XT1+(SX/KX)n现在考虑另外的情况,假设存在诱导物S.抑制子可以和诱导物结合为复合体[XSX].诱导物通过和抑制子结合,阻止抑制子抑制基因的表达,从而诱导基因的表达.此时,抑制子X可以有三种状态:自由的,与DNA位点结合,或者与诱导物结合:XT=X0+[XD]+[nSXX].(1.6.45){eq:mm5}这样,我们得到下面动力学方程:d[XD]=konX0D−koff[XD],(1.6.46)dtd[nSXX]n=jonX0SX−joff[nSXX].(1.6.47)dt这里,我们假设复合物[nSXX]不能与D结合,并且细胞内SX的数量很大,其数目的变化可以忽略.在平衡态时,可以得到关系nKX[nSXX]=X0SX,Kd[XD]=X0D(1.6.48){eq:mm8}10\n系统生物学数学基础这里KX=joff/jon为解离常数(forlacrepressor,KX∼1µM∼1000inducer(IPTG)molecules/cell).由上面关系(1.6.45)和(1.6.48)可以求解出X01=.XT1+D/Kd+Sn/KXX由自由DNA所占比例与自由抑制子的浓度的关系(1.6.34),可以得到promoteractivity(记为f=f(SX))和诱导物SX之间的关系βf=.(1.6.49){eq:mm19}1+(XT/Kd)/(1+fDT/(βKd)+Sn/KX)X当Sn/K≫D/K时,上面关系可以近似为XXTdβf=.(1.6.50){eq:mm20}1+(XT/Kd)/(1+Sn/KX)X上面关系给出了基因的活性与诱导物的浓度之间的关系.当SX=0时,有f(SX=0)≈β/(1+XT/Kd).这个也称为是是basalpromoteractivity,表示当没有诱导物时的promoter活性.当S=S∼(X/K)1/nKX1/2TdX时,基因的活性恢复达到最大值的一半(f=β/2).现在考虑激活子的情况:只有当X结合到基因位点D上时,相应的mRNA才会被转录.根据前面的讨论,基因的promoteractivity和激活子的浓度的关系可以通过Michaelis-Menten方程表示出来:βX∗promoteractivity=.(1.6.51){eq:mm10}Kd+X∗这里X∗表示具有活性(可以和DNA位点结合)的激活子的浓度.如果存在诱导物SX可以和激活子结合(假设激活子存在n个作用位点),激活子只有当与n个诱导物结合为复合体[nSXX]后,才有活性(这里忽略中间状态,即与X结合的诱导物的个数少于n个的情况).此时,有关系X+nS⇆X∗,X∗+D⇆D∗X通前面的讨论,有活性的激活子的浓度为XSnX∗=[nSX]=TX.XKn+SnXX因此,基因的活性和诱导物的浓度的关系为βX∗f(SX)=.(1.6.52){eq:mm22}Kd+X∗当S=S=(K/X)1/nKX1/2dTX时,基因的活性达到其最大值的一半.一般地,如果一个基因既有抑制子,又有激活子,基因的活性系数为Pβ(X/K)niiPiiif(X1,···,Xm)=(1.6.53){eq:mm22}1+(Xi/Ki)mii这里Xi表示抑制子或者是激活子的浓度,Ki表示相应的抑制或激活系数.11\n系统生物学数学基础§1.7补充阅读材料1.vanKampen,N.G.1992.Stochasticprocessinphysicsandchemistry.North-Holland,Amster-dam,1992.2.Gillespie,D.T.1977.Exactstochasticsimulationofcoupledchemicalreactions.J.Phys.Chem.81:2340-2361.3.Gillespie,D.T.2000.ThechemicalLangevinequation.J.Chem.Phys.113:297–306.12\n第二章基因表达§2.1引言我们知道,所有生物的遗传信息,都是以基因的形式储藏在细胞内的DNA(或RNA)分子中.随着个体的发育,DNA分子能有序地将其所承载的遗传信息,通过密码子-反密码子系统,转变成蛋白质分子,执行各种生理生化功能,完成生命的全过程.科学家把这个从DNA到蛋白质的过程称为基因表达(geneexpression),对这个过程的调控是现阶段分子生物学研究的中心课题.图2.1:中心法则{fig:dogma}基因表达调控主要表现在以下几个方面:1.转录水平上的调控(transcriptionalregulation);2.mRNA加工成熟水平上的调控(differentialprocessingofRNAtranscript);3.翻译水平上的调控(differentialtranslationofmRNA).§2.2实验事实§2.3数学模型常微分方程模型:13\n系统生物学数学基础PropertyE.coliYeast(S.cerevisae)Proteins/cell4×106∼4×109Timetotran-∼1min∼1minscribeageneTimetotrans-∼2min∼2minlateaproteinTypicalmRNA2−5min10mintoover1hlifetime∼30min(richCellgeneration∼2h(richmediummediumtoseveraltimetoseveralhourshoursTimescaleoftranscription∼1secfactorbindingtoDNAsite表2.1:TypicalparametervaluesfortheBacterialE.ColicellandSaccharonmycescerevisae(Yest)(Alon,2007){tab:1}dX1+−=λ1(n−X1)−λ1X1dtdX2=λ2X1−δ2X2dtdX3=λ3X2−δ3X3dt化学主方程:dP(X1,X2,X3)++=λ1(n−X1+1)P(X1−1,X2,X3)−λ1(n−X1)P(X1,X2,X3)dt−−+λ1(X1+1)P(X1+1,X2,X3)−λ1X1P(X1,X2,X3)+λ2X1P(X1,X2−1,X3)−λ2X1P(X1,X2,X3)+δ2(X2+1)P(X1,X2+1,X3)−δ2X2P(X1,X2,X3)+λ3X2P(X1,X2,X3−1)−λ3X2P(X1,X2,X3)+δ3(X3+1)P(X1,X2,X3+1)−δ3X3P(X1,X2,X3).(0≤X1≤n,X2,X3≥0)ChemicalLangevinequation14\n系统生物学数学基础图2.2:Intrinsicandextrinsicnoiseingeneexpression(Elowitzet.al.2002)dX1+−=λ1(n−X1)−λ1X1(2.3.1)dtqq+−+λ1(n−X1)ξ1(t)−λ1X1ξ2(t)+fλ+(n−X1)ηλ+(t)−fλ−X1ηλ−(t),1111dX2pp=λ2X1−δ2X2+λ2X1ξ3(t)−δ2X2ξ4(t)(2.3.2)dt+fλ2X1ηλ2(t)−fδ2X2ηδ2(t),dX3pp=λ3X2−δ3X3+λ3X2ξ5(t)−δ3X3ξ6(t)(2.3.3)dt+fλ3X2ηλ3(t)−fδ3X3ηδ3(t),图2.3:Amodeloftheexpressionofasinglegene.Eachsteprepresentsseveralbiochemicalreactions,whichareassociatedwithtransitionbetweenpromoterstates,productionanddecayofmRNAsandproteins.{fig:gene}15\n系统生物学数学基础60’md-1.dat’using1:3504030mRNA20100020406080100120140Time图2.4:模拟结果(Gillespie算法):[mRNA]vs.Time.{fig:ge1}60’md.dat’using1:3504030mRNA20100050100150200Time图2.5:模拟结果(求解Langevin方程):[mRNA]vs.Time.{eq:ge2}§2.4基因表达的随机性§2.5反馈控制§2.6补充阅读材料1.Orphanides,G.,Reinberg,D.,(2002)Aunifiedtheoryofgeneexpression,Cell108,439-451.2.Smolen,P.,Baxter,D.A.,Byrne,J.H.,(2000)Mathematicalmodelingofgenenetworks.Neuron26,567-580.3.Kærn,M.,Elston,T.C.,Blake,W.J.,Collins,J.J.,(2005)Stochasticityingeneexpression:fromtheoriestophenotypes.Nat.Rev.Genet.6,451-464.4.Paulsson,J.,(2005)Modelsofstochasticgeneexpression.Phy.LifeRev.2,157-175.16\n系统生物学数学基础5.Elowitz,M.B.,Levine,A.J.,Siggia,E.D.,Swain,P.S.,Stochasticgeneexpressioninasignlecell.Science297(2002),1183-1186.6.Swain,P.S.,Elowitz,M.B.,Siggia,E.D.,Intinsicandextrinsiccontributionstostochasticityingeneexpression.PNAS99(2002),12795-12800.17\n第三章基因调控§3.1ToggleSwitches§3.1.1Bistability正反馈可以产生双稳态.在这里,我们考虑下面的例子.White-opaqueswitchingisanepigeneticphenomenon,wheregeneticallyidenticalcellscanexistintwodistinctivecelltypes,whiteandopaque.Eachcelltypeisstablyinheritedformanygenerations,andswitchingbetweenthetwotypesofcellsoccursstochasticallyandrarely–roughlyoneswitchin104celldivisions.ThegeneWor1wasidentifiedasamaserregulatorofwhite-opaqueswitching[39,28].Inopaquecells,Wor1formsapositivefeedbackloop:itbindsitsownDNAregulatoryregionandactivatesitsowntranscriptionleadingtotheaccumulationofhighlevelsofWor1.上述的正反馈调控过程可以用图3.1表示.图3.1:基因表达的正反馈调控.{fig:3:bistability以GA表示被激活基因的个数(对单拷贝基因,也可以理解为被激活的概率),M表示细胞内相应基因的mRNA的个数,P细胞内基因所表达出来的蛋白质的个数.在很多情况下,蛋白质以二聚物或者高聚物的形式在表示调控基因的表达(与启动子结合).我们以Pn表示n-聚物的个数.这里,我们忽略中间产物,即蛋白质或者以分离的形式存在,或者以n-聚物的形式存在.则上述过程包括下列反应.其中蛋白质的聚合和与启动子的结合是快过程k+λ+1nP⇄Pn,Pn+GR⇄GA.(3.1.1){eq:3:fast1}k−λ−1基因的转录和mRNA的翻译是慢过程+−λ2λ2δ2λ3δ3GA−−→M,GR−−→M,M−→∅,M−→P−→∅.(3.1.2){eq:3:slow}18\n系统生物学数学基础上述过程可以用常微分方程描述为(这里只考虑单基因拷贝,因此GR=1−GA)dPn+n−=kP−kPn(3.1.3)dtdGA+−=λ1Pn(1−GA)−λ1GA(3.1.4)dtdM+−=λ2GA+λ2(1−GA)−δ2M(3.1.5)dtdP+n−=λ3M−δ3P−kP+knPn(3.1.6)dt假设快过程很快达到平衡,即上面方程(3.1.3)-(3.1.4)右边为零,则有PnP=KPn,G=(3.1.7){eq:3:5}nAAn+Pn其中K=k−/k+为n-聚物的解离常数,An=(λ−/λ+)/K.代入慢过程的方程,可以得到11dMPndt=λ21+aAn+Pn−δ2M(3.1.8)dP=λ3M−δ3P(3.1.9)dt−+−−这里λ2=λ2,a=(λ2−λ2)/λ2.令x=M/(δ3A/λ3),y=P/A,t˜=δ3t可以把上面方程组无量纲化dxnn=λ(1+ay/(1+y))−δx(3.1.10)dtdy=x−y(3.1.11)dt这里还以t表示无量纲化以后的时间,并且λ2λ3δ2λ=,δ=.Aδ2δ33由上面方程,系统的平衡态由代数方程g(y)=δy其中yng(y)=λ1+a1+yn所给出.当n=1时,上述方程有唯一的正解p−δ+aλ+λ+(−δ+aλ+λ)2+4δλy∗=.2δ当n>1时,上述方程可以有一个正解,或者三个,也可能有两个.解的个数与参数δ有关.存在临界值δ1<δ2,当δ=δ1或者δ=δ2时,系统有两个平衡点.而当δ<δ1或者δ>δ2时,系统有只有一个平衡点.当δ1<δ<δ2时,系统有三个平衡点.令(x∗,y∗)为某平衡点,该平衡点的稳定性由线性化矩阵−δλg′(y∗)A=1−119\n系统生物学数学基础y*gHyL152.01.5101.050.5y∆0.51.01.52.02.54.55.05.56.06.5(A)(B)图3.2:分岔图(λ=4,a=4,n=4){fig:3:bistability2的特征值所决定.矩阵A的特征值为1pλ1,2=−(1+δ)±(δ+1)2+4(λg′(y∗)−δ).2当特征值满足Re(λ)<0时,平衡点是稳定的.由此,容易得到以下结论:如果λg′(y∗)>δ,则平衡点1,2是不稳定的,如果λg′(y∗)<δ,则平衡点是不稳定的.当δ<δ时,中间值的平衡点是不稳定,另外两12个平衡点(分别对应与高表达水平和低表达水平).当δ<δ1或者δ>δ2时,系统只有一个平衡点,而且是稳定的.因此,当δ1<δ<δ2时,系统表现出双稳态.当δ从小于δ2变为大于δ2时,系统从低表达态变换到高表达态,从δ从大于δ1变为小于δ1时,系统从高表达态变换到低表达态.这样,双稳态为系统的switch创造了条件.§3.1.2Amodelforrepressorexpression[21]Inthecontextofthelysis-lysogenypathwayintheλvirus,theautoregulationofλrepressorexpressioniswellcharacterized.Inthesection,wepresenttwomodelsdescribingtheregulationofsuchanetwork.ThefullpromoterregioninλphagecontainsthethreeoperatorsitesknownasOR1,OR2andOR3.WefirstconsideramutantsystemwherebytheoperatorsiteOR1isabsentfromtheregion.Thebasicdynamicalpropertiesareasfollows:ThegenecIexpressesrepressor(CI),whichinturndimerizesandbindstotheDNAasatranscriptionfactor.Inthemutantsystem,thisbindingcantakeplaceatoneofthetwobindingsites,OR2orOR3.BindingatOR2enhancestranscription,whichtakesplacedownstreamofOR3,whereasbindingatOR3repressestranscription,effectivelyturningoffproduction.图3.3:DynamicalproperteisofλrepressorcI.{fig:sw1}这一调控系统的化学反应可以分为两类:快反应和慢反应.快反应包括分子的结合与分解,相应的反应常数的大概在几秒钟.因此,相对于慢反应(大概需要几分钟),可以认为快反应总是处于平衡态.20\n系统生物学数学基础令X,X2和D分别表示抑制子单体,抑制子dimer和DNA的promotersite的浓度,我们可以把相应的化学平衡方程表示为K12X⇆X2K2D+X2⇆DX2(3.1.12){eq:sw1}K3D+X⇆DX∗22K4DX2+X2⇆DX2X2其中DX和DX∗分别表示dimer结合到位点OR2和OR3上的情况,DXX表示同时结合到两个2222位点上的情况.这里Ki表示平衡常数,并且记K3=σ1K2,K4=σ2K2.则σ1和σ2分别表示相对于dimer-OR2的结合强度.mRNA的转录和降解一般是慢过程:ktDX2+P−→DX2+P+nX(3.1.13){eq:sw2}kdX−→A,这里P表示RNA聚合酶的浓度,n表示平均每个mRNA可以表达的蛋白质的个数.这些反应都是不可逆的.令X=[X],Y=[X],D=[D],U=[DX],V=[DX∗],Z=[DXX]分别为各反应物的浓度,22222则抑制子的浓度的变化方程为dX2=−2k1X+2k−1Y+nktP0U−kdX+r.(3.1.14){eq:sw3}dT这里,我们假设RNA聚合酶的浓度P0保持为常数.参数r表示蛋白质CI的basalrateofproduction,即在没有转录调控因子时基因cI的表达率.另外,变量Y,U和D与X的关系可以通过令快反应(3.1.12)达到平衡来给出:Y=KX21U=KDY=KKDX22122(3.1.15){eq:sw4}V=σ1K2DY=σ1K1K2DXZ=σKUY=σ(KK)2DX422212另外,DNA的总浓度是常数,记为dT:2224DT=D+U+V+Z=D(1+(1+σ1)K1K2X+σ2K1K2X).(3.1.16){eq:sw5}由(3.1.15)-(3.1.16)求解出Y,U,并且代入(3.1.14),得到下面方程(这里注意到关系K1=k1/k−1):dXnkKKPDX2t120TdT=1+(1+σ)KKX2+σK2K2X4−kdX+r.(3.1.17){eq:sw6}112212为简化分析,我们首先把上述方程无量纲化.分别以M和T表示浓度和时间的量纲,则所有浓度的量的量纲都是M,σ1,σ2是无量纲参数,其他参数的量纲如(??)所给出.−1−1−1−1−1[K1]=[K2]=M,[kt]=MT,[kd]=T,[r]=MT.(3.1.18){eq:sw7}√√因此,引进信的无量纲变量x=XK1K2和t=T(rK1K2),可以得到下面的无量纲化方程αx2x˙=−γx+1.(3.1.19){eq:sw8}1+(1+σ1)x2+σ2x4这里的导数表示对无量纲化时间t的导数.并且得到新的无量纲化参数pα=nktP0DT/r,γ=kd/(rK1K2).这里α表示抑制子所能提高的转录率和基础转录率之间的相对比值,γ表示蛋白质的降解率和基础表达率之间的比值.21\n系统生物学数学基础对于λphage的情况,我们有σ1∼1和σ2∼5.因此方程(3.1.19)中有两个参数α和γ.这两个参数可以决定平衡态时抑制子的浓度.下面我们来详细分析.当参数α和γ变化时,系统可以有一个平衡态,或者是三个平衡态.令αx2f(x)=+11+(1+σ1)x2+σ2x4则平衡态的数目取决于方程f(x)=γx的解的个数.对给定的α,存在两个临界值γ1<γ2,使得当γ<γ1或者γ>γ2时,只有一个平衡态;当γ1<γ<γ2时,有三个平衡态;当γ=γ1或者γ=γ2时,有两个平衡态.fHxL2015105x0.20.40.60.81.0图3.4:Bifurcationplotsforthevariablexandconcentrationofλrepressor(α=50,σ1=1,σ2=5){fig:sw2}如果x=x∗是系统的平衡点,则令y=x−x∗,在平衡点附近,我们有线性化方程′∗y˙=(f(x)−γ)y.根据微分方程的稳定性理论,如果f′(x∗)−γ<0,则对应的平衡点是稳定的,如果f′(x∗)−γ>0,则对应的平衡点是不稳定的.特别地,当γ<γ<γ时,系统有三个平衡点(x∗20).(B):随机系统(3.1.20)(σ=0.3).(C):随机系统(3.1.21)(σ=5).(D):随机系统(3.1.23)(σ=0.5).{fig:sw-sim}§3.1.4Noiseinduceswitches–intrinsicnoise上面的例子介绍了外部噪声诱导switches的例子,下面介绍内部噪声诱导switches的例子([26]).相互抑制的基因调控网络可以用化学反应方程描述为[A˙]=gA(1−[rB])−dA[A]−α0[A](1−[rA])+α1[rA],[B˙]=gB(1−[rA])−dB[B]−α0[B](1−[rB])+α1[rB],(3.1.24){eq:3:brepressor[˙rA]=α0[A](1−[rA])−α1[rA],[˙rB]=α0[B](1−[rB])−α1[rB],23\n系统生物学数学基础图3.6:Mutualrepressioncircuit.{fig:brepressor}这里gX,X=A,B为蛋白质X的最大产生率,dX是相应的降解率.为简单期间,我们忽略mRNA阶段,把转录和翻译过程统一看成是蛋白质的合成过程.我们以[rX]表示与蛋白质X结合的基因的相对浓度.则rA表示与蛋白质A结合的基因,控制蛋白质B的合成,rB表示与蛋白质B结合的基因,控制蛋白质A的合成.假定每个基因都是单拷贝的,因此0≤[rX]≤1.参数α0表示蛋白质与promoter的结合率,α1表示离解率.通常地,结合-离解过程相对于其他过程是很快的,因此有α0,α1≫dX,gX.因此,可以把[rX]的相对与时间的导数假定为零,由此可以得到[rX]与[X]的关系,得以下方程[A˙]=gA/(1+k[B])−dA[A],(3.1.25){eq:3:bre1}[B˙]=gB/(1+k[A])−dB[B],这里k=α0/α1是蛋白质得表达强度.由常微分方程描述的上述方程只有一个平衡点.例如,当gA=gB=g和dA=dB=d时,平衡解由p1+4kg/d−1[A]=[B]=2k给出.为了考虑随机效应,我们使用Master方程来描述上面系统.令P(NA,NB,rA,rB)表示细胞在时刻t有NX个自由的蛋白质X和rX个结合的抑制子的概率,这里NX=0,1,2,···,rX=0,1.则上面系统可以用Master方程描述为P˙(NA,N,rA,rB)=gAδrB,0P(NA−1,NB,rA,rB)+gBδrA,0P(NA,NB−1,rA,rB)+dA(NA+1)P(NA+1,NB,rA,rB)+dB(NB+1)P(NA,NB+1,rA,rB)−(gAδrB,0+gBδrA,0)P(NA,NB,rA,rB)−(dANA+dBNB)P(NA,NB,rA,rB)α0[(NA+1)δrA,1P(NA+1,NB,0,rB)+(NB+1)δrB,1P(NA,NB+1,rA,0)]α1[δrA,0P(NA−1,NB,1,rB)+δrB,0P(NA,NB−1,rA,1)]−α0(NAδrA,0+NBδrB,0P(NA,NB,rA,rB)−α1(δrA,1+δrB,1)P(NA,NB,rA,rB).(3.1.26){eq:3:breme}上述方程可以通过Gillespie算法求解.为了分析随机作用的影响,我们可以计算分布XP(NA,NB)=P(NA,NB,rA,rB)rA,rB与参数的关系.在计算模拟中,选取堆成参数g=g=g=0.05(s−1)和d=d=0.005(s−1).并ABAB且,我们比较弱抑制作用(α0=0.005,α1=1.0,k=0.005)和强抑制作用(α0=1.0,α1=0.02,k=50)的情况.计算结果如图3.7所示.计算结果表明,对于若反馈,系统只有一个稳定平衡态,不可能发生转态转移.对于强反馈,系统有三个可能的状态,分别对应于A占优,B占优,和互相抑制.因此,在一定条件下,可以发生状态之间的转移(Fig.3.8).这个就是内部噪声诱导转代转移的机制.24\n系统生物学数学基础(A)(B)图3.7:TheprobabilitiesP(NA,NB)fortheswitch,underconditionof(A)weakrepression(k=0.005)wherethereisonesymmetricpeakand(B)strongrepression(k=50)wherethreepeakappear,onedominatedbyA,theseconddominatedbyB,andthethirdinwhichbothspeciesaremutuallysuppressed.[26]{fig:intrinsicswitch25’md2.dat’using1:22015[A]1050050000100000150000200000250000Time(t)25’md2.dat’using1:32015[B]1050050000100000150000200000250000Time(t)图3.8:ThepopulationofunboundAandBproteinsvstime,obtainfromGillespiesimulationsoftheswitchwithparametersg=0.05,d=0.005,α0=1.0,α2=0.02.[26]{fig:intrinsicswitch225\n系统生物学数学基础§3.2生物振荡§3.2.1AtkinsonOscillatorAtkinsonOscillator的模型如图3.9所示.图3.9:Modulesandconnectivityforthegeneticclock.ThetopconstructcontainsthegInAp2pro-moterfusedtogInG.TranscriptionfromgInAp2requiresthephosphorylatedformoftheenhancerbindingproteinNRI(gInGproduct).ThispromoterisrepressedbyLac1bindingto2perfectlacoperatorsitesO*.ThebottomconstructcontainsthegInKpromoter.ThegInKprmoteralsorequiresNRIPforactivation,howevertheenhancerbindingsitesareles˜spetenthatthosatgInAp2(replottedfrom[5]).{fig:4:atkoscillator令[lacI],[LacI]分别为LacI的mRNA和蛋白质LacI的浓度,[nri],[NRI]分别为NRI的mRNA和蛋白质的浓度,并记磷酸化NRI-P的浓度为[NRI−P],则上述模型可用下面的常微分方程模型描述:d[lacI]=f1([NRI−P])−δ1[lacI]dtd[LacI]=λ2[lacI]−δ2[LacI]dtd[nri]=f2([NRI−P])f3([LacI])−δ3[nri](3.2.27){eq:4:atk1}dtd[NRI]=λ4[nri]−δ4[NRI]−k1[NRI]+k−1[NRI−P]dtd[NRI−P]=k1[NRI]−k−1[NRI−P]−δ5[NRI−P].dt这里λi表示mRNA翻译成蛋白质的反应速率,δi表示分子的降解和稀释的速率常数,k1和k−1表示NRI的磷酸化和去磷酸化的反应速率.函数fi表示蛋白质对基因活性的调控,分别定义如下:([NRI−P]/K)n11f1([NRI−P]=α1,0+α1,11+([NRI−P/K)n11([NRI−P]/K)n22f2([NRI−P])=α2,0+α2,11+([NRI−P/K)n221f3([LacI])=α3,11+(LacI/K)n3.3这里磷酸化过程相对与mRNA的转录和翻译是快过程,可以通过拟平衡假设,即令d[NRI−P]/dt=026\n系统生物学数学基础得到关系k1[NRI−P]=keq[NRI],keq=.k−1+δ5另外,我们引入下面无量纲化变量:[lacI][LacI][nri][NRI]x1=,x2=,x3=,x4=,t˜=δ2t(3.2.28){eq:4:atk2}δ2K3/λ2K3(δ4+δ5keq)(K1/keq)/λ4K1/keq并且引入参数δ1δ3(δ4+δ5keq)β1=,β3=,β4=δ2δ2δ2α1,1α2,1K2α1=,α2=,a=α1,0α2,0K1α1,0λ2λ4keqα2,0α3,1λ1=,λ3=,δ1δ2K3δ3(δ4+δ5keq)K1可以得到下面的无量纲化方程(这里还以t记无量纲化时间)n1dx1x4=β1λ11+α1n1−x1dt1+x4dx2=x1−x2dtn2(3.2.29){eq:4:atk3}dx3(x4/a)1dt=β3λ31+α21+(x/a)n21+xn3−x342dx4=β4(x3−x4).dt图3.10给出了模拟结果.可以看到,在一定条件下,我们得到了周期变化蛋白质表达水平.x2HtLx4HtL352.5302.0251.520151.0100.55tt2040608010020406080100(A)(B)图3.10:Atkinsonoscillator计算模拟结果.这里的无量纲参数为:β1=β3=30.0,β4=1.0,λ1=λ3=2.0,α1=α2=20.0,α3=1.0,a=1.0,n1=4,n2=5,n3=1,xi(0)=0.0.{fig:4:atk1}为了分析周期解的产生机制,我们把上述方程进一步简化.这里mRNA的转录都是快过程,因此上面方程中可以近似认为dx1/dt=dx3/dt=0.由此,可以把方程简化成二阶平面系统n1dx2x4=λ11+α1n1−x2dt1+x4(3.2.30){eq:atk:4}dx(x/a)n2134dt=λ31+α21+(x/a)n21+xn3−x4.42为方便起见,下面去n3=1.上述系统的平衡点由曲线n1x4x2=λ11+α1n11+x427\n系统生物学数学基础和λ(x/a)n234x2=1+α2−1x41+(x4/a)n2的交点给出.根据参数的不同取值,系统可以有一个,两个或者三个平衡点.平衡点的稳定性也和系统的参数有关.特别地,当系统只有一个参数,平且是不稳定时,就会出现稳定的周期解.下面的图3.11给出了几种情况下的解.xx2HtLx4HtL21.440201.2301.0150.820100.60.45100.2x412345tt2040608010020406080100(A)(A’)(A”)x2HtLx4HtLx2352.540302.025301.52020151.0100.5105x412345tt2040608010020406080100(B)(B’)(B”)x2x2HtLx4HtL4.0403.535303.0302.5202.0x4t1.51234520406080100t20406080100(C)(C’)(C”)xx2HtLx4HtL22.0400.201.5300.151.00.10200.50.05x412345tt2040608010020406080100(D)(D’)(D”)图3.11:不同的参数值所对应的动力学行为.(A):n2=4,α2=10.5;(B):α2=20,n2=5;(C)λ3=0.5,α2=100,n1=4,x2(0)=30.0,x4(0)=4.0,x1(0)x3(0)=0;(D):参数值同(C),初值为xi(0)=0.0其他参数同图3.10所给出.{fig:4:atk2}§3.2.2Asyntheticgene-metabolicoscillator基因代谢振子模型如图3.12所示.输入代谢物M1,则酶E1可以把M1转化为M2.随着代谢物M2的增加,酶E1被抑制,另外,酶E2被M2所激活.因为酶E2的作用,代谢物M2转变M1.这样就形成一个循环.这个振子可以通过Escherichiacoli来实现[15](图3.13).SuchaconceptualdesignwasrealizedusingtheacetatepathwayinE.coli(Fig.3.13).TheM1poolisacetylcoenzymeA(acetyl-CoA)andtheM2poolconsistsofacetylphosphate(AcP),acetate(OAc−)andtheprotonatedformacetate(HOAc).Acetyl-CoAisametabolicproductofsugar,fattyacidsandsomeaminoacids,andistheentrypointintothetricarboxylicacid(TCA)cycle(三羧酸循环).Acetyl-CoAisconcertedtoacetylphosphateinE.Colibyphosphateacetyltransferase(Pta),whichcorrespondstoenzymeE1inFig.3.12,andthentoacetatebyacetatekinase(Ack).28\n系统生物学数学基础Theprotonatedformofacetateispermeableacrossthecellmembrane.Underaerobicconditions,acetyl-CoAisfurtheroxidizedintheTCAcycle.TheremainingfluxgoestoproduceeitheracetateintheTCAcycle.Theremainingfluxgoestoproduceeitheracetateorethanol.Inwild-typeE.coli,theenzymeacetyl-CoAsynthetase(Acs)isinducedinthepresenceofacetate.However,suchinductionisundercataboliterepressionbyglucoseinthewild-typestrainsoastoavoidfutilecycling.Inourdesign,acetyl-CoAsynthetaseisusedasenzymeE2inFig.3.12.Thus,bothphosphateacetyltransferaseandacetyl-CoAsynthetaseneedtobere-wiredtorespondtotheM2pool.(A)(B)图3.12:Conceptualdiagramoftheoscillatorydynamics,highlightingthetwometabolitepools(M1andM2)andtheircontrols.Solidlinesindicatemetabolicfluxes.Dashedlinesindicatepositive(arrow)andnegative(bluntbar)transcriptionalortranslationalregulation(ref.[15]).{fig:4:metabolicoscillator−1上述模型可以用下面的常微分方程系统来描述.代谢物(AcCoA,AcP,OAc,HOAc)的浓度的变化满足方程dAcCoA=VAcs−VPta+Vgly−VTCAdtdAcP=VPta−VAckdt−1(3.2.31){eq:4:met1}dOAc=VAck−VAcE−VAcsdtdHOAc=VAcE−Voutdt其中Vi表示化学反应速率,如表3.1所示.这里glycolytic速率用常数表示.TCA循环和EtOH与HOAc的输出流都假设为关于AcCoA和HOAc是一阶的.酶催化流Pta,Ack和Acs都假定为是Michaelis-Menten形式的.另外,三种关键蛋白质LacI,Pta和Acs的浓度的变化由下面方程给出:dLacI=RLacI−Rd,LacIdtdPta=RPta−Rd,Pta(3.2.32){eq:4:met2}dtdAcs=RAcs−Rd,Acsdt这里R表示蛋白质的合成率,Rd表示降解率.这里蛋白质的合成包括mRNA转录和翻译,通过Hill方程描述,而降解(或者因为细胞生长因此的浓度稀释)用一阶动力学描述(表3.1).方程的解如图3.14所示.由计算结果可以看到,代谢物的周期振荡可以由糖醇解率Vgly的增加而激发.当Vgly小时,系统达到平衡态,没有振荡解,而当Vgly大时,系统的平衡态失去稳定性,出现周期振荡解.§3.2.3Mechanismsofnoise-resistanceingeneticoscillators这里介绍一种随机激励振子.振子的模型如图3.15所示.这个模型包含激活子和抑制子.这一基因调控网络包含两个基因,一个激活子A和一个抑制子R.激活子A与A和R的启动子结合,提高相应的基因转率效率.R可以和A结合称为复合物,从而分离激活子A.这样,R起抑制的作用.29\n系统生物学数学基础图3.13:BiologicalrealizationoftheconceptualdesigninFig.3.12.Theyellowboxeshighlightthetwometabolicpools,M1andM2.Ack,acetatekinase(醋酸盐激酶);Acp,acetylphosphate(乙酰磷酸盐);Acs,acetyl-Coasynthetase;OAc−,acetate(乙酰);Pta,phosphateacetyltransferase.(adoptfrom[15]){fig:4:ecolisugarsVgly=0.001Vgly=0.01MetabliteconcentrationsMetabliteconcentrationsAcPAcP0.10.10.0010.001AcCoAAcCoA10-510-510-710-7tt50010001500200025003000100200300400500Vgly=0.05Vgly=0.5MetabliteconcentrationsMetabliteconcentrations10AcP0.11AcP0.10.001AcCoAAcCoA0.0110-50.001tt100200300400500100200300400500图3.14:Computationalcharacterizationofthemetabolator.ThemetabolatorispronetooscillateatincreasingglycolyticratesVgly.Vglyinthefourpanels(fromlefttoright)are0.001,0.01,0.05and0.5.(replotfrom[15]){fig:4:meta2}30\n系统生物学数学基础RateexpressionParametersGlycolyticflux,VglyVgly=S0S0=0.001···0.5FluxtoTCAcycleandVTCA=kTCAAcCoAkTCA=10EtOHproduction,VTCAk1Pta·AcCoAPtaflux,VPtaVPta=k1=80,Km,1=0.06Km,1+AcCoA−1k2Acs·OAcAcsflux,VAcsVAcs=−1k2=0.8,Km,2=0.1Km,2+OAcFluxforthereactionAck−1k=1,k=1−1VAck=kAck,fAcP−kAck,rOAcAck,fAck,rAcPGGGGGGGGBOAcFGGGGGGGGAcid-baseequilibriumforV=C(AcP·H+−KOAc−)C=100,H+=10−7,K=10−4.5AcEeqeqaceticacid,VAcEHOAcintercellularVout=k3(HOAc−HOAcE)k3=0.01,HOAcE=0transportrate,Voutα(AcP/K)n1g,1LacIsynthesisrate,RLacIRLacI=n+α0α1=0.1,Kg,1=10,n=2,α0=01+(AcP/Kg,1)α2(AcP/Kg,2)nKg,2=10,n=2,α0=0,Acssynthesisrate,RAcsRAcs=n+α01+(AcP/Kg,2)α2=2α3Kg,3=0.001,n=2,α0=Patsynthesisrate,RPtaRPta=n+α01+(LacI/Kg,3)0,α3=2DegradationrateRd,X,Rd,X=kdXkd=0.06(X=LacI,Acs,Pta)表3.1:Rateexpressionandvaluesusedforparameters[15]{tab:4:met1}图3.15:Biochemicalnetworkofthecircadianoscillatormodel.(replotfrom[37]){fig:4.circadianoscillator31\n系统生物学数学基础上述系统可以用确定性动力学系统描述为dD/dt=θD′−γDAAAAAAdD/dt=θD′−γDARRRRRdD′/dt=γDA−θD′AAAAAdD′/dt=γDA−θD′RRRRRdM/dt=α′D′+αD−δMAAAAAMAA(3.2.33){eq:4:osc3}dA/dt=βM+θD′+θD′AAAARR−A(γADA+γRDR+γCR+δA)dM/dt=α′D′+αD−δMRRRRRMRRdR/dt=βRMR−γCAR+δAC−δRRdC/dt=γCAR−δAC.D′andDdenotethenumberofactivatorgeneswithandwithoutAAboundtoitspromoterAArespectively;similarly,D′andDreferetotherepressorpromoter;MandMdenotemRNAofRRARAandR;AandRcorrespondtotheactivatorandrepressorproteins;andCcorrespondstotheinactivatedcomplexformedbyAandR.Theconstantsαandα′denotethebasalandactivatedratesoftranscription,βtheratesoftranslation,δtheratesofspontaneousdegradation,γtheratesofbindingofAtoothercomponents,andθdenotestheratesofunbindingofAformthosecomponents.Thecellularvolumeisassumedtobetheunitysothatconcentrationandnumberofmoleculesareequivalent.NoticethatweassumethatthecomplexbreaksintoRbecauseofthedegradationofA,andtherefore,theparameterδAappearstwiceinthemodel.模型的计算模拟结果如图3.16所示.Concentrations15001000500t100200300400500(A)(B)图3.16:Oscillationsinrepressorandactivatorproteinnumbersobtainedfromnumericalsimulationsofthedeterministic(A)andstochastic(B)descriptionsofthemodel.Thevaluesofreac-tionsratesare:α=50h−1,α′=500h−1,α=0.01h−1,α′=50h−1,β=40h−1,β′=AARRAR5h−1,δ=10h−1,δ=0.5h−1,δ=1h−1,δ=0.2h−1,γ=1mol−1h−1,γ=MAMRARAR1mol−1h−1,γ=2mol−1h−1,θ=50h−1,θ=100h−1.TheinitialconditionsareCARD=D=1mol,D′=D=M=M=A=R=C=0,whichrequirethatthecellhasARARARsinglecopyoftheactivatorandrepressorgenes:D+D′=1molandD+D′=1mol.[37]{fig:4.circadianoscillatorAARR为研究上述振荡现象的形成机理,下面来把系统简化.我们所研究的系统包含有快和慢过程.把快慢过程分别处理,并且把快过程用拟平衡态来近似是常用的简化系统的方法.在本系统中,promoter的状态改变(通过promoter与激活子A的结合和解离)是快过程(θ=50h−1,θ=100h−1).基因的AR转录过程也是快过程(α=50h−1,α′=500h−1,α=0.01h−1,α′=50h−1).因此,可以把激活AARR的promoter的个数和mRNA的个数都认为是常数,即方程(3.2.33)中令dDdDdD′dD′dMdMA=R=A=R=A=R=0.dtdtdtdtdtdt32\n系统生物学数学基础由此,可以求解出一下关系:θAθRDA=,DR=θA+γAAθR+γRA′γAA′γRADA=,DR=(3.2.34){eq:4:osc3-equi1θA+γAAθR+γRA1α′γA+αθ1α′γA+αθM=AAAA,M=RRRR.ARδMAθA+γAAδMRθR+γRA另外,由于蛋白质A与promoter的结合是很快的过程,也可以认为A的浓度很快可以达到平衡,即dA=0.dt由此可以得到β1α′γA+αθAAAAA=.(3.2.35){eq:4:osc3-equi2γCR+δAδMAθA+γAA代入上面的关系,可以求解出q1′1′ρ(R)−K)2+4αρ(R)K(3.2.36)A=A˜(R)=(αAρ(R)−Kd)+(αAdAd{eq:4:osc3-equi322其中βAρ(R)=,Kd=θA/γA.δMA(γCR+δA)最后,我们可以得到两个慢变量的自治系统dRβRαRθR+α′γRA˜(R)=R−γA˜(R)R+δC−δRA˜(R)CARdtδMRθR+γR(3.2.37){eq:4:osc3-slow}dC=γCA˜(R)R−δACdt图3.17给出了简化方程的数值解,其中参数值和前面一样,我们可以看到,简化方程的解与原来的方程是一致的,除了每个循环的最大值和周期.因此,我们可以相信这个简化方程可以帮助我们了解循环产生的内部机制.ConcentrationsConcentrations15002000150010001000500500tt5010015020050100150200(A)(B)图3.17:计算模拟结果:(A)完整系统,(B)简化系统.{fig:4.simplify}方程(3.2.37)是一个二阶自治系统.对该方程的定性分析,可以使用常微分方程的通用方法.特别地,因为分子个数不可能是负的,也不可能是无穷大,上述方程的解一定是有界的(思考:如果在数学上严格证明这个结论?).另外,二阶自治系统的解曲线在相平面内不可能和自己相交.因此,改系统的解或者收敛到周期解,或者收敛高不动点.系统的不动点可以由下面条件给出:dR/dt=dC/dt=0,即求解方程βαθ+α′γA˜(R)RRRRR=δRR,C=(γC/δA)A˜(R)R.δMRθR+γRA˜(R)33\n系统生物学数学基础在上面参数条件下,上面的平衡方程有唯一的正解,记为(R∗,C∗).因此,如果该平衡点不稳定,则系统一定有稳定的极限环,即周期解.在我们所研究的参数条件下,平衡解为(R∗,C∗)=(66.7491,363.47).对应于平衡点的线性化矩阵的特征值为λ1,2=0.405797±0.5565i.因此,特征值有正实部,对应的平衡点是不稳定的.这样,我们证明了在一定条件下,可以出现稳定的振动状态.平衡点的稳定性可以通过平衡点处的线性化矩阵的秩来刻划.特别地,当线性化矩阵的秩大于零时,对应的特征值具有正实部.对上面的例子,平衡点处的线性化矩阵所对应的秩为∂A˜(R∗)β(α′−α)θγτ=RRRRR−γR∗−(γA˜(R∗)+δ+δ).A˜(R))CCRA∂RδMR(θR+γR2当τ>0时,相应的平衡点是不稳定的.对于我们所计算的参数值,有τ=0.811594>0.因此,平衡点是不稳定的.可以证明,保证平衡点不稳定的参数范围是比较大的.例如,如果同时改变θA和θR:(θA,θR)→(KθA,θR)当0.0240(图3.18).所有的转录率(α和α′)都乘K>0.08倍,或者蛋白质和mRNA的降解率都乘0.00090.因此,在这里所讨论的基因调控网络中,允许周期振荡的结果并不限于特殊选择的参数值,而是对比较大范围内选取的参数值都可以满足.系统的动力学行为的向量场分析如图??(A)所示.由此可以分析出振荡出现的内部机制.Τ1.00.5K2468101214-0.5-1.0图3.18:τ与K的依赖关系.{fig:4.ciroscirobust如果平衡点时稳定的,不能确定是否存在稳定的极限环.特别地,在一定条件下,确定性方程只有一个平衡点,因此不可能出现振荡解.但是,当有随机因素时,平衡态可以因为随机因素的干扰而变得不稳定.此时,平衡态可以被激发到失稳的状态,而出现由于随机激发所引起的周期振荡(如图3.20).随机激发出振荡的机制如图3.19(B)所示.§3.2.4振振振荡荡荡的的的同同同步步步前面介绍了几种单细胞基因调控网络引起振荡的例子.然而,在生物系统中,通常需要不同细胞的振子同步振动,例如生命节律的周期振动.这一节介绍振荡同步的机制.振荡的同步可以通过不同的机制实现,例如细胞间的通讯,外界激励等.细胞间的通讯已经在很多文献中介绍过.这一节中,我们介绍一个通过外界激励引起振荡同步的例子.34\n系统生物学数学基础(A)(B)图3.19:简化方程的向量场分析图示([37]).{fig:4:cirophaseR200015001000500t100200300400500(A)(B)图3.20:计算模拟结果:(A)确定性系统,(B)随机模型.这里δR=0.005,其他参数和前面一样.{fig:4:noise-osc35\n系统生物学数学基础考虑图3.21所考虑的例子.在这个例子中,每个细胞包含一个由三个抑制子组成的循环抑制调控网络,其中蛋白质LacI抑制基因tetR的promoter,而蛋白质TetR抑制基因cI的promoter,蛋白质CI由抑制lacI的promoter.这样的机制可以产生周期振动.在我们所考虑的系统中,有N个细胞,每个细胞包含相同的振子(基因调控网络和相应的参数都相同).所以,我们有N个振子.但是,这些振子的相位并不相同.在这个例子中,不同的细胞收到公共的外界刺激(AI),我们将看到,当外界的刺激满足一定的条件时,这些振子可以达到同步振动.图3.21:Theschematicdiagramofasyntheticgeneregulatorynetwork.{fig:5:syn1}以xk,yk,zk分别表示第k个细胞中基因lacI,tetR,cI所转录出来的mRNA的浓度,Xk,Yk,Zk第k个细胞中的蛋白质LacI,TetR和CI的浓度.以A表示外界刺激AI物的浓度.上面的系统可以用下面方程组表示dxkαγkAdXk=−xk+n+,=β(xk−Xk),dt1+Y1+AdtkdykαdYk=−yk+n,=β(yk−Yk),(3.2.38){eq:4:syn1}dt1+ZdtkdzkαdZk=−zk+n,=β(zk−Zk).dt1+Xdtk在这里,蛋白质对promoter的抑制作用都通过Hill函数表示,环境的刺激通过Michaels-Menten函数表示.所有mRNA的降解率都相同,这里通过选取适当的单位,可以设为1.另外,假定所有蛋白质的生成率和降解率都相同.刺激物AI的浓度的变化包括合成和降解两部分.另外,还有细胞外的刺激G(t):dA=λ−kAA+G(t).(3.2.39){eq:5:syn2}dt我们考虑三种形式的刺激:•周期激励:G(t)=σsin(ωt);•高斯白噪声激励:G(t)=σξ(t),满足hξ(t)i=0,hξ(t)ξ(t′)i=δ(t−t′);P∞•周期脉冲激励:G(t)=j=1σjδ(t−tj),其中tj=jτ,τ表示脉冲周期,σj以概率1/2取值σ或者−σ.36\n系统生物学数学基础如果没有外界的刺激,所有振子以相同的频率振动,但是振动的相位不相同(图3.22).为研究同步性,我们引入序参数XN1R(t)=wkexp(iθk(t)),(3.2.40){e33}Nk=1PN这里θk表示第k个振子的相位,ωk=|θ˙k(t)/n=1|θ˙k(t)|表示第k个振子的权重,i为虚数单位.这里引入权重因子wk是因为在数据采样时,因为振子在不同的相位具有不同的速度,因此数据并不是在所有相位上均匀的.根据上面的定义,R=1表示N个振子完全同步,R=0表示完全不同步,00且q>0,则平衡点是稳定的.3.如果p>0且q<0,则平衡点是不稳定的.39\n系统生物学数学基础(a)(b)0.910.80.80.70.60.6RR0.50.40.4s=0.70.20.3s=0.8s=0.90.2000.511.520500100015002000t(min)4t(min)x10(c)(d)18a(t)b(t)0.870.66R0.450.240300.20.40.60.81800100012001400160018002000st(min)图3.24:Synchronizationinthepresentofwhitenoisestimulus.(a)showsthetimeevlutionoforderparameterofthesystemunderthewhitenoisestimuluswithσ=0.2.Itsuggestthataweakweakisabletosynchronizatetheoscillators.(b)showsthetimeevlutionoforderparametersunderthewhitenoisestimuluswithdifferentvaluesofthevariance(σ=0.7,0.8,0.9).Itshowsthatthetimeittakesthesytemtoarchivethecompletesynchronizateddependsonthevarianceσ2.Fordifferentvalueofσthatvaryfrom0to1.0,theorderparameterofthesystemattimet=2000isshownat(c).Ingeneral,theorderparameteratafixtimeisincreasewithrespecttoσ.Thedashedlineshowsthatfittingcurvebyfunction(??)witha=6.53,b=5.03.TheorderparameterRσ(t)atanytimetcanbeapproximatedbyfunction(3.2.42),withcoefficientsa(t)andb(t)givenby(d).Whenthetimetislarge,a(t)approachconstantlimitsa≈6.10,andb(t)dependsonthetimetlinearly.Thedashedlinerepresentthefittingcurveofb(t)by(3.2.43).TheparametersusedinthissimulationaregivenatFigure3.22.{fig:4:f5}40\n系统生物学数学基础1s=1000s=20000.90.80.70.6R0.50.40.30.20.1056789101112131415t图3.25:Theorderparameterfor1000cellswithperiodτvariesfrom5to15,anddifferentvaluesofσ(redcirclesforσ=1000andbluesquaresforσ=2000).{fig:4:f7}1SinusoidalperiodicstimulusWhitenoisestimulus0.90.80.70.6R0.50.40.30.20.1000.10.20.30.40.50.60.70.80.91Db图3.26:Therelationbetweentheorderparameterandthevariabilityofthecells(measuredby∆β.Theparameterusedareσ=0.54andω=0.8forthesinusoidalperiodicstimulus.{fig:4:f8}41\n系统生物学数学基础§3.3生命节律节律是生命的普遍现象.很多器官表现出有规律的生理和行为上的以24小时为周期的周期现象.在相同的光照和温度条件下,这些现象的周期大概为一天,与日照时间基本同步.从基因调控的层次上来看,这些现象是和体内控制人体对外界光照的响应的时钟基因PER和TIM联系起来的.在这一节中,我们从基因调控的角度,建立节律控制的基因调控网络的数学模型,并且分析这些模型,讨论节律振荡的机制.§3.3.1DimerizationandproteolysisofPERandTIM首先,我们来讨论果蝇的节律控制的基因调控网络(Fig.3.27).在这个模型中,时钟基因所表达的蛋白质PER和TIM形成稳定的聚合物,调控他们本身的表达.PER单体可以迅速被磷酸化并且降解,而PER/TIM二聚体比较稳定,不容易降解.因此,当PER的数量比较少时,很快被降解.但是,当PER的数量比较多时,倾向于形成二聚体,而不容易降解.这样,PER的降解率与其数量之间满足非线性的关系.这里二聚物的作用相当于对PER的数量有一个正反馈调控.PER二聚体进一步和蛋白质TIM结合形成聚和物.而后,这个聚合物被传输到细胞核内,分别与两个基因的promoter结合,负调控这些基因的表达.图3.27:AsimplemolecularmechanismforthecircadianclockinDrosophila,adoptedfrom[35].PERandTIMproteinsaresynthesizedinthecytoplasm,wheretheymaybedestroyedbyprote-olysisortheymaycombinetoformrelativelystableheterodimers.Heteromericcomplexesaretransportedintothenucleus,wheretheyinhibittranscriptionofperandtimmRNA.WeassumethatPERmonomersarerapidlyphosphorylatedbyDBTandthendegraded.Dimers,weassume,arepoorersubstratedsforDBT.{fig:4:cir1}上面的机制可以用微分方程描述如下.首先,我们作一些简化.首先,PER和TIM的地位是类似的,并且在细胞内的变化也大概是同步的的.因此,把他们看成统一的时钟基因蛋白.另外,我们假定细胞核与细胞质内聚合物的浓度很快达到平衡,而不区分细胞核内外的聚合物浓度之间的差别.这样,我们的模型可以简化为三个变量,分别是mRNA的浓度M,单体的浓度P1,和集合物的浓度P2.模型可42\n系统生物学数学基础以用下面给的常微分方程描述dMvm=−kmM(3.3.46)dt1+(P2/Pcrit)2dP1kp1P12=vpM−−kp3P1−2kaP1+2kdP2(3.3.47)dtJp+P1+rP2dP22kp2P2=kaP1−kdP2−−kp3P2.(3.3.48)dtJp+P1+rP2这里,聚合物对时钟基因的表达的影响是有合作效应的(Hill系数为2).这里,单体和复合体都可以和DBT结合并且降解,但是单体的降解率比复合体的降解率大得多(k′≫k).单体与聚合体除p1p2了通过和DBT结合降解外,还有自己得慢降解过程(通过系数kp3表现出来).在这里,方程(3.3.47)和(3.3.48)中得Michaelis-Menten函数表示DBT催化得降解过程.可以通过下面得过程推导出来.这个过程可以描述为P1+DBT⇆DBT−P1→DBT,P2+DBT⇆DBT−P2→DBT因此,令D表示DBT的浓度,D1表示DBT−P1的浓度,D2表示DBT−P2的浓度,有dP1=−k1P1D+k−1D1(3.3.49)dtdP2=−k2P2D+k−2D2(3.3.50)dtdD1′=k1P1D−k−1D1−k1D1(3.3.51)dtdD2′=k2P2D−k−2D2−k2D2.(3.3.52)dt另外,有关系DT=D+D1+D2.因此,假定降解过程很快,有dD1/dt=dD2/dt=0即可以得到DTD1=keq,1P1D,D2=keq,2P2D,D=,1+keq,1P1+keq,2P2其中k1k2keq,1=′,keq,2=′.k−1+k1k−2+k2这样,我们有dPDk′P1T11=−(3.3.53)dt1/keq,1+P1+(keq,2/keq,1)P2dPDk′(k/k)P2T2eq,2eq,11=−(3.3.54)dt1/keq,1+P1+(keq,2/keq,1)P2比较上面的Michaelis-Menten函数,可以得到关系k′=Dk′,k=Dk′(k/k),J=1/k,r=k/k.p1T1p2T2eq,2eq,1peq,1eq,2eq,1因此,条件k′≫k即p1p2k′k≫kk,1eq,12eq,2也即kk′kk′1122≫.k−1+k1′k−2+k2′43\n系统生物学数学基础mRNA&ProteinLevelmRNA&ProteinLevel3.5Protein3.0Protein3.02.52.52.0mRNA2.01.51.51.0mRNA1.00.50.5TimeHhrLTimeHhrL2040608010020406080100(A)(B)图3.28:数值模拟结果.这里的参数为:v=1,k=0.1,v=0.5,k′=10,k=0.03,k=mmpp1p2p30.1,Pcrit=0.1,Jp=0.05,r=1.2,ka=800,kd=4.(A):完整模型(3.3.46)-(3.3.48).(B):简化模型(3.3.56)-(3.3.57).{fig:4:cirsol1}当k≫k时且k≪k′,k≪k′,即单体更加容易和DBT结合,并且结合都很容易降解,上面条件12−11−22是满足的.这样就证明了上面的方程和假设与模型的一致性.上述模型的数值解如图(3.28)(A)所示.我们看到,在一定的参数条件下,可以产生周期为24小时的节律振荡.为了分析节律振荡的机制,我们进一步简化上述模型.聚合反应相对于基因的表达和示快速化学反应过程(ka和kd都很大),这时,可以假设2P1→P2总能达到平衡.因此,令PT=P1+2P2为蛋白质的总数,则由平衡条件,可以得到关系P2=KeqPT,Keq=ka/kb.记P=qP,P=1(1−q)P,则有1T22T2q=p.(3.3.55){eq:4:cirq}1+1+8KeqPT这样,上面的方程可以简化为下面的两自由度的平面系统dMvm=−kmM(3.3.56)dt1+(PT(1−q)/(2Pcrit))2dPTkp1PTq+kp2PT=vpM−−kp3PT.(3.3.57)dtJp+qPT+(r/2)(1−q)PT这里q有关系(3.3.55)给出,并且k=k′−k≈k′.该系统的数值模拟结果见图3.28(B).p1p1p2p1下面对上述系统进行定性分析.为此,令vmkp1PTq+kp2PTf(PT)=2,g(PT)=+kp3PT.1+(PT(1−q)/(2Pcrit))Jp+qPT+(r/2)(1−q)PT则上述系统的平衡点由方程kmM=f(PT),vpM=g(PT)的解给出,即两条曲线M=f(P)/k和M=g(P)/v的交点.令(M∗,P∗)为响应的平衡点,则在平TmTp衡点处的线性化矩阵为−k−f′(P∗)mA=′∗vp−g(P)当−tr(A)=k+g′(P∗)<0时,对应的平衡点是不稳定的.由此,可以看到,要想得到不稳定的平衡m点,我们需要′∗g(P)<−km(3.3.58){eq:4:cond}即曲线M=g(PT)/vp在平衡点处的切线斜率比较小(小于零).不等式(3.3.58)给出了平衡点不稳定的充分条件.如果平衡点是唯一的(这个条件在很到的参数范围内是成立的),也给出了存在稳定极限环的充分条件.一般地,g′(P∗)依赖于系统的参数,图3.29(A)给出了与K的依赖关系.从这个关系可以eq看到,当参数如图3.28所给出时,如果K>8.3,则有关系g′(P∗)<−k,因此可以得到周期振荡解.eqm当Keq=200时,对应的解在相平面上如图3.29所示.44\n系统生物学数学基础另外,参数kp1的值也影响平衡点的稳定性(图3.30).周期振荡的周期也和参数有关.数值模拟的结果显示,振荡的周期随参数kp1的增加而增加,而与Keq的依赖关系并不明显.周期与参数Keq和kp1的关系如图3.31所示.在实验中,通常会通过基因突变技术或者环境的改变研究基因或者环境因素对节律振荡的影响.在模型研究中,可以通过改变相应参数对节律振荡的影响.例如,当温度改变时,通常会改变化学反应平衡常数Keq.对于某些基因突变,可能会影响蛋白质的降解率kp1,kp2等.图3.32给出了Tyson等人的数值模拟结果,可以看到在某种程度上,可以通过修改模型参数来模拟基因突变的实验结果[35].Mg'HP*L1.040.53Keq501001502002-0.5-1.01-1.5PT0.51.01.52.02.53.0(A)(B)图3.29:(A)g′(P∗)与参数K的依赖关系.(B)相平面图.eq{fig:4:circond}g'HP*LProteinLevel26kp1=30kp1=20154kp13kp1=10102030402-11-2TimeHhrL20406080100(A)(B)图3.30:(A)′g′(P∗)与参数k的依赖关系.(B)解与参数k的关系.这里k分别取值10,20,30,其p1p1p1它参数和前面一致.{fig:4:circond2}§3.3.2Circadianrhythmgenerator这一节介绍一个简单的节律振荡的数学模型.该模型可以看成是很多生物振荡的产生机制的简化版本.节律调控的一般模型如图3.33A所示.这里,参与节律调控的蛋白质在细胞核内被表达出来以后,进入到细胞质.在这里,蛋白质经过一系列的反应,包括磷酸化,生物多聚物,和蛋白质传输等过程,最后生成所谓的有效蛋白,参与调控与节律有关的基因的表达.另一方面,这些有效蛋白会被传输到细胞核内,负调控其基因的表达.因此,这个调控环路包括一个负反馈调控.另一方面,因为蛋白质的传输和生化反应等都需要时间,有时候还是相当长的时间,因此这个负反馈的环路存在时间上的滞后.图3.33B给出了相应的示意图.这里主要的变化量的mRNA和有效蛋白质的浓度.其中,从mRNA到蛋白质的过程有时间的延滞,并且有非线性的关系.这里的非线性关系通常表示多聚物的产生等过程.这里假设有效蛋白对mRNA的产生的负调控是快速的过程,没有延滞.这个调控是通过一个非线性函数,通常是Hill函数来描述的.45\n系统生物学数学基础图3.31:Two-parameterbifurcationdiagramforthemodel(adoptfrom[35]).{fig:4:cirbif}图3.32:Periodoftheendogenousrhythmsofwild-typeandmutantfiles.[35]{fig:4:cirmutant46\n系统生物学数学基础图3.33:(A).Schematicrepresentationofthebiologicalelementsoftheproteinsynthesiscascade,as-sumedtobeelementarytothecircadianrhythmgenerator.Theseincludetheautoinhibitionoftheproteinattranslationalortranscriptionallevelandposttranslationalprocessingsuchasphosphorylation,dimerization,andtransport.Protein∗denotestheeffectiveprotein,be-inginthemolecularstatecapableofinhibitingmRNAproduction,aswellasexpressingthecircadianrhythm.(B).ModelinterpretationofA,emphasizingthedelay(τ)andnonlinearityintheproteinproductioncascade,thenonlinearnegativefeedback,aswellasthemRNAandproteinproductionanddegradation.ThemRNAandproteinproduction(rM,rP)anddegra-dation(qM,qP)rateconstants,respectively,arealsousedastargestforexternalstimulation.(Adoptfrom[32]){fig:4:SCN1}47\n系统生物学数学基础上面的模型可以通过下面时滞微分方程来描述:dMrM=−qMM(3.3.59)dt1+(P/K)ndPm=rPM(t−τ)−qPP.(3.3.60)dt图3.34给出了上述模型的数值模拟结果.可以看到,对于适当的参数,可以产生24小时的节律振荡.−1图3.34:节律振荡模型(3.3.59)-(3.3.60)的模拟结果.这里的参数值为:rM=1.0hr,rP=−1−1−11.0hr,qM=0.21hr,qP=0.21hr,n=2.0,m=3.0,τ=4.0hr,k=1(Replotfrom[32]).{fig:4:scnsim}为了分析节律振荡的产生机制,同样地,我们可以分析平衡点的稳定性.首先,我们作无量纲化处理.令x=M/(rM/qM),y=P/K,˜t=qMt可以得到下面的无量纲化方程′1x=−x(3.3.61)1+yn′my=rx(t−τc)−δy.(3.3.62)d这里′表示,无量纲参数为:dt˜mr=(1/K)(rP/qM)(rM/qM),δ=qP/qM,τc=qMτ.因此,平衡点(x∗,y∗)由下面代数方程的解给出:x∗=1/(1+y∗n),y∗=(r/δ)x∗m.容易看到,该方程有唯一的正平衡点.然后,为得到平衡点的稳定性,可以考虑相应的方程在平衡点处的变分方程x˜′=−x˜−ay˜′(3.3.63){eq:SCN5}y˜=−δy˜+bx˜(t−τc)其中ny∗(n−1)a=,b=mrx∗(m−1).(1+y∗n)2设上面的变分方程有形如x˜=ceλt,y˜=ceλt(c,c6=0)121248\n系统生物学数学基础的解,则可以得到关系1+λac1τcλ=.−beδ+λc2由此,λ满足特征方程f(λ)=(1+λ)(δ+λ)+abe−τcλ=0.(3.3.64){eq:4:cireig}有下面的结论:系统(3.3.63)的零解是稳定的当且仅当特征方程(3.3.64)的所有(复)根都具有负实部.因此,要想得到不稳定的平衡点,只需要特征方程(3.3.63)具有正实部的根.容易看到,当τc=0(没有时滞)时,上面特征方程的根都具有负实部p−(δ+1)±(δ+1)2−4(δ+ab)λ1,2=.2因此,没有时滞,不可能得到不稳定的平衡点.此时,原系统的解不可能得到稳定的周期振荡解.为了得到特征方程有正实部的根的条件,我们来考察具有零实部特征根的条件.为此,令λ=iω,则f(λ)=0等价于−ω2+δ+abcosτω=0,(1+δ)ω−absinτω=0.cc由此可以得到,对给定的δ>0,ab和τc满足关系24222(ab)=ω+(1+δ)ω+δ,其中ω∈(0,π/τc)满足2(ω−δ)tanτcω=(1+δ)ω.上述条件给出了ab−τca平面上的曲线(依赖于δ),该曲线可以把平面分成两个区域,分别对应于平衡点的稳定区域和不稳定区域(图3.35).ab403020Unstable10StableΤc0.20.40.60.81.0图3.35:节律振荡模型(3.3.59)-(3.3.60)的分岔图.三条曲线(从上到下)分别对应于δ=1.5,1.0,0.5.{fig:4:scnstable§3.4胚胎发育生命有机体从单个细胞（受精卵）经过不断分裂、分化、生长和凋亡等形成了各种复合结构，如眼、鼻、耳、口、心和脑等，并由此成长为一个活生生的有机体。这是一个生命的奇迹。基因调控在这一奇迹的实现过程中起着至关重要的作用。了解基因调控在胚胎发育过程中的作用是当前发育生物学领域的研究前沿课题之一。从孟德尔的豌豆实验，摩尔根对果蝇的遗传规律的研究，到Waston和Crick发现脱氧核糖核酸(DNA)的双螺旋结构，以及以后的一系列重大发现，如遗传密码的破译、中心法则的建立、基因分离和克隆、基因的体外重组和基因工程质粒的构建等等，人们已经初步掌握了遗传信息在生命体中代代相传的自然规律。以动物为例，保存在DNA中的遗传信息通过父体和母体的受精作用传递到受精卵，受精卵发育为成熟的个体，再传递给下一代。在这一过程中，同一个体49\n系统生物学数学基础的所有细胞所包含的遗传信息都是一样的，都来源于受精卵。然而，这些具有相同遗传信息的细胞在生物体内却可以分化成上百种诸如骨骼细胞、肌肉细胞、表皮细胞、血细胞和神经细胞等不同类型的细胞。这些分化的各种类型的细胞并不是随机分布的，而是构成复杂的组织和器官，器官又按照一定的方式排列。相同的基因组如何产生不同类型的细胞？这些细胞又是如何形成恰当的排列？这些问题长期困惑发育生物学家。1978年，长期在美国加州理工学院从事果蝇遗传和发育研究的EdwardB.Lewis发表了他关于基因复合体如何控制体节发育的论文[25]。受其影响，欧洲分子生物学实验室两位年轻的发育生物学家ChristianeNusslein-Volhard和EricF.Wieshaus也以果蝇作为模式生物进行研究，试图搞清楚受精卵是如何发育成分节的胚胎的。他们鉴定出15种不同的由于突变引起体节缺陷的基因。当他们的研究结果在1980年发表后[31]，立即受到一批发育生物学家的关注。很快，人们在其他高等生物和人类细胞中发现了同样的或类似的基因，并证明这些基因在发育过程中执行了相似的功能。由于他们的出色工作开创了分子发育生物学研究，EdwardB.Lewis，ChristianeNusslein-Volhard和EricF.Wieshaus分享了1995年诺贝尔生理和医学奖。通过对果蝇体节发育的研究，人们发现了一些能控制胚胎发育的关键性基因，认识到胚胎发育过程中的分化和形态发生是受基因控制的。因此，为了了解胚胎发育过程，就必须了解相应的基因调控网络。该研究领域的中心问题是：控制胚胎发育的基因调控网络是怎么工作的？具体来讲，这个问题包括以下几个方面：1.控制胚胎发育的基因网络是如何组成的？2.它是如何工作的？3.它的可靠性如何得到保证？这一节以影响果蝇翅膀发育的基因调控网络为例子,介绍影响胚胎发育得重要因素:形态生成素(Morphogen)再胚胎得作用.将介绍想光得数学模型,并且对鲁棒性和稳定性等进行分析.果蝇因为其生命周期短、染色体简单等原因，一直是从事分子发育生物学研究的常用模式生物之一。而果蝇翅膀的发育因为和其他器官的发育过程相对独立，也是分子发育生物学家在研究胚胎发育时经常采用的研究对象。然而，即使是果蝇翅膀这一简单器官，其发育过程也是非常复杂的，涉及许多基因之间的相互调控。目前，生物学家通过大量的实验，已经初步了解了与果蝇翅膀发育有关的关键基因，相关信息可以参考果蝇的在线基因信息数据库FlyBase(http://www.flybase.org/)。这些实验事实使我们有可能深入研究相关的基因网络的工作机制。动物的发育模式虽然有很多变化，但发育过程中有许多事件带有普遍性。动物发育可细分为受精(fertilization)、卵裂(cleavage)、原肠形成(gastrulation)、器官发生(organogenesis)、变态(metamorphosis)和成熟(maturity)6个基本阶段，涉及细胞分裂(celldivision)、细胞分化(celldifferentiation)、模式形成(patternformation)、细胞迁移(cellmigration)和细胞凋亡(apoptosis)等一系列事件。果蝇翅膀的发育过程属于胚胎发育中的器官发生阶段。在这一阶段的初始时刻，果蝇胚胎里已经形成了翅膀成虫盘(wingimaginaldisc)，使得翅膀在后阶段的发育是一个相对独立的过程。在随后的发育过程中，细胞通过胚胎的对称轴和体内的形态发生素(morphogen)的浓度分布来确定其在整个胚胎中的位置信息，然后根据其位置确定其最终发育成的细胞形态。胚胎的对称轴包括前后轴(anteroposterioraxis，A-P轴)和背腹轴(dorsoventralaxis,D-V轴)。形态发生素是胚胎中通过形成浓度梯度来影响器官发生阶段细胞的运动和组织形式的物质。果蝇的翅膀成虫盘中的形态发生素包括Decapentaplegic(简称为Dpp)、Hedgehog(简称为Hh)和Wingless(简称为Wg)等。Dpp在成虫盘的A-P轴处合成，然后在细胞外沿垂直于A-P轴方向向两边扩散，形成形态发生素梯度(morphogengradient)。这些细胞外的Dpp和相应的受体(例如Thickveins，简称为Tkv)相结合形成复合体。该复合体被胞吞到细胞内部，作为信号控制目标蛋白Potomotorblind(简称为Omb)和Spalt(简称为Sal)的表达。通过这种方式，沿形态发生素梯度方向的不同细胞将会感受到不同浓度的信号，而根据不同的浓度阈值(concentrationthresholds)表现出不同的基因表达水平。目标蛋白的不同表达水平将决定相应细胞的最终命运。§3.5MorphogengradientMorphogengradientisanimportantconceptindevelopmentalbiology,becauseitdescribesamechanismsbywhichtheemissionofasignalfromonepartofanembryocandeterminethelocation,differentiationandfateofmanysurroundingcells[20].50\n系统生物学数学基础图3.36:Morphogengradient{fig:1}Morphogensaresecretedsignalingmoleculesthatorganizeafieldofsurroundingcellsintopat-terns.Theyformagradientofconcentrationemanatingfromalocalizedsource,anddeterminethearrangementandfateofrespondingcellsaccordingtodifferentconcentrationsofmorphogenper-ceivedbythecells.ThemorphogensassociatewiththedevelopmentofdrosophilawingarelistedatTab.3.3.Theideaofamorphogengradientisintimatelyassociatedwiththeconceptofpositionalinformation[38].Acellisbelievedtoreaditspositioninaconcentrationgradientofanextracellularsignalfactor,andtodetermineitsdevelopmentalfateaccordingly[20].Response†DevelopmentalMorphogenAnti-ResponseSignalsource∗Receptors(concentra-process(range)factor(time,h)tion)DrosophilaimagialAntero-posteriorDppTkvsal(high)‡Brk24-72wingdisc[?]boundary(long)Puntomb(low)neur(high)DrosophilaimagialDorso-ventralWgFzDill(middle)24-72wingdisc[?,?]boundary(long)vg(low)∗Shortrange,20µm;longrange,100µmormore.†Includesrepressionaswellasactivation.‡Onlyafter50hisOmbexpressionfurtherfromthesourcethansal.表3.3:Examplesofmorphogens{tab:1}Themajorfactorsshapingagradientarenotonlydifferentforeachmorphogen,butmaydifferforthesamemorphogenindifferentstagesofdevelopment.Therehasbeenmuchactivityinanalysingthemechanismoftransmissionofamorphogenacrossitsfield.Thereideasprevail:(1)diffusionintheextracellularmatrix;(2)relaybysequentialinternalizationandre-emissionfromcelltocell;and(3)cytoplasmiccontactbythreadsofcytoplasmconnectingdistantcells[20].Thetimingofgradientformationislikelytoberapid.Inlaterdevelopment,asintheDrosophilawingdisc,aDppgradientisnormallyformedslowly,extendingover25celldiametersin3days.Nevertheless,thesameDppgradientcanbereformed,aftertemperatureinterruption,astherateof4cellsin1hour[14,34].InXenopusembryos,however,anactivingradientcanbeformedexperimentallyover100µmin1hour[19],andnaturalgradientsarenormallyformedinXenopusandDrosophilain2hoursorless.Thecellsresponsetothemorphogenthroughdifferentresponsethresholdsofmorphogenconcen-tration.Accordingtothisidea,eachcellwouldhaveonlyabinarychoice:respondtothemorphogenornot.Somecellscanmakemoresophisticateresponsethananon/offswitch[20].Anotherimportantquestionconcerningthecellularbasisofmorphogenperceptionaskswhetheracellneedsitsneighbourstodetermineitspositioninagradient,orwhetheritcanmeasureconcen-trationonitsown.Theexperimentsarguethatcellsinterpretpositioninaconcentrationgradientindependentlyoftheirneighbours[20].Acellmayrespondtomorphogenconcentrationthrougitsreceptorsintwoways.Oneistobearmedwithreceptorshavingdifferentbindingcharacteristics(typeIreceptor);forexample,high-andlow-affinityreceptorsandtheirtransductionpathwayscouldoperateatlowandhighconcentrations51\n系统生物学数学基础ofligand,respectively.Theotheristovarytheoccupancyofonetypeofreceptor(typeIIreceptor),andhenceitssignalingactivity,accodingtoligandconcentration.Wethereforeneedtoknowwhetherdifferentmorphogenresponsesaretransmittedbyoneormorekindsofreceptor.Experimentalresultsshowthatthechoiceofgeneresponsedependsontheabsolutenumberofoccupiedreceptors,entirelyindependentlyofhowmanyunoccupiedreceptorsarepresent[13,20]∗.Herearethreeideasonhowcellsmakedirectresponsetomorphogengradients[20]1.Theavailabilityofligand(morphogen)isthelimitingfactorindeterminingthelevelofresponsetoconcentration.2.Cellsrespondtoligandconcentrationaccordingtotheabsolutenumberofreceptorscooupiedatanytime.3.Acellwithaparticularnumberofoccupiedreceptorswillcontinuetoexpressthesamegeneuntileithertheoccupancyofreceptorsgoesuportheperiodofcompetenceterminates.Anunderstandingofmorphogengradientsrequiresanswerstotwodifferentquestions.Thefirstaskshowadesiredconcentrationgradientisformed.Thesecondquestionaskshowcellsinterpretamorphogenconcentration[20].§3.6模型在果蝇的翅膀盘器官芽中,MorphogenDpp在A-P轴上产生,然后分别向两边扩散,形成从高到低的浓度分布(图??).器官芽的形状如图3.38所示.图3.37:Dppgradient[14].{fig:21}在这里,扩散是主要的机制,使新合成的Dpp转移到胚胎的其他地方.这一机制可以通过扩散方程描述出来.下面给出几种模型(图3.39).模型一:DiffusionandReversibleBinding(Kerszberg&Wolpert[22]).d[L]∂2[L]=D2−konRtot[L](1−[LR])+koff[LR](3.6.65)dtdXd[LR]=konRtot[L](1−[LR])−koff[LR].(3.6.66)dt∗Theoverexpression,byuptotenfold,oftheactivintypeIIreceptorinXenopusembryo.52\n系统生物学数学基础图3.38:WingimaginaldiscofDrosophila[33].{fig:dpp21}这里Rtot表示自由的受体和结合的受体的总和:Rtot=[R]+[LR].在这里和下面的所有模型中,满足下面的边值条件:∂[L]D=v,[L]|X=Xmax≡0.(3.6.67){eq:4:bvc}∂XX=0这个边值条件表示在X=0处又持续的Dpp产生,而在器官芽的边界X=Xmax,多余的Dpp会被其他分子降解掉.模型二:Diffusion,ReversibleBinding,andDegradation(Lander,Nie&Wan[23]).d[L]∂2[L]=D2−konRtot[L](1−[LR])+koff[LR](3.6.68)dTdXd[LR]=konRtot[L](1−[LR])−koff[LR]−kdeg[LR].(3.6.69)dt模型三:Diffusion,Reversible,Binding,ReversibleInternalization,Degradation(Lander,Nie&Wan[23]).∂[L]∂2[L]=D2−kon[L][R]out+koff[LR]out(3.6.70)∂T∂X∂[LR]out=kon[L][R]out−koff[LR]out−kin[LR]out+kout[LR]in(3.6.71)∂T∂[LR]in=kin[LR]out−kout[LR]in−kdeg[LR]in(3.6.72)∂T∂[R]out=−kon[L][R]out+koff[LR]out−kp[R]out+kq[R]in(3.6.73)∂T∂[R]in=ωR−kg[R]in+kp[R]out−kq[R]in(3.6.74)∂T模型四:Diffusion,reversiblebindingwithreceptorandnon-receptor,reversibleinternalization,degradation(Lander,Nie&Wan,2007[24]).53\n系统生物学数学基础(A)(B)(C)(D)图3.39:Dpp扩散模型:(A)Diffusionandreversiblebinding.(B)Diffusion,reversiblebindinganddegradation.(C)Diffusion,reversiblebinding,reversibleinternalization,degradation.(D)Diffusion,reversiblebindingwithreceptorandnon-receptor,reversibleinternalization,degra-dation.{fig:dppmodels}∂[L]∂2[L]=D2−kon[L][R]out+koff[LR]out(3.6.75)∂T∂X−jon[L][N]out+joff[LN]out∂[LR]out=kon[L][R]out−koff[LR]out−kin[LR]out+kout[LR]in(3.6.76)∂T∂[LR]in=kin[LR]out−kout[LR]in−kdeg[LR]in(3.6.77)∂T∂[R]out=−kon[L][R]out+koff[LR]out−kp[R]out+kq[R]in(3.6.78)∂T∂[R]in=ωR−kg[R]in+kp[R]out−kq[R]in(3.6.79)∂T∂[LN]out=jon[L][N]out−joff[LN]out−jin[LN]out+jout[LN]in(3.6.80)∂T∂[LN]in=jin[LN]out−jout[LN]in−jdeg[LN]in(3.6.81)∂T∂[N]out=−jon[L][N]out+joff[LN]out−jp[N]out+jq[N]in(3.6.82)∂T∂[N]in=ωN−jg[N]in+jp[N]out−jq[N]in(3.6.83)∂T对于Dpp,在其开始合成后,很快可以形成定态的浓度分布.而这个定态的浓度分布使确定细胞的发育的关键信息.因此,为了通过数学方法描述上面的定态分布,我们令上面的方程中所有关于时间T的微分为零,从而可以得到一套包含常微分方程和一系列代数方程的系统.通过简化,可以得到含边值54\n系统生物学数学基础条件(3.6.67)的二阶常微分方程.然而,这个方程通常是非线性的,一般不能求解出精确解,我们需要在一定的条件下求出近似解.下面以模型二为例进行简单的分析.在模型二中,令∂[L]∂[R]==0∂TdT可以得到方程∂2[L]0=D2−konRtot[L](1−[LR])+koff[LR]dX0=konRtot[L](1−[LR])−koff[LR]−kdeg[LR].从第二个方程可以求解出[LR]:konRtot[L]koff+kdeg[LR]=,konRtot1+[L]koff+kdeg代入第一个方程,可以得到边值问题konRtot[L]∂2[L]k+koffdegD−kdeg=0∂X2konRtot1+[L](3.6.84)koff+kdeg∂[L]D=v,[L]|X=Xmax=0.∂XX=0令koff+kdeg[L]XK=,u=,x=konRtotKXmaxkX2vXdegmaxmaxψ=,β=KDKD可以得到无量纲化的边值问题′′u′u−ψ=0,u(0)=β,u(1)=0.(3.6.85){eq:4:bvp1}1+u下面,我们来求解方程(3.6.85).两边乘以u′,再对x积分,可以得到ZxZx′′′u′0=uudx−ψudx001+uZxZx′′u=udu−ψdu001+u1′2xx=u|0−ψ(u−log(1+u))|021′22=(u−β)−ψ(u−log(1+u)−(u(0)−log(1+u(0)))2令β2a2=−(u(0)−log(1+u(0)),2ψ可以得到一阶方程pu′=−2ψ(u−log(1+u)+a2),u(1)=0.因此,该方程的解u=u(x)由Zupdup=2ψ(1−x)(3.6.86){eq:4:bvp2}0(u−log(1+u)+a2)55\n系统生物学数学基础给出,其中u(0)由超越方程Zu(0)dupp=2ψ(3.6.87){eq:4:bvp3}0(u−log(1+u)+a2)给出.上面方程的严格表达式不可能得到,只能通过数值方法求解.下面我们在特殊的近似条件下求解上面的方程.为此,令12u−log(1+u)≈u.2这个近似在很多情况下是可以接受的,特别是当u<1时.此时,上面的方程(3.6.86)-(3.6.87)变为Zw(0)dwpZwdwp√=ψ,√=ψ(1−x)0w2+10w2+1u其中w=√.通过上面的方程,可以求解出2apw(x)=sinh(ψ(1−x)).√分布u(x)和w(x)指差一个常数倍2a.例如,通过u(0)2p=w(0)2=sinh2ψ2a2可以求解出√(β2/ψ)sinh2ψu(0)2=.√21+sinhψ由此,可以得到β212a=√.2ψ1+sinh2ψ因此,原来方程的解可以近似为βpu(x)=qsinh(ψ(1−x)).√2ψ(1+sinhψ)√这个解,特别时由sinh(ψ(1−x))所刻划出了的指数衰减,跟实验观察到的结果是符合的.√由上面的解的形式也可以看到,Dpp的浓度衰减特征长度由ψ所刻划出来,这个参数是和[LR]的降解率联系起来的,表征信号分子的衰减速度.当ψ比较很大时,衰减很快,不能形成长距离的Morphogen的分布.当ψ很小时,衰减太慢,不足以产生多种命运的细胞.因此,在生命系统中,ψ的值不能太大也不能太小.根据实验的观察,一般地,Dpp在器官芽的边界处的浓度与产生区域的浓度相比,大概是3%,即[L](Xmax)≈3%[L](0)由此可以得到ψ∼5.前面的近似的条件是u比较小,这个可以通过u(0)来刻划.一般地,Dpp的浓度在器官芽的的边缘处的变化很小,即u′(1)=−a≈0.由此,要想u比较小,只需要β2/ψ很小.在ψ固定的情况下,只需要β小,即Dpp的合成速度很小.这里影响细胞发育的信号是符合体[LR]的浓度,其无量纲化的值由β√qsinh(ψ(1−x))√2u(x)ψ(1+sinhψ)y(x)==1+u(x)β√1+qsinh(ψ(1−x))√2ψ(1+sinhψ)56\n系统生物学数学基础表示.实验发现,果蝇的翅膀的发育对Dpp的基因突变是鲁棒型好的,即使其合成速率明显提到,也不影响翅膀的发育.我们看需要什么样的条件才能满足这样的要求.如果Dpp的合成速率很低(β≪1),则βpy(x)≈u(x)=qsinh(ψ(1−x)).√2ψ(1+sinhψ)可以看到,y(x)和β成线性关系.即生成率增加一倍(这个通过基因突变是很容易实现的),信号的强度也增加一倍,这个结果和试验现象是不符合的.从上面的分析,Dpp的合成速率不应该太低.但是,也不能太高.因为如果Dpp的浓度太高,则在整个器官芽,信号的浓度都很高,无法实现调控细胞的不同命运(图3.40).yHxLyHxL0.8Β=200.4Β=100.60.3Β=2Β=10.20.40.10.2xx0.20.40.60.81.00.20.40.60.81.0(A)(B)图3.40:DppMorphogenGradient(ψ=5).{fig:4:dppgrad}最近的研究结果表明,上面的机制不能够同时保证有生物学意义的信号分布和好的鲁棒性.为了满足这样的条件,需要其他的分子的作用.最近发现,有一种分子,主要是分布在细胞表面,可以和Morpohgen结合,帮助实现Dpp的浓度分布.理论研究表明,这些分子对于实现好的鲁棒性是很关键的.这方面的研究工作还在进行中.§3.7二阶常微分方程边值问题的数学基础在讨论Morphogen的扩散问题时,经常会遇到微分方程的边值问题.这一节介绍相关的基础知识,详细内容可以参考相关的专著.这里的边值问题的一般形式为w′′−f(w,x)=0′′(3.7.88){eq:4s:bvp1}(aw+bw)|x=0=h0,(cw+dw)|x=1=h1.§3.7.1解解解的的的存存存在在在唯唯唯一一一性性性对于边值问题,边界条件的提法非常关键,如果提得不好,解不一定存在,或者存在却并不唯一.上下解的方法是证明解的存在性的常用方法.该方法的主要内容如下.首先,我们可以把上面的方程改写为更加一般的形式.定义H=C2([0,1],R)表示所有p阶可微p的函数,定义线性算子Li:H27→H0,把二阶可微函数映射为连续函数.特别地,定义d2L1u=udx2dL2u=±(a+b)udxx=0dL3u=±(c+d)udxx=157\n系统生物学数学基础则Li满足上面的条件.这里L2和L3中的符号根据可以根据边界条件来确定,需要满足一定的条件.我们在下面详细讨论.下面记线性泛函L:H→H3为L=(L,L,L).另外定义泛函f:H→H32012220为f(u)=(−f(u(x),x),±h0,±h1),∀u∈H2则上面的边值问题变为Lw+f(w)=0.为了求解上面的方程,我们可以通过迭代运算w=(L−1◦f)(w)来求解.但是,L通常不可逆.为此,我们把上面的方程改写为(L−Λ)w+(f(w)+Λw)=0其中Λ=diag{λ1,λ2,λ3}.通过选取适当的Λ,可以使(L−Λ)是可逆的,并且其逆算子满足一定的单调性条件,这样,就可以通过迭代−1w=((L−Λ)◦(f+Λ))w:=T(w)和适当的初始解u0,v0来求解上面的方程.特别地,如果u0T(vi),则可以得到序列u0v0;2.存在常数λi>0(i=1···,k)使得对任意满足u0≥ϕ1≥ϕ2≤v0的函数ϕ1,ϕ2∈H2,fi(ϕ1)−fi(ϕ2)>−λi(ϕ1−ϕ2),3.定义算子Λ:H7→(H)k为Λu=(λu,···,λu).则其逆算子(L−Λ)−1是负算子.221k则方程(3.7.89)至少有一个解u(x)∈H2,并且v0≤u≤u0.证证证明:改写(??)为(L−Λ)u+g(u)=0其中g(u)=f(u)+Λu.58\n系统生物学数学基础则g是单调的,即g(ϕ1)>g(ϕ2)对任意u0≥ϕ1≥ϕ2≥v0.这是因为g(ϕ1)−g(ϕ2)=f(ϕ1)+Λϕ1−f(ϕ2)−Λϕ2>−Λ(ϕ1−ϕ2)+Λ(ϕ1−ϕ2)=0.定义影射T:H0→H2,表示β=Tα满足方程(L−Λ)β+g(α)=0.事实上,我们有Tα=−(L−Λ)−1g(α).下面证明T满足下面的单调性条件.首先,T是单调的,即如果v0≤α1≤α2≤u0,则Tα1≤Tα2.事实上,我们有Tα−Tα=−(L−Λ)−1g(α)+(L−Λ)−1g(α)1211=(L−Λ)−1(g(α)−g(α))21≤0.其次,如果α是上解,则α>Tα.为证明,令β=Tα,则(L−Λ)(β−α)≥0因为(L−Λ)−1是负算子,可以得到β−α≤0.类似地,如果α是下解,则α0.11而β′(1)≤0,因此,一定存在一点x0,β′′(x)=0.但是,这和β′′(x)=12222α(x2)+λβ(x2)>0矛盾.(2).对于问题w′′−f(w,x)=0,w′(0)=h,w′(1)=h(3.7.93){eq:4s:b1}01这里,定义L使得d2Lu=(,−u′(0),u′(1)).dx2则对任意Λ=diag(λ,0,0)59\n系统生物学数学基础算子T=(L−Λ)−1映β=Tα为方程(L−Λ)β=α这里α∈H3.可以验证上面所定义的算子T是负算子.事实上,上面的方程即0β′′−λβ=α(x),−β(0)=α,β′(1)=α.123可以证明,如果α2,α3≥0,并且α1(x)≥0,(∀0τ−tS),状态P中的所有细胞都是被标记的,因此没有被标记的细胞都是状态G0的细胞,其数目是衰减的,衰减率为b=β+δ.由(5.2.3),(5.2.7)和(5.2.8)可以得到2fLγfPγtSb=1−(1−e)(5.2.9){eq:5:12}fN1−e−γtSfL由(5.2.9),(5.2.8),(5.2.7),(5.2.3)可以通过fL,fP,fN=1−fN和tS计算出γ,β,τ,δ.下面来估计tS的范围.为此,首先来估计γ的范围.显然有γ≥0,当γ=0时,可以求解出fLfPβ(γ=0)=δ(γ=0)=,τ(γ=0)=tS.(5.2.10){eq:5:15}tSfNfL由(5.2.3)和δ>0,我们可以得到2eγτ−1>0,或者γτ0和b<0,并且b=µ1(a+δ).因此,可以预见,当δ越大时,(a,b)会超出S的范围,系统的定态解变得不稳定(图5.3(B)).系统的分岔图如图5.4所示.可以看到,当δ增大(干细胞的死亡率增加)时,系统可以出现周期解,并且周期随δ的增加而增加.图5.4:Bifurcationdiagram[11]{fig:5:ss2}§5.4整体模型模型见图5.5方程:dQ−γSτS=−β(Q)Q−(κN(N)+κR(R)+κP(P))Q+2eβ(QτS)QτS(5.4.16)dtdN=−γNN+ANκN(NτN)QτN(5.4.17)dtdR−γRτRS=−γRR+AR{κR(RτRM)QτRM−eκR(RτRM+τRS)QτRM+τRS}(5.4.18)dtdP−γPτPS=−γPP+AP{κP(PτPM)QτPM−eκP(PτPM+τPS)QτPM+τPS}(5.4.19)dt72\n系统生物学数学基础图5.5:造血系统的模型{fig:5:full}whereθsθn21β(Q)=k0θ2+QsκN(N)=f0θn+Nn21κ¯pκ¯rκP(P)=κR(R)=1+KpPm1+KrRr73\n第六章细胞调亡74\n第七章蛋白质折叠与随机动力学75\n参考文献[1]Abkowitz,J.L.,Holly,R.D.,Hammond,W.P.1988.Cyclichematopoiesisindogs:studiesoferythroidburstformingcellsconfirmanearlystemcelldefect,Exp.Hematol.16:941-945.[2]AbkowitzJ,CatlinS,GuttorpP.1996.Evidentthathaematopoiesismaybeastochasticprocessinvivo,NatureMed.2:190-197.[3]AbkowitzJ,GolinelliD,HarrisonD,GuttorpP.2000.Theinvivokineticsofmurinehemopoieticstemcells,Blood96:3399-3405.[4]Alon,U.2006.Anintroductiontosystemsbiology,designprinciplesofbiologicalcircuits.Chap-man&Hall/CRC.London.[5]Atkinson,M.R.,Savageau,M.A.,Myers,J.T.,Ninfa,A.J.2003.DevelopmentofgeneticcircuitryexhibitingtoggleswitchoroscillatorybehaviorinEscherichiacoli.Cell113:597-607.[6]Bernard,S.,B`elair,J.,Mackey,M.C.2003.Oscillationsincyclicalneutropenia:newevidencebasedonmathematicalmodeling,J.Theo.Bio.223:283-298.[7]BeutlerE,LichtmanMA,CollerBS,KippsTJ.1995.Williamshematology,McGraw-Hill,NewYork,1995.[8]Boggs,D.,Boggs,S.,Saxe,D.,Gress,L.,Canfield,D.1982.Hematopoieticstemcellswithhighproliferativepotential:assayoftheirconcentrationinmarrowbythefrequencyanddurationofcureofW/Wvmice,J.Clin.Invest.70:242-253.[9]BradfordG,WilliamsB,RossiR,BertoncelloI.1997.Quiescence,cycling,andturnoverintheprimitivehaematopoieticstemcellcompartment,Exper.Hematol.25:445-453.[10]CheshierS,MorrisonS,LiaoX.WeissmanI.1999.Invivoproliferationandcellcyclekineticsoflongtermselfrenewinghaematopoieticstemcells,Proc.Natl.Acad.Sci.USA96:3120-3125.[11]Colijn,C.,Mackey,M.C.2005.Bifurcationandbistabilityinamodelofhematopoieticregulation.preprintsubmitted.[12]DanceyJT,DeubelbeissKA,HarkerLA,FinchCA,1976.Neutrophilkineticsinman,J.Clin.Invest58:705-715.[13]Dyson,S.,Gurdon,J.B.1998.Theinterpretationofpositioninamorphogengradientasrevealedbyoccupancyofactivinreceptors.Cell93:557-568.[14]Entchev,E.V.,Schwabedissen,A.,Gonzalez-Gaitan,M.2000.GradientformationoftheTGF-βhomologDpp.Cell103:981-991.[15]Fung,E.,Wong,W.W.,Suen,J.K.,Bulter,T.,Lee,S.,Liao,J.C.2005.Asyntheticgene-metabolicoscillator.Nature435:118-122.[16]Gillespie,D.T.1977.Exactstochasticsimulationofcoupledchemicalreactions.J.Phys.Chem.81:2340-2361.[17]Gillespie,D.T.2000.ThechemicalLangevinequation.J.Chem.Phys.113:297–306.76\n系统生物学数学基础[18]Gillespie,D.T.2002.ThechemicalLangevinandFokker-Planckequationsforthereversibleisomerizationreaction.J.Phys.Chem.A106:5063-5071.[19]Gurdon,J.B.,Harger,P.,Mitchell,A.,Lemaire,P.1997.Activinhasadirectlongrangesignalingactivityandcanformaconcentrationgradientbydiffusion.Curr.Biol.7:671-681.[20]Gurdon,J.B.,Bourillot,P.Y.2001.Morphogengradientinterpretation.Nature413:797-803.[21]Hasty,J.,Pradines,J.,Dolnik,M.,Collins,J.J.2000.Noise-basedswitchesandamplifiersforgeneexpression.Proc.Natl.Acad.Sci.USA97:2075-2080.[22]Kerszberg,M.,Wolpert,L.1998.Mechanismsforpositionalsignallingbymorphogentransport:atheoreticalstudy.J.Theor.Biol.191:103-114.[23]Lander,A.D.,Nie,Q.,Wan,F.Y.M.2002.Domorphogengradientsarisebydiffusion?Dev.Cell,2:785–796.[24]Lander,A.D.,Nie,Q.,Wan,F.Y.M.2007.Membrane-AssociatedNon-ReceptorsandMor-phogenGradients.Bull.Math.Biol.69:33-54.[25]Lewis,E.B.1978.AgenecomplexcontrollingsegmentationinDrosophila.Nature,276:565-570.[26]Lipshtat,A.,Loinger,A.,Balaban,N.Q.,Biham,O.2006.Genetictoggleswitchwithoutcoop-erativebinding.Phy.Rev.Lett.96,188101.[27]Kloeden,P.E.,andE.Platen.1992.NumericalSolutionofStochasticDifferentialEquations.Springer-Verlag,NewYork.[28]Huang,G.,WangH.,ChouS.,Nie,X.,Chen,J.,Liu,H.P.2006.BistableexpressionofWOR1,amasterregulatorofwhite-opaqueswitchinginCandidaalbicans.Proc.Natl.Acad.Sci.USA,103,12813-12818.[29]Mackey,M.C.2000.Cellkineticstatusofhaematopoieticstemcells,CellProlif.34:71-83.[30]Micklem,H.,Lennon,J.,Ansell,J.,Gray,R.A.1987.Numbersanddispersionofrepopulatinghematopoieticcellclonesinradiationchimerasasfunctionsofinjectedcelldose,Exp.Hematol15(3):251-257.[31]Nusslein-Volhard,C.,Wieshaus,E.1980.MutationsaffectingsegmentnumberandpolarityinDrosophila,Nature,287:795–801.[32]oldeScheperT.,Klinkenberg,D.,Pennartz,C.,vanPelt,J.1999.Amathematicalmodelfortheintracellularcircadianrhythmgenerator.J.Neurosci.,19,40-47.[33]Tabata,T.2001.Geneticsofmorphogengradients.Nat.Rev.Genet.2:620-630.[34]Teleman,A.A.,Cohen,S.M.2000.DppgradientformationintheDrosophilawingimagianldisc.Cell103:971-980.[35]Tyson,J.J.,Hong,C.I.,Thron,C.D.,Novak,B.1999.AsimplemodelofcircadianrhythmsbasedondimerizationandproteolysisofPERandTIM.Biophys.J.,77,2411-2417.[36]vanKampen,N.G.1992.Stochasticprocessinphysicsandchemistry.North-Holland,Amster-dam,1992.[37]Vilar,J.M.G.,Kueh,H.Y.,Barkai,N.,Leibler,S.2002.Mechanismsofnoise-resistanceingeneticoscillators.Proc.Natl.Acad.Sci.USA,99,5988-5992.[38]Wolpert,L.1989.Positionalinformationrevisited.Development107(Suppl.):3-12.[39]Zordan,R.E.,Galgoczy,D.J.,Johnson,A.D.2006.Epigeneticpropertiesofwhite-opaqueswitchinginCandidaalbicansarebasedonaself-sustainingtranscriptionalfeedbackloop.Proc.Natl.Acad.Sci.USA,103,12807-12812.[40]http://nobelprize.org/nobel_prizes/medicine/laureates/2001/press.html.77\n插图2.1中心法则...........................................132.2Intrinsicandextrinsicnoiseingeneexpression(Elowitzet.al.2002)..........152.3Amodeloftheexpressionofasinglegene.Eachsteprepresentsseveralbiochemicalreactions,whichareassociated2.4模拟结果(Gillespie算法):[mRNA]vs.Time........................162.5模拟结果(求解Langevin方程):[mRNA]vs.Time.....................163.1基因表达的正反馈调控....................................183.2分岔图(λ=4,a=4,n=4).................................203.3DynamicalproperteisofλrepressorcI............................203.4Bifurcationplotsforthevariablexandconcentrationofλrepressor(α=50,σ1=1,σ2=5)223.5模拟结果:α=50,γ=15,σ1=1,σ2=5.(A):确定性系统,γ=15(t<20)和γ=10(t>20).(B):随机系统(3.1.203.6Mutualrepressioncircuit...................................243.7TheprobabilitiesP(NA,NB)fortheswitch,underconditionof(A)weakrepression(k=0.005)wherethereisone3.8ThepopulationofunboundAandBproteinsvstime,obtainfromGillespiesimulationsoftheswitchwithparameters3.9Modulesandconnectivityforthegeneticclock.ThetopconstructcontainsthegInAp2promoterfusedtogInG.3.10Atkinsonoscillator计算模拟结果.这里的无量纲参数为:β1=β3=30.0,β4=1.0,λ1=λ3=2.0,α1=α2=20.03.11不同的参数值所对应的动力学行为.(A):n2=4,α2=10.5;(B):α2=20,n2=5;(C)λ3=0.5,α2=100,n1=4,3.12Conceptualdiagramoftheoscillatorydynamics,highlightingthetwometabolitepools(M1andM2)andtheircon3.13BiologicalrealizationoftheconceptualdesigninFig.3.12.Theyellowboxeshighlightthetwometabolicpools,3.14Computationalcharacterizationofthemetabolator.Themetabolatorispronetooscillateatincreasingglycolyt3.15Biochemicalnetworkofthecircadianoscillatormodel.(replotfrom[37])........313.16Oscillationsinrepressorandactivatorproteinnumbersobtainedfromnumericalsimulationsofthedeterministi3.17计算模拟结果:(A)完整系统,(B)简化系统.........................333.18τ与K的依赖关系.......................................343.19简化方程的向量场分析图示([37])...............................353.20计算模拟结果:(A)确定性系统,(B)随机模型.这里δR=0.005,其他参数和前面一样..353.21Theschematicdiagramofasyntheticgeneregulatorynetwork..............363.22Whennostimulusispresented,thecellsoscillatewiththesamefrequency(a)withdifferentphases(b).Without3.23Thesynchronizationinthepresentofsinusoidalperiodicstimulus.Whenthestimuluswithresonancefrequency3.24Synchronizationinthepresentofwhitenoisestimulus.(a)showsthetimeevlutionoforderparameterofthesystem3.25Theorderparameterfor1000cellswithperiodτvariesfrom5to15,anddifferentvaluesofσ(redcirclesforσ3.26Therelationbetweentheorderparameterandthevariabilityofthecells(measuredby∆β.Theparameterused3.27AsimplemolecularmechanismforthecircadianclockinDrosophila,adoptedfrom[35].PERandTIMproteins3.28数值模拟结果.这里的参数为:v=1,k=0.1,v=0.5,k′=10,k=0.03,k=0.1,P=0.1,J=0.05,rmmpp1p2p3critp3.29(A)g′(P∗)与参数K的依赖关系.(B)相平面图......................45eq3.30(A)′g′(P∗)与参数k的依赖关系.(B)解与参数k的关系.这里k分别取值10,20,30,其它参数和前面一致.p1p1p13.31Two-parameterbifurcationdiagramforthemodel(adoptfrom[35])...........463.32Periodoftheendogenousrhythmsofwild-typeandmutantfiles.[35]..........463.33(A).Schematicrepresentationofthebiologicalelementsoftheproteinsynthesiscascade,assumedtobeelementa−1−1−13.34节律振荡模型(3.3.59)-(3.3.60)的模拟结果.这里的参数值为:rM=1.0hr,rP=1.0hr,qM=0.21hr,qP=03.35节律振荡模型(3.3.59)-(3.3.60)的分岔图.三条曲线(从上到下)分别对应于δ=1.5,1.0,0.5.493.36Morphogengradient.....................................5178\n系统生物学数学基础3.37Dppgradient[14]........................................523.38WingimaginaldiscofDrosophila[33]............................533.39Dpp扩散模型:(A)Diffusionandreversiblebinding.(B)Diffusion,reversiblebindinganddegradation.(C)Diffusion,3.40DppMorphogenGradient(ψ=5).............................575.1CtemCellmodel.......................................675.2细胞分裂过程.........................................685.3参数:稳定区域........................................725.4Bifurcationdiagram[11]...................................725.5造血系统的模型........................................7379