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

上海大學(xué)碩士學(xué)位論文 摘要 基因調(diào)控網(wǎng)絡(luò)是由一組基因、蛋白質(zhì)、r n a 、小分子以及它們之間的相互 作用所構(gòu)成的一種生化系統(tǒng)。對基因調(diào)控網(wǎng)絡(luò)的算法研究在理論和實踐上有著 重要的意義。 本文首先給出化學(xué)主方程，并介紹幾種精確隨機模擬算法和快速近似算 法。其次，對剛性生化系統(tǒng)的混合模擬算法進行了介紹，并給出多尺度隨機模 擬算法，然后在其基礎(chǔ)上發(fā)展了改進多尺度隨機模擬算法。數(shù)值模擬結(jié)果證實 了改進多尺度隨機模擬算法的優(yōu)點。最后介紹時滯生化系統(tǒng)的精確模擬算法， 并提出一種模擬時滯生物化學(xué)反應(yīng)系統(tǒng)的自適應(yīng)f l e a p 算法( d a r l e a p ) 。 該算法能夠較大地提高模擬速度。數(shù)值實驗證實了d a r l e a p 算法的優(yōu)點。 關(guān)鍵詞：基因調(diào)控網(wǎng)絡(luò)；多尺度算法；改進多尺度算法；時滯生化系統(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)下進行的研究工作。 除了文中特別加以標注和致謝的地方外，論文中不包含其他人已發(fā) 表或撰寫過的研究成果。參與同一工作的其他同志對本研究所做的 任何貢獻均已在論文中作了明確的說明并表示了謝意。 簽名：j 亟垂蘧日期：衛(wèi)邋 本論文使用授權(quán)說明 本人完全了解上海大學(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)，力求闡述作為一個系統(tǒng)的生物系統(tǒng)。系統(tǒng)生物學(xué)是在細胞、組織、器官 和生物體整體水平上研究結(jié)構(gòu)和功能各異分子之間的相互作用，并通過計算機 模擬來定量地描述和預(yù)測生物功能、表型和行為。系統(tǒng)生物學(xué)將在基因組序列 的基礎(chǔ)上完成由生命密碼到生命過程的研究。這是一個逐步整合的過程：由對 生物體內(nèi)各種分子的鑒別及其相互作用的認識到對途徑、網(wǎng)絡(luò)、模塊的研究并 最終完成整個生命活動的路線圖。系統(tǒng)生物學(xué)的研究大致有以下四個方面乜3 1 ： 1 ) 系統(tǒng)結(jié)構(gòu)( s y s t e ms t r u c t u r e s ) 的識別和研究這包括基因相互作用、生 化代謝途徑以及這些相互作用以何種機制調(diào)節(jié)生物系統(tǒng)的研究； 2 ) 系統(tǒng)的動態(tài)特征( s y s t e md y n a m i c s ) 研究生物系統(tǒng)在不同條件下的行 為，我們可以對系統(tǒng)進行代謝分析( 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 ) 、動態(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)的特定行為的機理，分叉分析又 可以研究系統(tǒng)在多維空間中行為的時域變化( t i m e - v a r y i n gc h a n g e s ) ，其中每一 維代表系統(tǒng)中一個生化成分的特定濃度； 3 ) 系統(tǒng)的控制方法( t h ec o n t r o lm e t h o d ) 系統(tǒng)控制一個細胞的狀態(tài)能夠 使得細胞功能受損最小，從而為疾病治療提供潛在的藥物靶標； 4 ) 系統(tǒng)的設(shè)計方法( t h ed e s i g nm e t h o d ) 將原有的生物系統(tǒng)修改和重建成 具有所需功能的系統(tǒng)，其方法和策略可以從確定的設(shè)計原則、仿真模擬來發(fā)展， 這樣可以避免盲目的反復(fù)試驗( 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)與開發(fā)等方面將具有巨大 應(yīng)用潛力。系統(tǒng)生物學(xué)還可應(yīng)用于傳統(tǒng)的生物技術(shù)產(chǎn)業(yè)，開拓能源生物技術(shù)、 上海人學(xué)碩士學(xué)位論文 材料生物技術(shù)和環(huán)境生物技術(shù)等新領(lǐng)域。系統(tǒng)生物學(xué)的研究將不僅對生命科學(xué) 和生物技術(shù)產(chǎn)生重大影響，而且對整個國民經(jīng)濟、社會和人類本身產(chǎn)生重大影 響。 1 2 基因調(diào)控網(wǎng)絡(luò)的動力學(xué) 系統(tǒng)生物學(xué)的研究內(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 ) 的動力學(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)控 開了基因調(diào)控網(wǎng)絡(luò)研究的先 聲。由于細胞中的眾多組件可以通過復(fù)雜的途徑相互作用，所以大多數(shù)基因調(diào) 控網(wǎng)絡(luò)龐大而且復(fù)雜，如圖1 1 中所示的一個基因調(diào)控網(wǎng)絡(luò)圖晦1 。因此對于如 此龐大并且復(fù)雜的調(diào)控過程的動力學(xué)性質(zhì)，僅僅靠觀察和描述的方法是很難理 解的，因此需要引入數(shù)學(xué)工具通過對其建立模型來進行分析和模擬。隨著計算 機的普及，利用計算機進行分析和模擬也日益普遍。 牛g 骶2 圖1 1 基因調(diào)控網(wǎng)絡(luò)圖 生命的基本單元一細胞是由許許多多不同種類的生物大分子所組成，其中 2 上海大學(xué)碩士學(xué)位論文 每一種生物大分子，如一種m r n a 或一種蛋白質(zhì)，在細胞內(nèi)都存在著若干個相同 的分子。例如，參與細胞構(gòu)造的蛋白質(zhì)如肌動蛋白( a c t i n s ) 的分子數(shù)目通常 都比較大，可能達到數(shù)百萬個，而負責(zé)細胞活動調(diào)控的如轉(zhuǎn)錄因子則分子數(shù)目 一般較少，可能只有幾十個。但不管是多還是少，在細胞內(nèi)不存在只有單個分 子的生物大分子種類。也就是說，過去人們從定性角度談?wù)摰哪衬郴蚧蚰衬?蛋白質(zhì)，實際上都是某某基因或某某蛋白質(zhì)的若干個分子的集合體。 這種生物分子集合體的行為與無機世界里的分子集合體的行為有著巨大的 不同，其動力學(xué)性質(zhì)也有著明顯的區(qū)別。例如，大氣中的各種氣體分子的運動， 可以用玻爾茲曼方程等統(tǒng)計物理學(xué)的理論進行描述，而生物分子集合體則不會 遵守同樣的方程。我們知道，新陳代謝是生命存在的必要形式，所有的生物分 子都有著一定的壽命，在生命的每時每刻，總是有一些分子在“死亡( 降解) ， 有一些分子在“誕生一( 合成) 。因此，在細胞內(nèi)的每一種分子集合體，都有若 干個分子不停地被降解，若干個分子不停地被合成。除了分子數(shù)量的波動外， 這種集合體內(nèi)每個分子的活性也不會是均一的，有的高，有的低。這種在數(shù)量 和活性上均具有隨機波動的分子，在參與細胞的活動時會引入一種擾動，稱為 外在噪聲( e x t r i n s i cn o i s e ) 。另一方面，即使在分子的數(shù)量和活性完全一致 的細胞群體中，反應(yīng)的速率和進行的時間在不同的細胞中仍然存在著隨機差異。 這種隨機性質(zhì)被稱為內(nèi)在噪聲( i n t r i n s i cn o i s e ) 。現(xiàn)有的研究表明訓(xùn)，隨機 噪聲廣泛地存在于細胞的d n a 復(fù)制、基因轉(zhuǎn)錄和蛋白質(zhì)合成。 噪聲的存在對細胞的正?；顒佑兄匾挠绊?，它有不利的一面。比如會影 響細胞的正常生理過程，對基因表達造成影響，產(chǎn)生突變，引起疾病。但它也 有有利的一面。例如，生物體內(nèi)許多活動都具有節(jié)律性，也就是通常所說的“生 物鐘 研究人員發(fā)現(xiàn)，噪聲在生物鐘的調(diào)控方面扮演著重要的角色。此外，在 真核生物中，噪聲導(dǎo)致基因表達的不一致可以作為一種手段，用來參與細胞分 化的控制。當前，對生命過程中的隨機性的認識及其對生命活動影響的研究正 在成為一個新的科學(xué)前沿口 引州1 1 幻。更重要的是，由于在細胞內(nèi)的生命活動過 程中存在著各種隨機性，使得細胞的特性和行為表現(xiàn)出種復(fù)雜性和多樣性， 這使得生命體在應(yīng)付各種環(huán)境的變化上具有更強的生存機會。這一點不只是體 上海大學(xué)碩上學(xué)位論文 現(xiàn)在整個物種群體應(yīng)對環(huán)境的變化上，對于單個細胞也是如此。 當然，在細胞內(nèi)，不同的噪聲源所引起的噪聲對細胞的影響是不一樣的。 同樣，細胞內(nèi)也存在著不同的機制可以對噪聲進行抑制、放大甚至加以利用。 由于在生命過程中噪聲無處不在，它不僅是生命過程的一部分，也為我們探索 一些生命活動提供了一個途徑。這些噪聲不僅僅與參與生化反應(yīng)的各種分子數(shù) 目有關(guān)，也與濃度、溫度等因素有關(guān)。從通訊的角度來看，這些噪聲也是一種 通訊手段，可以傳達各種信息，不但使若干個反應(yīng)之間的聯(lián)系變得更為微妙， 也在細胞間的通訊上起著重要的作用。因此，研究生命系統(tǒng)中的隨機行為，對 于更好地了解生命過程以及在計算機中重構(gòu)細胞和生命方面有著重要的意義。 另外由于基因調(diào)控網(wǎng)絡(luò)中生物化學(xué)反應(yīng)的時間多尺度性，如d n a 與蛋白質(zhì) 結(jié)合以及蛋白質(zhì)聚合等快速反應(yīng)，轉(zhuǎn)錄翻譯和降解等慢速反應(yīng)，以及時滯特性 ( t i m ed e l a y ) 等使得我們研究基因調(diào)控網(wǎng)絡(luò)中的動力學(xué)性質(zhì)變得尤為困難。 已有一些這方面的進展 1 4 - 1 9 o 但這些理論都基于這樣一個假設(shè)前提：所有的生 物化學(xué)反應(yīng)過程均為服從泊松統(tǒng)計的馬爾科夫過程，即不考慮現(xiàn)實中的時滯特 性。而在實際的基因調(diào)控網(wǎng)絡(luò)中，轉(zhuǎn)錄和翻譯等反應(yīng)過程不僅相對較慢，而且 常常由多個基元反應(yīng)組成，這種過程服從帶有時滯的高斯統(tǒng)計，研究這種時滯 的調(diào)控網(wǎng)絡(luò)的解析方法和數(shù)值模擬方法也就尤為重要乜1 刪。 1 3 本文的工作及安排 本文重點研究基因調(diào)控網(wǎng)絡(luò)的隨機模擬算法，給出其數(shù)學(xué)模型以及常見的 數(shù)值模擬算法晦刎。第二章介紹連續(xù)確定性模型和離散隨機性模型，給出幾種精 確的蒙特卡羅算法和快速近似算法。針對生化反應(yīng)系統(tǒng)的多尺度現(xiàn)象，第三章 介紹幾種有效克服多尺度的算法，給出c a o 等人提出的多尺度隨機算法，并且 在其基礎(chǔ)上發(fā)展了改進多尺度隨機算法。該算法比多尺度算法能更有效地模擬 剛性系統(tǒng)。第四章介紹時滯基因調(diào)控網(wǎng)絡(luò)的精確時滯隨機模擬算法。針對模擬 復(fù)雜的時滯系統(tǒng)的低效性，提出了模擬時滯生化系統(tǒng)的自適應(yīng)f l e a p 算法 ( d a r l e a p ) 。該算法能夠自動調(diào)節(jié)f ，使得在時間區(qū)間【f ，t + r ) 內(nèi)發(fā)生多次反 4 上海人學(xué)碩士學(xué)位論文 應(yīng)事件并且能夠精確地滿足l e a p 條件。數(shù)值試驗證實了d a r l e a p 算法能夠有 效地模擬復(fù)雜的時滯生化網(wǎng)絡(luò)。 論文安排如下： 本章闡述了系統(tǒng)生物學(xué)的背景，簡要介紹了基因調(diào)控網(wǎng)絡(luò)的相關(guān)背景，闡 述了本文的研究目的和意義，并介紹了目前國內(nèi)外在這方面的研究進展。 第二章分別介紹了基因調(diào)控網(wǎng)絡(luò)的兩種動力學(xué)模型：連續(xù)確定性模型和離 散隨機模型以及常用的算法。 第三章介紹剛性生化網(wǎng)絡(luò)的特性和一些常用克服剛性的隨機算法。給出多 尺度算法，并且在其基礎(chǔ)上提出了改進多尺度算法。 第四章介紹時滯生化反應(yīng)系統(tǒng)的精確算法并且提出了o a r l e a p 算法。 第五章對本文進行總結(jié)與展望。 上海大學(xué)碩? l 學(xué)位論文 第二章基因調(diào)控網(wǎng)絡(luò)的動力學(xué)模型及其模擬算法 從計算角度出發(fā)，系統(tǒng)生物學(xué)中基因調(diào)控網(wǎng)絡(luò)的研究大致可以分為以下兩 條路線：一支是基于數(shù)據(jù)挖掘的基因調(diào)控網(wǎng)絡(luò)的重構(gòu)，該領(lǐng)域主要是從大量的 試驗數(shù)據(jù)( 如相互作用數(shù)據(jù)) 中提取信息，如從g e n b a n k 、p d b 等數(shù)據(jù)庫中采用 神經(jīng)網(wǎng)絡(luò)、隱馬氏模型、支持向量機等方法重構(gòu)基因調(diào)控網(wǎng)絡(luò)，這些方法都是 一種靜態(tài)的研究手段，主要研究基因調(diào)控網(wǎng)絡(luò)的拓撲性質(zhì)。而另一支研究路線 則是從整體出發(fā)研究基因調(diào)控網(wǎng)絡(luò)( 如反應(yīng)通道、單細胞和多細胞通信等) 的 動力學(xué)性質(zhì)，這是一種動態(tài)的研究方法。因此需要建立設(shè)計能夠準確描述細胞 內(nèi)生物化學(xué)網(wǎng)絡(luò)的動力學(xué)模型。生化網(wǎng)絡(luò)中的化學(xué)反應(yīng)呈現(xiàn)遠離熱力學(xué)平衡的 現(xiàn)象啪瑚1 ，并且分子復(fù)制數(shù)目有時候非常低，化學(xué)反應(yīng)具有明顯的隨機性。但 是當分子數(shù)目巨大時，隨機不明顯可以忽略不計，可以用連續(xù)的確定的微分方 程來描述。一般的生物化學(xué)反應(yīng)系統(tǒng)中參加反應(yīng)的分子數(shù)目通常不是很大，而 隨機性則在反應(yīng)中起著重要的作用，因此可用化學(xué)主方程刻畫其動力學(xué)行為。 本文主要討論基因調(diào)控網(wǎng)絡(luò)中的動力學(xué)性質(zhì)的數(shù)值模擬。由于對基因調(diào)控 網(wǎng)絡(luò)的刻畫的不同，其數(shù)值模擬大致可以分為以下2 種不同的方法：( 1 ) 以連 續(xù)的方式處理，變量為反應(yīng)分子的濃度，表示為確定性方程( 組) ，在2 1 小節(jié) 給出；( 2 ) 以離散的方式處理，變量為反應(yīng)分子的個數(shù)，反應(yīng)過程表示為隨機 過程，這一方法在處理所含分子數(shù)目較低的系統(tǒng)時就顯得尤為重要，是本文研 究的重點，在剩余的小節(jié)中給出。 2 1 連續(xù)確定性動力學(xué)模型 當生物反應(yīng)系統(tǒng)的反應(yīng)分子的數(shù)目足夠大，其隨機性可以忽略時，其數(shù)目 可以認為是與時間有關(guān)的連續(xù)變量。我們可以用反應(yīng)物的濃度作為變量的微分 方程來刻畫基因調(diào)控網(wǎng)絡(luò)的動力學(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)分子的個數(shù)，q ( t ) 為反應(yīng)物濃度。求解方程組( 2 i 1 ) 可以得 到調(diào)控網(wǎng)絡(luò)中各反應(yīng)分子的濃度隨時間的變化情流。關(guān)于微分方程( 組) 的數(shù) 值方法相對比較成熟，如歐拉方法( e u l e r ) 、龍格一庫塔方法( r u n g e - k u t t a ) 虛盤 守。 由于生物化學(xué)反應(yīng)的時間多尺度性，反應(yīng)系統(tǒng)中經(jīng)常存在有些反應(yīng)很慢， 而有些反應(yīng)則很快。在準平衡態(tài)假設(shè)下，原微分方程組( 2 1 1 ) 簡化成由微分 方程和代數(shù)方程組成的耦合系統(tǒng)，通過求解這種耦合系統(tǒng)，我們同樣可以得到 調(diào)控網(wǎng)絡(luò)中各反應(yīng)分子的濃度隨時陽l 的變化情況。 2 2 離散隨機動力學(xué)模型 當細胞中的生物化學(xué)反應(yīng)的隨機性不能忽略時，連續(xù)確定性系統(tǒng)不能正確地 預(yù)測生物化學(xué)反應(yīng)系統(tǒng)的動力學(xué)性質(zhì)。因此必須需要能描述其隨機性的精確數(shù) 學(xué)模型。本小節(jié)首先介紹生物化學(xué)反應(yīng)網(wǎng)絡(luò)，然后給出其動力學(xué)模型。 考慮溫度恒定、體積不變和空間均勻的非時滯生物化學(xué)反應(yīng)系統(tǒng)。該系統(tǒng) 具有種化學(xué)分子s = s l ，s ) ，通過m 個化學(xué)反應(yīng)通道r = r l ，一，r m ) 作 用。系統(tǒng)在時刻t 的狀態(tài)向量為x ( t ) = ( x l ( f ) ，x ( f ”，其中x f ( f ) 為第i 種化 學(xué)分子在時刻t 的分子數(shù)目，f = l ，n 。收= ( v 小，) 為通道尺i 發(fā)生一次 反應(yīng)的分子改變向量，其中k = l ，m 。這m 個生物化學(xué)反應(yīng)形成一個相互作 用網(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 為隨機反應(yīng)速率常數(shù)，它只與分子的物理性質(zhì)和系統(tǒng)的反應(yīng) 溫度以及體積有關(guān)。若記僅為反應(yīng)通道盈中反應(yīng)分子的表達式，a k ( x ( f ) ) 為通 道尺i 的反應(yīng)傾向函數(shù)，則有表達式吼( x ( f ) ) - c a ( x ( t ) ) 。此式的含義將在后面 給出。反應(yīng)通道咒的反應(yīng)類型如下表所示： 表2 1 反應(yīng)通道的類型以及吃的表達形式 其中五為分子最的分子數(shù)目。 給定x ( j ) = 石，在無窮小時間區(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 ) ) 出+ 口( 凼) 并且兩個及以上的通道發(fā)生反應(yīng)的概率為d ( 凼) 。從而有下式成立 p ( x ，f + 出) = p ( x ) + 口，( x 一吩娥x 一_ ，t ) d t - 巳( x ) 尸( x ，t ) d t + o ( d t ) 當d t o 時，這就變威化學(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 的直接求解非常困難，因為方程的維數(shù)等于可能的分子狀態(tài)數(shù)。于是 我們通過逼近方法來研究主方程的性質(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é)位論文 對方程進行泰勒展開，然后從二階導(dǎo)數(shù)后截斷，從而得到不同狀態(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 ) 對狀態(tài)截斷來逼近c m e 。對于c m e 的維數(shù)不是 很大時，f s p 方法能夠提供精確的解。當c m e 的維數(shù)很大時，f s p 方法能夠提供 保證精度的解。雖然該算法能夠有效地解一些簡單的生物化學(xué)反應(yīng)系統(tǒng)，但是 對于復(fù)雜的生物化學(xué)反應(yīng)系統(tǒng)仍然是不可解的。許多蒙特卡羅算法被提出用來 克服解c m e 的困難。下面兩小節(jié)給出一些重要的算法。 2 3 精確隨機模擬算法 為了模擬生物化學(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ā)展了改進下一反應(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 ) 。這些隨機算法的解等價于c m e 的解，因而是精 確隨機模擬算法。我們首先給出4 種算法，然后對4 種算法做出點評。 對于任何精確隨機模擬算法，在每一個時間步長上，系統(tǒng)只有一個狀態(tài)。 在給定狀態(tài)x ( t ) = 工下，必須通過回答下述兩個問題完成反應(yīng)系統(tǒng)的模擬：1 ) 什么時候發(fā)生下一個反應(yīng)? 2 ) 下一個反應(yīng)是哪個通道? 具體地，若f 與_ ，分別 表示下一反應(yīng)發(fā)生的時間和發(fā)生反應(yīng)的通道，精確隨機模擬算法通過描述下述 事件e 而完成模擬整個過程： e ：時間區(qū)間 t , t + r ) 內(nèi)不發(fā)生任何反應(yīng)，且r ，將在無窮小區(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 ) ，分別對_ ，求和以及對f 積分，可以得到下一反應(yīng)尺，發(fā)生的概 9 上海大學(xué)碩士學(xué)位論文 率分布和下一反應(yīng)發(fā)生的時間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)獨立產(chǎn)生2 個隨機數(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)初始時刻= 0 ，結(jié)束時間0 ； 步驟2 根據(jù)系統(tǒng)狀態(tài)x ( f ) 計算每個反應(yīng)的傾向函數(shù)( = l ，m )a o ； 步驟3 從u ( o ，1 ) 產(chǎn)生隨機數(shù)吒，匕； 步驟4 根據(jù)式( 2 3 2 ) 和( 2 3 3 ) 計算下一反應(yīng)通道及其發(fā)生時間間 隔a t ； 步驟5 更新反應(yīng)時間和系統(tǒng)狀態(tài)：t = t + a t ，x = x + 匕； 步驟6 如果f 砭( 瓦) ) 。假定系統(tǒng)的狀態(tài)為x ( f ) = 工時，假設(shè) 知最，互。通道r 下一次發(fā)生的時間間隔為似：墨號孚。選擇最小的時間間隔 吼l 叫 出，- - - m z a , , ，則系統(tǒng)下一反應(yīng)的時間間隔為越，發(fā)生反應(yīng)的通道為r 。更 新系統(tǒng)： f = f + 馘互= 正+ 口( 石) 她和x = 工+ 和對于通道r 更新 只= 露+ l n ( 1 r ) ，其中廠u ( o ，1 ) 。重復(fù)上述過程，得到m n r m 算法。 刖r m 算法如下： 步驟1 初始化系統(tǒng)狀態(tài)，初始反應(yīng)時刻乞= o ，結(jié)束時間為f ， 最= 互= o ，七= 1 ，m ； 步驟2 根據(jù)系統(tǒng)狀態(tài)石( f ) 計算每個反應(yīng)的傾向函數(shù)q ( 后= 1 ，m ) ； 1 2 上海大學(xué)碩士學(xué)位論文 步驟3 從u ( o ，1 ) 產(chǎn)生m 個獨立變量氣； 步驟4 對每一個后，設(shè)最：蚍馬； ，七 步驟5 對每一個后，設(shè) ：墨二互； “七 步驟6 i 發(fā)a t i = 卿 ，更新系統(tǒng)：t = t + a t ，瓦= 瓦+ 口( x ) a t 。和z = 石+ 和對于通道局更新只= 只+ 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 是最快的，因為n r m 和m n r m 雖然每次迭代只需用一個隨機數(shù)，但是它查找 最小時間，也需要消耗一定的c p u 時間，具體討論參見。 3 2 。 2 4 快速的隨機模擬算法 上述算法為精確模擬算法，由于精確地模擬每一個反應(yīng)事件導(dǎo)致計算量很 大，因而對于模擬復(fù)雜的生化系統(tǒng)是低效的。為了加快模擬的速度，g i l l e s p i e 等人提出了f l e a p 算法棚1 。在給定時間區(qū)間 t ，f + f ) 內(nèi)系統(tǒng)狀態(tài)改變很小 時，該算法利用泊松變量生成每一個反應(yīng)通道在該區(qū)間內(nèi)的反應(yīng)次數(shù)。該類算 法能在給定的時間區(qū)間【f ，f + r ) 上發(fā)生多次反應(yīng)，故對于分子數(shù)目很大的系統(tǒng)是 比較有效的。對某類多通道分子反應(yīng)系統(tǒng)，在整個反應(yīng)過程中，其中一些反應(yīng) 通道發(fā)生反應(yīng)的頻率與另外一些反應(yīng)通道發(fā)生反應(yīng)的頻率相比，要大幾個數(shù)量 級，并且一些分子參加反應(yīng)頻率比另外一些分子參加反應(yīng)的頻率要大得多。系 統(tǒng)的狀態(tài)在短時間內(nèi)，主要由反應(yīng)頻率高的反應(yīng)通道刻畫，在長時間內(nèi)，主要 由反應(yīng)頻率低的反應(yīng)通道來刻畫。這種現(xiàn)象叫做多時間尺度現(xiàn)象( 簡稱多尺度 現(xiàn)象) ，在數(shù)學(xué)上稱做“剛性 。許多算法被提出用來克服系統(tǒng)的剛性，將在下 一章給出。這一小節(jié)，我們介紹模擬分子數(shù)目多的系統(tǒng)的算法。首先我們介紹 g i l l e s p i e 的f l e a p 算法，該類算法允許反應(yīng)通道在給定時間區(qū)間【f ，f + f ) 內(nèi) 上海大學(xué)碩士學(xué)位論文 發(fā)生多次反應(yīng)，但必須滿足l e a p 條件：系統(tǒng)的反應(yīng)傾向函數(shù)幾乎不變。精確 的表述如下：給定時間區(qū)間 f ，f + f ) 和系統(tǒng)狀態(tài)x ( f ) = 工，r s 通道反應(yīng)的次數(shù)由 泊松隨機變量巧( r ；功= p o i s s i o n ( a ，( x ) f ) 逼近。在時間區(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 利用上述選取時間步長_ r 的方法，我們可以得到r - l e a p 算法 步驟l 初始化系統(tǒng)狀態(tài)，初始反應(yīng)時刻= o ，結(jié)束時間為0 步驟2 在時刻t 系統(tǒng)狀態(tài)為x ，計算所有反應(yīng)的傾向函數(shù)以及它們的和 肼 口o ( x ) = 乃； 步驟3 根據(jù)式( 2 4 3 ) 計算f ； 步驟4 如果步驟2 中得到的時間步長f 小于l 和某個數(shù)( 不妨假設(shè)是1 0 的乘積，則放棄這個時間步長，轉(zhuǎn)而接連執(zhí)行一定個數(shù)( 不妨設(shè)為 1 0 0 ) 的精確隨機模擬( s s a ) 過程，然后再嘗試獲取新的時間步長f 。 1 4 上海大學(xué)碩上學(xué)位論文 反之，如果步驟2 中得到的時間步長f 大于1 a o 和那個數(shù)的乘積， 則接受這個時間步長r ，向下進行到步驟4 ； 步驟5 對于每個反應(yīng)_ ，( j = l ，m ) ，產(chǎn)生一個均值為a i ( 工) f 的p o i s s o n 隨 機變量k ，代表反應(yīng)，在時間步長f 內(nèi)發(fā)生的次數(shù)； 步驟6 更新f = f + f ，x = x q - 巧_ ； 步驟7 如果t t ，返回到步驟2 循環(huán)，否則結(jié)束。 當m ，很大時，式( 2 4 5 ) 和( 2 4 6 ) 的計算量為o ( m 2 + 2 m 2 ) ，為 了有效地降低其計算量，c a o 等人3 7 1 提出了對于f l e a p 算法的有效的步長選 擇算法。該算法是基于反應(yīng)分子的相對變化的。對于上述的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)分子集合，詳細請參考文獻 3 7 ，( 2 4 8 ) 式的計算量為2 m 2 。故 當很大時，能夠有效地降低計算量，提高速度。 上述的算法都是假定在給定時間區(qū)間【t ，f + f ) 內(nèi)系統(tǒng)狀態(tài)改變很小時，該算 法利用泊松隨機數(shù)確定每一個反應(yīng)通道在該區(qū)間內(nèi)的反應(yīng)次數(shù)。由于泊松數(shù)的 無界性，有可能產(chǎn)生的反應(yīng)次數(shù)導(dǎo)致某些分子變成負的，但是負的分子數(shù)不符 合現(xiàn)實。為了避免負分子數(shù)的出現(xiàn)，t i a n 等人口1 和c h a t t e r j e e 等人嘲獨立地提 出了二項f l e a p 算法。因為二項分布產(chǎn)生的隨機數(shù)是有界的，因此二項 f l e a p 算法分布可以避免負分子數(shù)目的產(chǎn)生。c a i 等人翰1 提出的k - l e a p 算 法能夠減少產(chǎn)生負分子數(shù)目的概率。c a o 等人洶1 發(fā)展了一種避免負分子數(shù)產(chǎn)生 的算法。該算法將反應(yīng)通道分成容易產(chǎn)生負分子數(shù)的關(guān)鍵通道和不容易產(chǎn)生負 分子數(shù)的非關(guān)鍵通道。用標準的g i l l e s p i e 隨機算法模擬關(guān)鍵通道，用f l e a p 算法模擬非關(guān)鍵通道。這樣大大地降低了產(chǎn)生負分子數(shù)的概率，但該類算法在 確定了時間步長f 后，就不能再進行調(diào)整了，屬于先驗f l e a p 算法。最近d a v i d 1 5 上海人學(xué)碩士學(xué)位論文 提出了將后驗的l e a p 檢查點植入f l e a p 算法1 ，稱之為后驗r - l e a p 算法。 該算法利用單位速率泊松過程可以自動地調(diào)節(jié)時間步長f ，使得精確地滿足 l e a p 條件，從而產(chǎn)生負分子數(shù)的概率為o ，從而可以保證精度。下面我們介紹 后驗f l e a p 算法。 單位速率泊松過程k 具有以下兩種性質(zhì)： a ) 對于r i m 互瓦( f ) ，增量砭( 乃) 一砭( 五) = k ( 互一石) 是與系統(tǒng)狀態(tài)x ( f ) 是獨立 的參數(shù)為互一互的泊松分布； b ) 給定o s 彳 露= 互。對于步長f ，通道r 的內(nèi)部時間為疋+ 吼( x ( f ) ) r 時， 發(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 乞) 則在時間區(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 弦按升序插入到已存儲的 1 6 蔓 紫 上海大學(xué)碩士學(xué)位論文 信息里。用更小的f 重復(fù)上述過程，直到找到符合條件的f 。 后驗f - l e a p 算法： 步驟l 初始化：初始時間及初始分子數(shù)目t = o ，x ( f ) = x ，結(jié)束時間t f ，對每 一個k ，正= q = 0 ，幺= 0 ，o 】，計算，設(shè)o q r ：a o ( x ) ； ，- l 步驟4 如果t + r t f ，結(jié)束。否則，更新時