（計(jì)算數(shù)學(xué)專業(yè)論文）基因調(diào)控網(wǎng)絡(luò)的隨機(jī)算法研究.pdf

上海大學(xué)碩士學(xué)位論文 摘要 基因調(diào)控網(wǎng)絡(luò)是由一組基因、蛋白質(zhì)、r n a 、小分子以及它們之間的相互 作用所構(gòu)成的一種生化系統(tǒng)。對(duì)基因調(diào)控網(wǎng)絡(luò)的算法研究在理論和實(shí)踐上有著 重要的意義。 本文首先給出化學(xué)主方程，并介紹幾種精確隨機(jī)模擬算法和快速近似算 法。其次，對(duì)剛性生化系統(tǒng)的混合模擬算法進(jìn)行了介紹，并給出多尺度隨機(jī)模 擬算法，然后在其基礎(chǔ)上發(fā)展了改進(jìn)多尺度隨機(jī)模擬算法。數(shù)值模擬結(jié)果證實(shí) 了改進(jìn)多尺度隨機(jī)模擬算法的優(yōu)點(diǎn)。最后介紹時(shí)滯生化系統(tǒng)的精確模擬算法， 并提出一種模擬時(shí)滯生物化學(xué)反應(yīng)系統(tǒng)的自適應(yīng)f l e a p 算法( d a r l e a p ) 。 該算法能夠較大地提高模擬速度。數(shù)值實(shí)驗(yàn)證實(shí)了d a r l e a p 算法的優(yōu)點(diǎn)。 關(guān)鍵詞：基因調(diào)控網(wǎng)絡(luò)；多尺度算法；改進(jìn)多尺度算法；時(shí)滯生化系統(tǒng)；自 適應(yīng)f l e a p 算法 n 上海大學(xué)碩士學(xué)位論文 a b s t r a c t t h eg e n er e g u l a t o r yn e t w o r k ( g r n ) i sab i o c h e m i c a ls y s t e mw h i c hi sc o m p o s e d o fg e n e s ，p r o t e i n s ，r n a s ，s m a l lm o l e c u l e st o g e t h e rw i t hm e i ri n t e r a c t i v ee f f e c t t h e r e s e a r c ho fa l g o r i t h mf o r ( g r n ) i sv e r ym e a n i n g f u li nt h e o r ya n dp r a c t i c e f i r s t , t h i sp a p e rg i v e sc h e m i c a lm a s t e re q u a t i o nw h i c hd e p i c t st h es t o c h a s t i c s y s t e ma n dp r e s e n t s s o m ee x a c ts t o c h a s t i c s i m u l a t i n g a l g o r i t h m s a n df a s t a p p r o x i m a t ea l g o r i t h m s t h e n , w ei n t r o d u c et h eh y b r i da l g o r i t h mf o rs t i f f n e s s b i o c h e m i c a ls y s t e ma n dd e v e l o pm u l t i s c a l es t o c h a s t i cs i m u l a t i o na l g o r i t h m ( m s s a ) i m p r o v e dm u l t i s c a l e s t o c h a s t i cs i m u l a t i o na l g o r i t h m ( i m s sr ) i sd e v e l o p e d n u m e r i c a le x p e r i m e n t sd e m o n s t r a t et h ea d v a n t a g eo fi m s sr a tl a s t ，t h ee x a c t s t o c h a s t i ca l g o r i t h mf o rs i m u l a t i n gb i o c h e m i c a ls y s t e m sw i t hd e l a yi si n t r o d u c e da n d a l la d a p t i v er _ l e a pa l g o r i t h mf o rs i m u l a t i n gc h e m i c a l l yr e a c t i n gs y s t e mw i t hd e l a y s i sp r o p o s e d t h i sa l g o r i t h ml a r g e l ya c c e l e r a t e st h es p e e do fs i m u l a t i n g n u m e r i c a l e x p e r i m e n t sd e m o n s t r a t et h ea d v a n t a g eo fd ar - l e a pa l g o r i t h m k e yw o r d s ：g r n ；m s s a ；i m s s f ：b i o c h e m i c a ls y s t e m w i t h d e l a y s ；d a r - l e a p i i i 上海大學(xué)碩士學(xué)位論文 原創(chuàng)性聲明 本人聲明：所呈交的論文是本人在導(dǎo)師指導(dǎo)下進(jìn)行的研究工作。 除了文中特別加以標(biāo)注和致謝的地方外，論文中不包含其他人已發(fā) 表或撰寫過(guò)的研究成果。參與同一工作的其他同志對(duì)本研究所做的 任何貢獻(xiàn)均已在論文中作了明確的說(shuō)明并表示了謝意。 簽名：j 亟垂蘧日期：衛(wèi)邋 本論文使用授權(quán)說(shuō)明 本人完全了解上海大學(xué)有關(guān)保留、使用學(xué)位論文的規(guī)定，即： 學(xué)校有權(quán)保留論文及送交論文復(fù)印件，允許論文被查閱和借閱；學(xué) ?？梢怨颊撐牡娜炕虿糠謨?nèi)容。 ( 保密的論文在解密后應(yīng)遵守此規(guī)定) 簽名：避導(dǎo)師簽名弓娩日期：( 壟塑蘭竺 上海大學(xué)碩上學(xué)位論文 1 1 系統(tǒng)生物學(xué) 第一章緒論 系統(tǒng)生物學(xué)n 矗3 4 1 是生物學(xué)的新領(lǐng)域，其目的在于從系統(tǒng)水平上理解生物 系統(tǒng)，力求闡述作為一個(gè)系統(tǒng)的生物系統(tǒng)。系統(tǒng)生物學(xué)是在細(xì)胞、組織、器官 和生物體整體水平上研究結(jié)構(gòu)和功能各異分子之間的相互作用，并通過(guò)計(jì)算機(jī) 模擬來(lái)定量地描述和預(yù)測(cè)生物功能、表型和行為。系統(tǒng)生物學(xué)將在基因組序列 的基礎(chǔ)上完成由生命密碼到生命過(guò)程的研究。這是一個(gè)逐步整合的過(guò)程：由對(duì) 生物體內(nèi)各種分子的鑒別及其相互作用的認(rèn)識(shí)到對(duì)途徑、網(wǎng)絡(luò)、模塊的研究并 最終完成整個(gè)生命活動(dòng)的路線圖。系統(tǒng)生物學(xué)的研究大致有以下四個(gè)方面乜3 1 ： 1 ) 系統(tǒng)結(jié)構(gòu)( s y s t e ms t r u c t u r e s ) 的識(shí)別和研究這包括基因相互作用、生 化代謝途徑以及這些相互作用以何種機(jī)制調(diào)節(jié)生物系統(tǒng)的研究； 2 ) 系統(tǒng)的動(dòng)態(tài)特征( s y s t e md y n a m i c s ) 研究生物系統(tǒng)在不同條件下的行 為，我們可以對(duì)系統(tǒng)進(jìn)行代謝分析( m e t a b o l i ca n a l y s i s ) 、敏感性分析( s e n s i t i v i t y a n a l y s i s ) 、動(dòng)態(tài)分析( d y n a m i ca n a l y s i s ) 如相圖( p h a s ep o r t r a i t ) 和分叉分析 ( b i f u r c a t i o na n a l y s i s ) 等等，從而發(fā)現(xiàn)生物系統(tǒng)的特定行為的機(jī)理，分叉分析又 可以研究系統(tǒng)在多維空間中行為的時(shí)域變化( t i m e - v a r y i n gc h a n g e s ) ，其中每一 維代表系統(tǒng)中一個(gè)生化成分的特定濃度； 3 ) 系統(tǒng)的控制方法( t h ec o n t r o lm e t h o d ) 系統(tǒng)控制一個(gè)細(xì)胞的狀態(tài)能夠 使得細(xì)胞功能受損最小，從而為疾病治療提供潛在的藥物靶標(biāo)； 4 ) 系統(tǒng)的設(shè)計(jì)方法( t h ed e s i g nm e t h o d ) 將原有的生物系統(tǒng)修改和重建成 具有所需功能的系統(tǒng)，其方法和策略可以從確定的設(shè)計(jì)原則、仿真模擬來(lái)發(fā)展， 這樣可以避免盲目的反復(fù)試驗(yàn)( b l i n dt r a i l a n d - e r r o r ) 。 系統(tǒng)生物學(xué)不僅是生命科學(xué)理論的重大發(fā)展，而且具有極其廣闊的應(yīng)用背 景。系統(tǒng)生物學(xué)在人類疾病診斷與治療、藥物的發(fā)現(xiàn)與開(kāi)發(fā)等方面將具有巨大 應(yīng)用潛力。系統(tǒng)生物學(xué)還可應(yīng)用于傳統(tǒng)的生物技術(shù)產(chǎn)業(yè)，開(kāi)拓能源生物技術(shù)、 上海人學(xué)碩士學(xué)位論文 材料生物技術(shù)和環(huán)境生物技術(shù)等新領(lǐng)域。系統(tǒng)生物學(xué)的研究將不僅對(duì)生命科學(xué) 和生物技術(shù)產(chǎn)生重大影響，而且對(duì)整個(gè)國(guó)民經(jīng)濟(jì)、社會(huì)和人類本身產(chǎn)生重大影 響。 1 2 基因調(diào)控網(wǎng)絡(luò)的動(dòng)力學(xué) 系統(tǒng)生物學(xué)的研究?jī)?nèi)容之一就是研究基因調(diào)控網(wǎng)絡(luò)( g e n er e g u l a t o r y n e t w o r k ) 的動(dòng)力學(xué)性質(zhì)?；蛘{(diào)控網(wǎng)絡(luò)是由一組或多組基因、蛋白質(zhì)、小分子 以及它們之間的相互調(diào)控作用所構(gòu)成的一種復(fù)雜的生化系統(tǒng)。 最早的關(guān)于基因調(diào)控網(wǎng)絡(luò)的模型可以追溯到分子生物學(xué)研究的早期。j a c o b 和m o n o d 于1 9 6 3 年發(fā)表的“乳糖操縱子的調(diào)控 開(kāi)了基因調(diào)控網(wǎng)絡(luò)研究的先 聲。由于細(xì)胞中的眾多組件可以通過(guò)復(fù)雜的途徑相互作用，所以大多數(shù)基因調(diào) 控網(wǎng)絡(luò)龐大而且復(fù)雜，如圖1 1 中所示的一個(gè)基因調(diào)控網(wǎng)絡(luò)圖晦1 。因此對(duì)于如 此龐大并且復(fù)雜的調(diào)控過(guò)程的動(dòng)力學(xué)性質(zhì)，僅僅靠觀察和描述的方法是很難理 解的，因此需要引入數(shù)學(xué)工具通過(guò)對(duì)其建立模型來(lái)進(jìn)行分析和模擬。隨著計(jì)算 機(jī)的普及，利用計(jì)算機(jī)進(jìn)行分析和模擬也日益普遍。 牛g 骶2 圖1 1 基因調(diào)控網(wǎng)絡(luò)圖 生命的基本單元一細(xì)胞是由許許多多不同種類的生物大分子所組成，其中 2 上海大學(xué)碩士學(xué)位論文 每一種生物大分子，如一種m r n a 或一種蛋白質(zhì)，在細(xì)胞內(nèi)都存在著若干個(gè)相同 的分子。例如，參與細(xì)胞構(gòu)造的蛋白質(zhì)如肌動(dòng)蛋白( a c t i n s ) 的分子數(shù)目通常 都比較大，可能達(dá)到數(shù)百萬(wàn)個(gè)，而負(fù)責(zé)細(xì)胞活動(dòng)調(diào)控的如轉(zhuǎn)錄因子則分子數(shù)目 一般較少，可能只有幾十個(gè)。但不管是多還是少，在細(xì)胞內(nèi)不存在只有單個(gè)分 子的生物大分子種類。也就是說(shuō)，過(guò)去人們從定性角度談?wù)摰哪衬郴蚧蚰衬?蛋白質(zhì)，實(shí)際上都是某某基因或某某蛋白質(zhì)的若干個(gè)分子的集合體。 這種生物分子集合體的行為與無(wú)機(jī)世界里的分子集合體的行為有著巨大的 不同，其動(dòng)力學(xué)性質(zhì)也有著明顯的區(qū)別。例如，大氣中的各種氣體分子的運(yùn)動(dòng)， 可以用玻爾茲曼方程等統(tǒng)計(jì)物理學(xué)的理論進(jìn)行描述，而生物分子集合體則不會(huì) 遵守同樣的方程。我們知道，新陳代謝是生命存在的必要形式，所有的生物分 子都有著一定的壽命，在生命的每時(shí)每刻，總是有一些分子在“死亡( 降解) ， 有一些分子在“誕生一( 合成) 。因此，在細(xì)胞內(nèi)的每一種分子集合體，都有若 干個(gè)分子不停地被降解，若干個(gè)分子不停地被合成。除了分子數(shù)量的波動(dòng)外， 這種集合體內(nèi)每個(gè)分子的活性也不會(huì)是均一的，有的高，有的低。這種在數(shù)量 和活性上均具有隨機(jī)波動(dòng)的分子，在參與細(xì)胞的活動(dòng)時(shí)會(huì)引入一種擾動(dòng)，稱為 外在噪聲( e x t r i n s i cn o i s e ) 。另一方面，即使在分子的數(shù)量和活性完全一致 的細(xì)胞群體中，反應(yīng)的速率和進(jìn)行的時(shí)間在不同的細(xì)胞中仍然存在著隨機(jī)差異。 這種隨機(jī)性質(zhì)被稱為內(nèi)在噪聲( i n t r i n s i cn o i s e ) ?，F(xiàn)有的研究表明訓(xùn)，隨機(jī) 噪聲廣泛地存在于細(xì)胞的d n a 復(fù)制、基因轉(zhuǎn)錄和蛋白質(zhì)合成。 噪聲的存在對(duì)細(xì)胞的正?；顒?dòng)有重要的影響，它有不利的一面。比如會(huì)影 響細(xì)胞的正常生理過(guò)程，對(duì)基因表達(dá)造成影響，產(chǎn)生突變，引起疾病。但它也 有有利的一面。例如，生物體內(nèi)許多活動(dòng)都具有節(jié)律性，也就是通常所說(shuō)的“生 物鐘 研究人員發(fā)現(xiàn)，噪聲在生物鐘的調(diào)控方面扮演著重要的角色。此外，在 真核生物中，噪聲導(dǎo)致基因表達(dá)的不一致可以作為一種手段，用來(lái)參與細(xì)胞分 化的控制。當(dāng)前，對(duì)生命過(guò)程中的隨機(jī)性的認(rèn)識(shí)及其對(duì)生命活動(dòng)影響的研究正 在成為一個(gè)新的科學(xué)前沿口 引州1 1 幻。更重要的是，由于在細(xì)胞內(nèi)的生命活動(dòng)過(guò) 程中存在著各種隨機(jī)性，使得細(xì)胞的特性和行為表現(xiàn)出種復(fù)雜性和多樣性， 這使得生命體在應(yīng)付各種環(huán)境的變化上具有更強(qiáng)的生存機(jī)會(huì)。這一點(diǎn)不只是體 上海大學(xué)碩上學(xué)位論文 現(xiàn)在整個(gè)物種群體應(yīng)對(duì)環(huán)境的變化上，對(duì)于單個(gè)細(xì)胞也是如此。 當(dāng)然，在細(xì)胞內(nèi)，不同的噪聲源所引起的噪聲對(duì)細(xì)胞的影響是不一樣的。 同樣，細(xì)胞內(nèi)也存在著不同的機(jī)制可以對(duì)噪聲進(jìn)行抑制、放大甚至加以利用。 由于在生命過(guò)程中噪聲無(wú)處不在，它不僅是生命過(guò)程的一部分，也為我們探索 一些生命活動(dòng)提供了一個(gè)途徑。這些噪聲不僅僅與參與生化反應(yīng)的各種分子數(shù) 目有關(guān)，也與濃度、溫度等因素有關(guān)。從通訊的角度來(lái)看，這些噪聲也是一種 通訊手段，可以傳達(dá)各種信息，不但使若干個(gè)反應(yīng)之間的聯(lián)系變得更為微妙， 也在細(xì)胞間的通訊上起著重要的作用。因此，研究生命系統(tǒng)中的隨機(jī)行為，對(duì) 于更好地了解生命過(guò)程以及在計(jì)算機(jī)中重構(gòu)細(xì)胞和生命方面有著重要的意義。 另外由于基因調(diào)控網(wǎng)絡(luò)中生物化學(xué)反應(yīng)的時(shí)間多尺度性，如d n a 與蛋白質(zhì) 結(jié)合以及蛋白質(zhì)聚合等快速反應(yīng)，轉(zhuǎn)錄翻譯和降解等慢速反應(yīng)，以及時(shí)滯特性 ( t i m ed e l a y ) 等使得我們研究基因調(diào)控網(wǎng)絡(luò)中的動(dòng)力學(xué)性質(zhì)變得尤為困難。 已有一些這方面的進(jìn)展 1 4 - 1 9 o 但這些理論都基于這樣一個(gè)假設(shè)前提：所有的生 物化學(xué)反應(yīng)過(guò)程均為服從泊松統(tǒng)計(jì)的馬爾科夫過(guò)程，即不考慮現(xiàn)實(shí)中的時(shí)滯特 性。而在實(shí)際的基因調(diào)控網(wǎng)絡(luò)中，轉(zhuǎn)錄和翻譯等反應(yīng)過(guò)程不僅相對(duì)較慢，而且 常常由多個(gè)基元反應(yīng)組成，這種過(guò)程服從帶有時(shí)滯的高斯統(tǒng)計(jì)，研究這種時(shí)滯 的調(diào)控網(wǎng)絡(luò)的解析方法和數(shù)值模擬方法也就尤為重要乜1 刪。 1 3 本文的工作及安排 本文重點(diǎn)研究基因調(diào)控網(wǎng)絡(luò)的隨機(jī)模擬算法，給出其數(shù)學(xué)模型以及常見(jiàn)的 數(shù)值模擬算法晦刎。第二章介紹連續(xù)確定性模型和離散隨機(jī)性模型，給出幾種精 確的蒙特卡羅算法和快速近似算法。針對(duì)生化反應(yīng)系統(tǒng)的多尺度現(xiàn)象，第三章 介紹幾種有效克服多尺度的算法，給出c a o 等人提出的多尺度隨機(jī)算法，并且 在其基礎(chǔ)上發(fā)展了改進(jìn)多尺度隨機(jī)算法。該算法比多尺度算法能更有效地模擬 剛性系統(tǒng)。第四章介紹時(shí)滯基因調(diào)控網(wǎng)絡(luò)的精確時(shí)滯隨機(jī)模擬算法。針對(duì)模擬 復(fù)雜的時(shí)滯系統(tǒng)的低效性，提出了模擬時(shí)滯生化系統(tǒng)的自適應(yīng)f l e a p 算法 ( d a r l e a p ) 。該算法能夠自動(dòng)調(diào)節(jié)f ，使得在時(shí)間區(qū)間【f ，t + r ) 內(nèi)發(fā)生多次反 4 上海人學(xué)碩士學(xué)位論文 應(yīng)事件并且能夠精確地滿足l e a p 條件。數(shù)值試驗(yàn)證實(shí)了d a r l e a p 算法能夠有 效地模擬復(fù)雜的時(shí)滯生化網(wǎng)絡(luò)。 論文安排如下： 本章闡述了系統(tǒng)生物學(xué)的背景，簡(jiǎn)要介紹了基因調(diào)控網(wǎng)絡(luò)的相關(guān)背景，闡 述了本文的研究目的和意義，并介紹了目前國(guó)內(nèi)外在這方面的研究進(jìn)展。 第二章分別介紹了基因調(diào)控網(wǎng)絡(luò)的兩種動(dòng)力學(xué)模型：連續(xù)確定性模型和離 散隨機(jī)模型以及常用的算法。 第三章介紹剛性生化網(wǎng)絡(luò)的特性和一些常用克服剛性的隨機(jī)算法。給出多 尺度算法，并且在其基礎(chǔ)上提出了改進(jìn)多尺度算法。 第四章介紹時(shí)滯生化反應(yīng)系統(tǒng)的精確算法并且提出了o a r l e a p 算法。 第五章對(duì)本文進(jìn)行總結(jié)與展望。 上海大學(xué)碩? l 學(xué)位論文 第二章基因調(diào)控網(wǎng)絡(luò)的動(dòng)力學(xué)模型及其模擬算法 從計(jì)算角度出發(fā)，系統(tǒng)生物學(xué)中基因調(diào)控網(wǎng)絡(luò)的研究大致可以分為以下兩 條路線：一支是基于數(shù)據(jù)挖掘的基因調(diào)控網(wǎng)絡(luò)的重構(gòu)，該領(lǐng)域主要是從大量的 試驗(yàn)數(shù)據(jù)( 如相互作用數(shù)據(jù)) 中提取信息，如從g e n b a n k 、p d b 等數(shù)據(jù)庫(kù)中采用 神經(jīng)網(wǎng)絡(luò)、隱馬氏模型、支持向量機(jī)等方法重構(gòu)基因調(diào)控網(wǎng)絡(luò)，這些方法都是 一種靜態(tài)的研究手段，主要研究基因調(diào)控網(wǎng)絡(luò)的拓?fù)湫再|(zhì)。而另一支研究路線 則是從整體出發(fā)研究基因調(diào)控網(wǎng)絡(luò)( 如反應(yīng)通道、單細(xì)胞和多細(xì)胞通信等) 的 動(dòng)力學(xué)性質(zhì)，這是一種動(dòng)態(tài)的研究方法。因此需要建立設(shè)計(jì)能夠準(zhǔn)確描述細(xì)胞 內(nèi)生物化學(xué)網(wǎng)絡(luò)的動(dòng)力學(xué)模型。生化網(wǎng)絡(luò)中的化學(xué)反應(yīng)呈現(xiàn)遠(yuǎn)離熱力學(xué)平衡的 現(xiàn)象啪瑚1 ，并且分子復(fù)制數(shù)目有時(shí)候非常低，化學(xué)反應(yīng)具有明顯的隨機(jī)性。但 是當(dāng)分子數(shù)目巨大時(shí)，隨機(jī)不明顯可以忽略不計(jì)，可以用連續(xù)的確定的微分方 程來(lái)描述。一般的生物化學(xué)反應(yīng)系統(tǒng)中參加反應(yīng)的分子數(shù)目通常不是很大，而 隨機(jī)性則在反應(yīng)中起著重要的作用，因此可用化學(xué)主方程刻畫其動(dòng)力學(xué)行為。 本文主要討論基因調(diào)控網(wǎng)絡(luò)中的動(dòng)力學(xué)性質(zhì)的數(shù)值模擬。由于對(duì)基因調(diào)控 網(wǎng)絡(luò)的刻畫的不同，其數(shù)值模擬大致可以分為以下2 種不同的方法：( 1 ) 以連 續(xù)的方式處理，變量為反應(yīng)分子的濃度，表示為確定性方程( 組) ，在2 1 小節(jié) 給出；( 2 ) 以離散的方式處理，變量為反應(yīng)分子的個(gè)數(shù)，反應(yīng)過(guò)程表示為隨機(jī) 過(guò)程，這一方法在處理所含分子數(shù)目較低的系統(tǒng)時(shí)就顯得尤為重要，是本文研 究的重點(diǎn)，在剩余的小節(jié)中給出。 2 1 連續(xù)確定性動(dòng)力學(xué)模型 當(dāng)生物反應(yīng)系統(tǒng)的反應(yīng)分子的數(shù)目足夠大，其隨機(jī)性可以忽略時(shí)，其數(shù)目 可以認(rèn)為是與時(shí)間有關(guān)的連續(xù)變量。我們可以用反應(yīng)物的濃度作為變量的微分 方程來(lái)刻畫基因調(diào)控網(wǎng)絡(luò)的動(dòng)力學(xué)性質(zhì) 6 上海大學(xué)碩士學(xué)位論文 掣釩啉啉 c ( m d c 芴2 ( 一t ) = 正( q ( f ) ，c 2 ( f ) ，c ( f ” ( 2 1 1 ) 1 d c n 廠( t ) = 厶( c l ( f ) ，c 2 ( f ) ，c a f ) ) 其中n 為反應(yīng)分子的個(gè)數(shù)，q ( t ) 為反應(yīng)物濃度。求解方程組( 2 i 1 ) 可以得 到調(diào)控網(wǎng)絡(luò)中各反應(yīng)分子的濃度隨時(shí)間的變化情流。關(guān)于微分方程( 組) 的數(shù) 值方法相對(duì)比較成熟，如歐拉方法( e u l e r ) 、龍格一庫(kù)塔方法( r u n g e - k u t t a ) 虛盤 守。 由于生物化學(xué)反應(yīng)的時(shí)間多尺度性，反應(yīng)系統(tǒng)中經(jīng)常存在有些反應(yīng)很慢， 而有些反應(yīng)則很快。在準(zhǔn)平衡態(tài)假設(shè)下，原微分方程組( 2 1 1 ) 簡(jiǎn)化成由微分 方程和代數(shù)方程組成的耦合系統(tǒng)，通過(guò)求解這種耦合系統(tǒng)，我們同樣可以得到 調(diào)控網(wǎng)絡(luò)中各反應(yīng)分子的濃度隨時(shí)陽(yáng)l 的變化情況。 2 2 離散隨機(jī)動(dòng)力學(xué)模型 當(dāng)細(xì)胞中的生物化學(xué)反應(yīng)的隨機(jī)性不能忽略時(shí)，連續(xù)確定性系統(tǒng)不能正確地 預(yù)測(cè)生物化學(xué)反應(yīng)系統(tǒng)的動(dòng)力學(xué)性質(zhì)。因此必須需要能描述其隨機(jī)性的精確數(shù) 學(xué)模型。本小節(jié)首先介紹生物化學(xué)反應(yīng)網(wǎng)絡(luò)，然后給出其動(dòng)力學(xué)模型。 考慮溫度恒定、體積不變和空間均勻的非時(shí)滯生物化學(xué)反應(yīng)系統(tǒng)。該系統(tǒng) 具有種化學(xué)分子s = s l ，s ) ，通過(guò)m 個(gè)化學(xué)反應(yīng)通道r = r l ，一，r m ) 作 用。系統(tǒng)在時(shí)刻t 的狀態(tài)向量為x ( t ) = ( x l ( f ) ，x ( f ”，其中x f ( f ) 為第i 種化 學(xué)分子在時(shí)刻t 的分子數(shù)目，f = l ，n 。收= ( v 小，) 為通道尺i 發(fā)生一次 反應(yīng)的分子改變向量，其中k = l ，m 。這m 個(gè)生物化學(xué)反應(yīng)形成一個(gè)相互作 用網(wǎng)絡(luò)： 足：v l l s + v 1 2 逆+ + h 凡山q l s + q 2 是+ + 峨凡 v 2 1 s + 屹口叉+ + 吃凡山t l s + 吃曼+ + 囈 ：v u l s l + v u 墨+ + & 山屹l 墨+ 屹2 蔓+ + 7 上海人學(xué)碩士學(xué)位論文 其中q ，k - - 1 9 , , , 9 m 為隨機(jī)反應(yīng)速率常數(shù)，它只與分子的物理性質(zhì)和系統(tǒng)的反應(yīng) 溫度以及體積有關(guān)。若記僅為反應(yīng)通道盈中反應(yīng)分子的表達(dá)式，a k ( x ( f ) ) 為通 道尺i 的反應(yīng)傾向函數(shù)，則有表達(dá)式吼( x ( f ) ) - c a ( x ( t ) ) 。此式的含義將在后面 給出。反應(yīng)通道咒的反應(yīng)類型如下表所示： 表2 1 反應(yīng)通道的類型以及吃的表達(dá)形式 其中五為分子最的分子數(shù)目。 給定x ( j ) = 石，在無(wú)窮小時(shí)間區(qū)間p ，s + 出) 內(nèi)，通道凰發(fā)生反應(yīng)的概率為 a k ( x o ) ) 凼+ d ( 凼) 。即 。 p ( 艘1 0 + 凼) 一眠( s ) = 1 f x ( f ) ，f j ) = a k ( x o ) ) 出+ 口( 凼) 并且兩個(gè)及以上的通道發(fā)生反應(yīng)的概率為d ( 凼) 。從而有下式成立 p ( x ，f + 出) = p ( x ) + 口，( x 一吩娥x 一_ ，t ) d t - 巳( x ) 尸( x ，t ) d t + o ( d t ) 當(dāng)d t o 時(shí)，這就變威化學(xué)主方程( c h e m i c a lm a j = s i t e re q u a t i o n ，c m e ) ： 掣= 姜儼一v j ) e ( x 刊一姜咿m f ) c m e 的直接求解非常困難，因?yàn)榉匠痰木S數(shù)等于可能的分子狀態(tài)數(shù)。于是 我們通過(guò)逼近方法來(lái)研究主方程的性質(zhì)。其中?？艘黄绽士吮平?( f o k k e r p l a n c ka p p r o x i m a t i o n ) 是一種常用的方法。該方法是在狀態(tài)x 附近 上海大學(xué)碩士學(xué)位論文 對(duì)方程進(jìn)行泰勒展開(kāi)，然后從二階導(dǎo)數(shù)后截?cái)?，從而得到不同狀態(tài)之間概率流 的一種連續(xù)逼近。另外b r i a n 等人啪1 提出了有限狀態(tài)影射算法( f i n i t es t a t e p r o j e c t i o na l g o r i t h m ，f s p ) 對(duì)狀態(tài)截?cái)鄟?lái)逼近c(diǎn) m e 。對(duì)于c m e 的維數(shù)不是 很大時(shí)，f s p 方法能夠提供精確的解。當(dāng)c m e 的維數(shù)很大時(shí)，f s p 方法能夠提供 保證精度的解。雖然該算法能夠有效地解一些簡(jiǎn)單的生物化學(xué)反應(yīng)系統(tǒng)，但是 對(duì)于復(fù)雜的生物化學(xué)反應(yīng)系統(tǒng)仍然是不可解的。許多蒙特卡羅算法被提出用來(lái) 克服解c m e 的困難。下面兩小節(jié)給出一些重要的算法。 2 3 精確隨機(jī)模擬算法 為了模擬生物化學(xué)反應(yīng)系統(tǒng)的軌跡，g i l l e s p i e 啪瑚3 提出了直接反應(yīng)法 ( d i r e c t i o nr e a c t i o nm e t h o d ，d m ) 和第一反應(yīng)法( f i r s tr e a c t i o nm e t h o d ， f r m ) 。在2 0 0 0 年，g i b s o n 等人脅3 在第一反應(yīng)法的基礎(chǔ)上給出了下一反應(yīng)法( n e x t r e a c t i o nm e t h o d ，n r m ) 。在2 0 0 7 年，d a v i d m l 發(fā)展了改進(jìn)下一反應(yīng)法( m o d i f i e d n e x tr e a c t i o nm e t h o d ，洲r(nóng) m ) 。這些隨機(jī)算法的解等價(jià)于c m e 的解，因而是精 確隨機(jī)模擬算法。我們首先給出4 種算法，然后對(duì)4 種算法做出點(diǎn)評(píng)。 對(duì)于任何精確隨機(jī)模擬算法，在每一個(gè)時(shí)間步長(zhǎng)上，系統(tǒng)只有一個(gè)狀態(tài)。 在給定狀態(tài)x ( t ) = 工下，必須通過(guò)回答下述兩個(gè)問(wèn)題完成反應(yīng)系統(tǒng)的模擬：1 ) 什么時(shí)候發(fā)生下一個(gè)反應(yīng)? 2 ) 下一個(gè)反應(yīng)是哪個(gè)通道? 具體地，若f 與_ ，分別 表示下一反應(yīng)發(fā)生的時(shí)間和發(fā)生反應(yīng)的通道，精確隨機(jī)模擬算法通過(guò)描述下述 事件e 而完成模擬整個(gè)過(guò)程： e ：時(shí)間區(qū)間 t , t + r ) 內(nèi)不發(fā)生任何反應(yīng)，且r ，將在無(wú)窮小區(qū)間 【t + t ，t + r + d r ) 發(fā)生一次反應(yīng)。 p ( f ，j ) d r ：事件e 發(fā)生的概率； g i l l e s p i e 給出了概率密度p ( f ，j ) ： 即棚= 吩e x p 卜軍q ) 億3 - ， 由方程( 2 3 1 ) ，分別對(duì)_ ，求和以及對(duì)f 積分，可以得到下一反應(yīng)尺，發(fā)生的概 9 上海大學(xué)碩士學(xué)位論文 率分布和下一反應(yīng)發(fā)生的時(shí)間f 的概率密度： g ( r ) = a o ( x ) e x p ( 一a o ( 工) r ) ，f 0 ，弓( 肛麗a j ( x ) 其中a o ( 力= 口j ( x ) 。從u ( o ，1 ) 內(nèi)獨(dú)立產(chǎn)生2 個(gè)隨機(jī)數(shù)，吃，根據(jù)上述概率密 度函數(shù)，可以得出 f = 一衄r t ) a o ( 工) ， ( 2 3 2 ) 而為滿足下述不等式的最小正整數(shù) 圭q r 2 a o ( x ) ( 2 3 3 ) i = l d m 算法總結(jié)如下： 步驟1 初始化系統(tǒng)狀態(tài)x ( 0 ) ，設(shè)定反應(yīng)初始時(shí)刻= 0 ，結(jié)束時(shí)間0 ； 步驟2 根據(jù)系統(tǒng)狀態(tài)x ( f ) 計(jì)算每個(gè)反應(yīng)的傾向函數(shù)( = l ，m )a o ； 步驟3 從u ( o ，1 ) 產(chǎn)生隨機(jī)數(shù)吒，匕； 步驟4 根據(jù)式( 2 3 2 ) 和( 2 3 3 ) 計(jì)算下一反應(yīng)通道及其發(fā)生時(shí)間間 隔a t ； 步驟5 更新反應(yīng)時(shí)間和系統(tǒng)狀態(tài)：t = t + a t ，x = x + 匕； 步驟6 如果f 砭( 瓦) ) 。假定系統(tǒng)的狀態(tài)為x ( f ) = 工時(shí)，假設(shè) 知最，互。通道r 下一次發(fā)生的時(shí)間間隔為似：墨號(hào)孚。選擇最小的時(shí)間間隔 吼l 叫 出，- - - m z a , , ，則系統(tǒng)下一反應(yīng)的時(shí)間間隔為越，發(fā)生反應(yīng)的通道為r 。更 新系統(tǒng)： f = f + 馘互= 正+ 口( 石) 她和x = 工+ 和對(duì)于通道r 更新 只= 露+ l n ( 1 r ) ，其中廠u ( o ，1 ) 。重復(fù)上述過(guò)程，得到m n r m 算法。 刖r m 算法如下： 步驟1 初始化系統(tǒng)狀態(tài)，初始反應(yīng)時(shí)刻乞= o ，結(jié)束時(shí)間為f ， 最= 互= o ，七= 1 ，m ； 步驟2 根據(jù)系統(tǒng)狀態(tài)石( f ) 計(jì)算每個(gè)反應(yīng)的傾向函數(shù)q ( 后= 1 ，m ) ； 1 2 上海大學(xué)碩士學(xué)位論文 步驟3 從u ( o ，1 ) 產(chǎn)生m 個(gè)獨(dú)立變量氣； 步驟4 對(duì)每一個(gè)后，設(shè)最：蚍馬； ，七 步驟5 對(duì)每一個(gè)后，設(shè) ：墨二互； “七 步驟6 i 發(fā)a t i = 卿 ，更新系統(tǒng)：t = t + a t ，瓦= 瓦+ 口( x ) a t 。和z = 石+ 和對(duì)于通道局更新只= 只+ i n ( 1 r ) ，其中reu ( o ，1 ) ； 步驟7 更新傾向函數(shù)吼( k = l 9 * * o9 m ) ； 步驟8 如果t f ，返回到步驟2 循環(huán)，否則結(jié)束。 通常情況下n r m 和心r m 要比其它兩種算法要快一些，然而有些情況使用 d m 是最快的，因?yàn)閚 r m 和m n r m 雖然每次迭代只需用一個(gè)隨機(jī)數(shù)，但是它查找 最小時(shí)間，也需要消耗一定的c p u 時(shí)間，具體討論參見(jiàn)。 3 2 。 2 4 快速的隨機(jī)模擬算法 上述算法為精確模擬算法，由于精確地模擬每一個(gè)反應(yīng)事件導(dǎo)致計(jì)算量很 大，因而對(duì)于模擬復(fù)雜的生化系統(tǒng)是低效的。為了加快模擬的速度，g i l l e s p i e 等人提出了f l e a p 算法棚1 。在給定時(shí)間區(qū)間 t ，f + f ) 內(nèi)系統(tǒng)狀態(tài)改變很小 時(shí)，該算法利用泊松變量生成每一個(gè)反應(yīng)通道在該區(qū)間內(nèi)的反應(yīng)次數(shù)。該類算 法能在給定的時(shí)間區(qū)間【f ，f + r ) 上發(fā)生多次反應(yīng)，故對(duì)于分子數(shù)目很大的系統(tǒng)是 比較有效的。對(duì)某類多通道分子反應(yīng)系統(tǒng)，在整個(gè)反應(yīng)過(guò)程中，其中一些反應(yīng) 通道發(fā)生反應(yīng)的頻率與另外一些反應(yīng)通道發(fā)生反應(yīng)的頻率相比，要大幾個(gè)數(shù)量 級(jí)，并且一些分子參加反應(yīng)頻率比另外一些分子參加反應(yīng)的頻率要大得多。系 統(tǒng)的狀態(tài)在短時(shí)間內(nèi)，主要由反應(yīng)頻率高的反應(yīng)通道刻畫，在長(zhǎng)時(shí)間內(nèi)，主要 由反應(yīng)頻率低的反應(yīng)通道來(lái)刻畫。這種現(xiàn)象叫做多時(shí)間尺度現(xiàn)象( 簡(jiǎn)稱多尺度 現(xiàn)象) ，在數(shù)學(xué)上稱做“剛性 。許多算法被提出用來(lái)克服系統(tǒng)的剛性，將在下 一章給出。這一小節(jié)，我們介紹模擬分子數(shù)目多的系統(tǒng)的算法。首先我們介紹 g i l l e s p i e 的f l e a p 算法，該類算法允許反應(yīng)通道在給定時(shí)間區(qū)間【f ，f + f ) 內(nèi) 上海大學(xué)碩士學(xué)位論文 發(fā)生多次反應(yīng)，但必須滿足l e a p 條件：系統(tǒng)的反應(yīng)傾向函數(shù)幾乎不變。精確 的表述如下：給定時(shí)間區(qū)間 f ，f + f ) 和系統(tǒng)狀態(tài)x ( f ) = 工，r s 通道反應(yīng)的次數(shù)由 泊松隨機(jī)變量巧( r ；功= p o i s s i o n ( a ，( x ) f ) 逼近。在時(shí)間區(qū)間i t ，t + r ) 內(nèi)的系統(tǒng)狀 態(tài)改變量為： m x ( t + r ) - x = ( f ；x ) = 髟( r ；x ) v s ( 2 4 1 ) l = l 滿足l e a p 條件為： i a j ( 石+ ( f ；功) 一a ( 石) i 翻o ( x ) ，( = 1 ，m ) ( 2 4 2 ) 其中s 為預(yù)設(shè)誤差參數(shù)( 0 s 1 ) 。g i l l e s p i e 等人h 1 給出了滿足( 2 4 2 ) 式 的最大f 值： 仁嘶 踹，篇 仫4 舢【i 一( x ) i 7 i 巧( x ) i j 其中 ，( 砷= 厶( x ) 口，( 石) k 1 ( 2 4 5 ) m 巧( 曲= 廠2 ( 石) 口( 功 硝力= 善攀乩，m 億4 利用上述選取時(shí)間步長(zhǎng)_ r 的方法，我們可以得到r - l e a p 算法 步驟l 初始化系統(tǒng)狀態(tài)，初始反應(yīng)時(shí)刻= o ，結(jié)束時(shí)間為0 步驟2 在時(shí)刻t 系統(tǒng)狀態(tài)為x ，計(jì)算所有反應(yīng)的傾向函數(shù)以及它們的和 肼 口o ( x ) = 乃； 步驟3 根據(jù)式( 2 4 3 ) 計(jì)算f ； 步驟4 如果步驟2 中得到的時(shí)間步長(zhǎng)f 小于l 和某個(gè)數(shù)( 不妨假設(shè)是1 0 的乘積，則放棄這個(gè)時(shí)間步長(zhǎng)，轉(zhuǎn)而接連執(zhí)行一定個(gè)數(shù)( 不妨設(shè)為 1 0 0 ) 的精確隨機(jī)模擬( s s a ) 過(guò)程，然后再嘗試獲取新的時(shí)間步長(zhǎng)f 。 1 4 上海大學(xué)碩上學(xué)位論文 反之，如果步驟2 中得到的時(shí)間步長(zhǎng)f 大于1 a o 和那個(gè)數(shù)的乘積， 則接受這個(gè)時(shí)間步長(zhǎng)r ，向下進(jìn)行到步驟4 ； 步驟5 對(duì)于每個(gè)反應(yīng)_ ，( j = l ，m ) ，產(chǎn)生一個(gè)均值為a i ( 工) f 的p o i s s o n 隨 機(jī)變量k ，代表反應(yīng)，在時(shí)間步長(zhǎng)f 內(nèi)發(fā)生的次數(shù)； 步驟6 更新f = f + f ，x = x q - 巧_ ； 步驟7 如果t t ，返回到步驟2 循環(huán)，否則結(jié)束。 當(dāng)m ，很大時(shí)，式( 2 4 5 ) 和( 2 4 6 ) 的計(jì)算量為o ( m 2 + 2 m 2 ) ，為 了有效地降低其計(jì)算量，c a o 等人3 7 1 提出了對(duì)于f l e a p 算法的有效的步長(zhǎng)選 擇算法。該算法是基于反應(yīng)分子的相對(duì)變化的。對(duì)于上述的s 有 f ：m i n ， m a x e x , ( t ) g , , 1 ，塑墮娶幽 ( 2 4 7 ) 把k ll 忍( x ( f ) ) l 7 茸( x ( f ) ) i 其中 露( x ( f ) ) = a j ( x ( t ) ) v v ，i l ，= 1 ( 2 4 8 ) 辭( x ( f ) ) = 吩( 彳( f ) 鶘， f l l 為反應(yīng)分子集合，詳細(xì)請(qǐng)參考文獻(xiàn) 3 7 ，( 2 4 8 ) 式的計(jì)算量為2 m 2 。故 當(dāng)很大時(shí)，能夠有效地降低計(jì)算量，提高速度。 上述的算法都是假定在給定時(shí)間區(qū)間【t ，f + f ) 內(nèi)系統(tǒng)狀態(tài)改變很小時(shí)，該算 法利用泊松隨機(jī)數(shù)確定每一個(gè)反應(yīng)通道在該區(qū)間內(nèi)的反應(yīng)次數(shù)。由于泊松數(shù)的 無(wú)界性，有可能產(chǎn)生的反應(yīng)次數(shù)導(dǎo)致某些分子變成負(fù)的，但是負(fù)的分子數(shù)不符 合現(xiàn)實(shí)。為了避免負(fù)分子數(shù)的出現(xiàn)，t i a n 等人口1 和c h a t t e r j e e 等人嘲獨(dú)立地提 出了二項(xiàng)f l e a p 算法。因?yàn)槎?xiàng)分布產(chǎn)生的隨機(jī)數(shù)是有界的，因此二項(xiàng) f l e a p 算法分布可以避免負(fù)分子數(shù)目的產(chǎn)生。c a i 等人翰1 提出的k - l e a p 算 法能夠減少產(chǎn)生負(fù)分子數(shù)目的概率。c a o 等人洶1 發(fā)展了一種避免負(fù)分子數(shù)產(chǎn)生 的算法。該算法將反應(yīng)通道分成容易產(chǎn)生負(fù)分子數(shù)的關(guān)鍵通道和不容易產(chǎn)生負(fù) 分子數(shù)的非關(guān)鍵通道。用標(biāo)準(zhǔn)的g i l l e s p i e 隨機(jī)算法模擬關(guān)鍵通道，用f l e a p 算法模擬非關(guān)鍵通道。這樣大大地降低了產(chǎn)生負(fù)分子數(shù)的概率，但該類算法在 確定了時(shí)間步長(zhǎng)f 后，就不能再進(jìn)行調(diào)整了，屬于先驗(yàn)f l e a p 算法。最近d a v i d 1 5 上海人學(xué)碩士學(xué)位論文 提出了將后驗(yàn)的l e a p 檢查點(diǎn)植入f l e a p 算法1 ，稱之為后驗(yàn)r - l e a p 算法。 該算法利用單位速率泊松過(guò)程可以自動(dòng)地調(diào)節(jié)時(shí)間步長(zhǎng)f ，使得精確地滿足 l e a p 條件，從而產(chǎn)生負(fù)分子數(shù)的概率為o ，從而可以保證精度。下面我們介紹 后驗(yàn)f l e a p 算法。 單位速率泊松過(guò)程k 具有以下兩種性質(zhì)： a ) 對(duì)于r i m 互瓦( f ) ，增量砭( 乃) 一砭( 五) = k ( 互一石) 是與系統(tǒng)狀態(tài)x ( f ) 是獨(dú)立 的參數(shù)為互一互的泊松分布； b ) 給定o s 彳 露= 互。對(duì)于步長(zhǎng)f ，通道r 的內(nèi)部時(shí)間為疋+ 吼( x ( f ) ) r 時(shí)， 發(fā)生的反應(yīng)次數(shù)k 概率分布分兩種情況： 如果存在f o ，一，e 耐一1 ) 使得t t k 五+ ( x ( f ) 弦 硭成立，根據(jù)性質(zhì)b 知 砭( 瓦+ q ( x ( f ) ) f ) = b i n o m i a t ( y , ( t l * + i ) 一k ( 互) ，) + k ( 乃) 其中 如果疋+ 吼( x ( f ) 弦吃 k ( 正+ 口i ( x o ) ) f ) = p o i s s o n ( t t + 口i ( x ( ，) ) f 一互乞) + k ( z 乞) 則在時(shí)間區(qū)間 f ，f + f ) 發(fā)生的反應(yīng)次數(shù)m = k ( 瓦+ ( x ( f 弦) 一q 代入狀態(tài)方程 ( 2 4 1 ) 得到狀態(tài)x ( t + r ) ，檢查是否滿足l e a p 條件。如滿足則接受f ，更新 系統(tǒng)。如不滿足，則將k ( 疋+ 口。( x ( f ) f ) ，互+ 噸( x ( f 弦按升序插入到已存儲(chǔ)的 1 6 蔓 紫 上海大學(xué)碩士學(xué)位論文 信息里。用更小的f 重復(fù)上述過(guò)程，直到找到符合條件的f 。 后驗(yàn)f - l e a p 算法： 步驟l 初始化：初始時(shí)間及初始分子數(shù)目t = o ，x ( f ) = x ，結(jié)束時(shí)間t f ，對(duì)每 一個(gè)k ，正= q = 0 ，幺= 0 ，o 】，計(jì)算，設(shè)o q r ：a o ( x ) ； ，- l 步驟4 如果t + r t f ，結(jié)束。否則，更新時(shí)