本發(fā)明屬于系統(tǒng)生物學(xué)技術(shù)領(lǐng)域，主要涉及生物信息學(xué)和生物數(shù)據(jù)挖掘，具體涉及一種基于模塊化因子圖的信號通路機(jī)制確認(rèn)方法。

背景技術(shù)：

多發(fā)性骨髓瘤(multiple myeloma，MM)是一種以惡性漿細(xì)胞克隆性增殖為特點的惡性腫瘤，是目前血液系統(tǒng)第二大惡性腫瘤。目前MM的治療主要包括傳統(tǒng)化療、新藥靶向治療及免疫治療等。雖然新的靶向治療明顯提高了MM的療效,但患者中位生存時間仍在3～5年,其發(fā)病機(jī)制尚不明確。因此,進(jìn)一步研究影響MM細(xì)胞生長的相關(guān)機(jī)制,尋建立模型研究方法，是亟需解決的問題。

無論從生物學(xué)還是從臨床表現(xiàn)看,MM細(xì)胞的特性不僅僅決定于其遺傳學(xué)特性(如染色體重排,缺失,擴(kuò)增或某些特定基因的突變)。相反,該疾病的病理生理學(xué)表現(xiàn)明顯受MM細(xì)胞與其所處的骨髓微環(huán)境間雙向相互作用的影響。Virginia Hughes指出,骨髓微環(huán)境對MM細(xì)胞的存活、生長以及耐藥等重要環(huán)節(jié)有著息息相關(guān)的作用。骨髓基質(zhì)細(xì)胞(bone marrow stromal cells,BMSCs)作為骨髓微環(huán)境的主要成員,與MM的發(fā)生、發(fā)展有著密切的關(guān)聯(lián)。

隨著對MM生物學(xué)研究的不斷深入,人們發(fā)現(xiàn)在腫瘤發(fā)生發(fā)展過程中，信號通路控制著眾多至關(guān)重要的細(xì)胞生物學(xué)過程。目前，對MM信號通路靶點的研究十分廣泛，主要靶向的通路包括這些信號通路包括PI3K/Akt/mTOR/P70S6K信號通路,IKK-αF/NF-κB信號通路,Ras/Raf/MAPK信號通路和JAK/STAT3通路,它們都可以通過以下途徑被激活:上游的細(xì)胞因子與相應(yīng)受體的結(jié)合,或通過粘附啟動的激酶途徑直接由細(xì)胞粘附誘導(dǎo)增殖、抗凋亡信號通路的激活。

眾所周知，傳統(tǒng)的生物實驗非常昂貴并且要花費大量的時間，所以近年來越來越多的人在用生物模型去模擬生物生長狀況，從模擬的層面上去分析藥物影響或者提取關(guān)鍵蛋白質(zhì)。Huiming Peng等人用系統(tǒng)生物學(xué)的方法研究p38MAPK異型的抗藥性，確定生物模型之后利用設(shè)置參數(shù)值的方式分別去探索p38的五種異型的抗藥性；Xiaoqiang Sun等人基于細(xì)胞內(nèi)的信號通路利用微分方程建模的方法研究在組織骨再生的過程中細(xì)胞因子的組合預(yù)測，對人體組織骨再生成時不同細(xì)胞因子對成骨細(xì)胞和破骨細(xì)胞的刺激作用進(jìn)行了探索，并且篩選出較好的細(xì)胞因子組合。但是，上面的研究方法存在一定的局限性，并沒有指出如何預(yù)測未知的致病因子的影響。

遺傳學(xué)改變和。生化條件引起的通路激活經(jīng)常發(fā)生在腫瘤惡變早期和進(jìn)展期，同時它們也是患者預(yù)后的重要指標(biāo)。系統(tǒng)生物學(xué)方法期望通過建立細(xì)胞信號轉(zhuǎn)導(dǎo)過程的模型,找到參與此過程的各種分子之間相互作用的網(wǎng)絡(luò),闡明其在基因調(diào)控、疾病發(fā)生中的作用。近幾年,對信號轉(zhuǎn)導(dǎo)網(wǎng)絡(luò)的定量分析逐日升溫,通常采用一系列的方程模型描述信號轉(zhuǎn)導(dǎo)通路的內(nèi)部變化過程。

現(xiàn)有技術(shù)中，采用的主要方案包括以下幾種：

(1)臨床實驗。在對于多發(fā)性骨髓瘤的研究中，目前絕大多數(shù)還是處于利用實驗方法去觀察腫瘤細(xì)胞的生長發(fā)育，尤其治療期間，基本上是靠醫(yī)生的經(jīng)驗去判斷。實驗成本比較昂貴。

(2)常微分方程(ordinary differential equations,ODE)。這是是描述動力學(xué)系統(tǒng)的常用方法，應(yīng)用微分方程組可以構(gòu)建一個復(fù)雜的數(shù)學(xué)模型，用以代表一系列生化反應(yīng)的相互作用模式，并且模擬生物系統(tǒng)中各組分的時序性動態(tài)變化。常微分方程(ODE)是質(zhì)量反應(yīng)動力學(xué)過程的數(shù)學(xué)代表，可以用來描述連續(xù)時間范圍內(nèi)生物系統(tǒng)各組分的動態(tài)變化。對于那些不考慮空間大小，并且反應(yīng)速度和反應(yīng)底物的濃度成一定的比例關(guān)系的生化反應(yīng)系統(tǒng)比較適用。Chen利用ODEs來描述ErbB信號通路的輸入輸出對細(xì)胞分化和增值的影響,文中用299個ODE方程表示828個級聯(lián)反映，共有229個參數(shù)，計算規(guī)模很大。

(3)Petri Net的不確定性、并行性、異步性，以及對分布式系統(tǒng)的描述和分析能力使其在描述生物系統(tǒng)特性時有很大的優(yōu)勢。Chen和等用Petri nets建立了鳥氨酸循環(huán)的代謝模型。

上述方案主要存在的缺陷有以下幾點：

(1)人為經(jīng)驗判斷，準(zhǔn)確率不高。

(2)ODE等數(shù)學(xué)方法僅僅是描述生物量的變化，并不能直接以圖形的方式展現(xiàn)生物系統(tǒng)的結(jié)構(gòu)特性。若可實現(xiàn)生物定量數(shù)據(jù)與圖形的結(jié)合與自動轉(zhuǎn)化，將能更好地刻畫生物系統(tǒng)結(jié)構(gòu)與動態(tài)性質(zhì)之間的關(guān)系。

(3)標(biāo)準(zhǔn)的Petri Net常用來定性分析生物網(wǎng)絡(luò)的結(jié)構(gòu)性質(zhì)，不能用于生物計算。

以下對本發(fā)明所涉及到的技術(shù)詞匯/技術(shù)術(shù)語注釋如下：

1、多發(fā)性骨髓瘤(multiple myeloma，MM)

2、骨髓基質(zhì)細(xì)胞(bone marrow stromal cells,BMSCs)

3、常微分方程(ordinary differential equations,ODE)

4、反向蛋白質(zhì)陣列(reverse phase protein arrays,RPPA)

技術(shù)實現(xiàn)要素：

有鑒于此，本發(fā)明在總結(jié)前人的研究基礎(chǔ)上，提出建立一個多層次的計算系統(tǒng)生物學(xué)模型來研究多發(fā)性骨髓瘤細(xì)胞的生長機(jī)制，利用現(xiàn)有RPPA數(shù)據(jù)，結(jié)合常微分方程組和Petri網(wǎng)來描述信號通路，并且對骨髓基質(zhì)細(xì)胞和腫瘤細(xì)胞通路之間相互作用的組合影響進(jìn)行量化的探索，并且使用模塊化的計算模型，降低計算成本。

具體而言，本發(fā)明所提出的技術(shù)方案如下：

本發(fā)明提供了一種基于模塊化因子圖的骨髓瘤信號通路機(jī)制確認(rèn)方法，該方法包括：

步驟1、獲取RPPA數(shù)據(jù)；

步驟2、對所述RPPA數(shù)據(jù)進(jìn)行預(yù)處理，粗粒度篩選關(guān)鍵蛋白質(zhì)；

步驟3、基于所述粗粒度篩選的關(guān)鍵蛋白質(zhì)，構(gòu)建在細(xì)胞剛性環(huán)境下，細(xì)胞內(nèi)蛋白質(zhì)的相互作用通路，形成信號通路；

步驟4、采用常微分方程描述所述信號通路，并將所述信號通路分解成多個小模塊，針對每個所述小模塊，進(jìn)行參數(shù)優(yōu)化，建立系統(tǒng)生物學(xué)模型；

步驟5、對所述系統(tǒng)生物學(xué)模型進(jìn)行參數(shù)分析，所述參數(shù)分析包括穩(wěn)定性分析和敏感性分析。

優(yōu)選地，所述步驟1中，RPPA數(shù)據(jù)的獲取，通過以下方式：

步驟1.1、用壓強(qiáng)為100pa和400pa的細(xì)胞膠體模擬正常細(xì)胞和腫瘤細(xì)胞，記錄不同時間點細(xì)胞內(nèi)蛋白質(zhì)的濃度；

步驟1.2、利用蛋白質(zhì)芯片，獲得正常細(xì)胞和腫瘤細(xì)胞的兩組RPPA數(shù)據(jù)。

優(yōu)選地，所述步驟2具體包括：

步驟2.1、對粗粒度篩選出的全部關(guān)鍵蛋白質(zhì)數(shù)據(jù)，以t＝0min為標(biāo)準(zhǔn)，進(jìn)行規(guī)范化，所述規(guī)范化方法為：

其中t₀表示t＝0min，表示第i個蛋白質(zhì)在t_j時刻的濃度，表示第i個蛋白質(zhì)在t₀時刻的濃度，為規(guī)范化后的蛋白質(zhì)濃度；

步驟2.2、計算正常細(xì)胞、腫瘤細(xì)胞內(nèi)蛋白質(zhì)濃度變化率，具體方式為：

將其中濃度變化大于50％的蛋白質(zhì)作為有意義的表達(dá)的蛋白質(zhì)，作為粗粒度篩選出的關(guān)鍵蛋白質(zhì)。

優(yōu)選地，所述步驟3具體包括：

步驟3.1、基于粗粒度篩選出的關(guān)鍵蛋白質(zhì)，通過IPA數(shù)據(jù)庫，搜索相互作用的通路；

步驟3.2、在所述步驟3.1中搜索出的通路中，選擇p≤0.05的通路，作為信號通路，其中p表示某蛋白質(zhì)在該通路中出現(xiàn)的誤差率。

優(yōu)選地，所述步驟4具體包括：

步驟4.1、使用常微分方程組描述信號通路，并使用RPPA數(shù)據(jù)中，高水平表達(dá)的蛋白質(zhì)在不同時間點的采樣數(shù)據(jù)，確定信號通路中的關(guān)鍵參數(shù)；使用Petri網(wǎng)描述整個信號通路；

步驟4.2、基于所述Petri網(wǎng)描述的整個信號通路，將整個信號通路分解成多個小模塊；使用粒子群優(yōu)化方法優(yōu)化各個所述小模塊參數(shù)，獲得相對較小的參數(shù)范圍，所述參數(shù)優(yōu)化的目標(biāo)函數(shù)是：

其中，表示蛋白質(zhì)濃度時間序列數(shù)據(jù)，表示通過常微分方程獲得的模擬的蛋白質(zhì)濃度時間序列，i表示蛋白質(zhì)索引,t_j表示時間點，Θ表示常微分方程中的參數(shù)，M表示蛋白質(zhì)數(shù)量，N表示時間點數(shù)量；

步驟4.3、把整個信號通路分解成兩個子網(wǎng)通路，并使用因子圖表示每個子網(wǎng)通路，為因子圖中的每個因子節(jié)點構(gòu)造適應(yīng)函數(shù)，所述適應(yīng)函數(shù)為：

步驟4.4、使用置信度傳播方法，調(diào)和所述兩個子網(wǎng)通路中共享的蛋白質(zhì)參數(shù)，并以所述步驟4.2中獲得的相對較小的參數(shù)范圍中的參數(shù)，作為置信度傳播方法的輸入?yún)?shù)，得到一個更小的參數(shù)范圍；

步驟4.5、以所述更小的參數(shù)范圍中的參數(shù)，作為粒子群優(yōu)化方法的輸入，對所述系統(tǒng)生物學(xué)模型進(jìn)行參數(shù)優(yōu)化，獲得最終的系統(tǒng)生物學(xué)模型。

優(yōu)選地，所述步驟4.2中，整個信號通路分解成多個小模塊，可依據(jù)如下分解原則：a.每個子模塊中蛋白質(zhì)數(shù)據(jù)盡可能少；b.每個子模塊中至少有一個蛋白質(zhì)濃度是有實際數(shù)據(jù)依據(jù)。

優(yōu)選地，所述步驟4.3中，分解成子網(wǎng)通路可依據(jù)如下規(guī)則：a.從細(xì)胞表型出發(fā)，依次找出促進(jìn)細(xì)胞增長或細(xì)胞死亡的蛋白質(zhì)；b.如果按兩條表型找出的兩類蛋白質(zhì)中共享蛋白質(zhì)數(shù)量超過90％，則把兩條子網(wǎng)合并為一個大網(wǎng)；c.如果其中一個大的子網(wǎng)的蛋白質(zhì)數(shù)量是另外一個子網(wǎng)蛋白質(zhì)數(shù)量的2倍或以上，則重新分解這個大的子網(wǎng)。

優(yōu)選地，所述粒子群優(yōu)化方法的具體方式如下：

v_i(t+1)＝wv_i(t)+c₁·rand()·(p_i(t)-x_i(t))+c₂·Rand()·(p_g(t)-x_i(t))

x_i(t+1)＝x_i+v_i(t+1)。

優(yōu)選地，所述穩(wěn)定性分析通過計算變異系數(shù)對參數(shù)進(jìn)行分析，所述變異系數(shù)的計算方法如下：

C·V＝(標(biāo)準(zhǔn)偏差SD/平均值Mean)×100％。

優(yōu)選地，所述敏感性分析中，參數(shù)的敏感性計算方法如下：

其中，ΔP_i是第i個參數(shù)值的變化量,表示一個很小的變化，例如可以是增加或減少1％，[ProteinName]表示系統(tǒng)輸出蛋白質(zhì)，例如Casp3、p90RSK、CyclinD1、p21、p7056k，P_i表示優(yōu)化的參數(shù)，s_i表示敏感性參數(shù)。

與現(xiàn)有技術(shù)相比，本發(fā)明技術(shù)方案具有以下的有益效果：

(1)基于腫瘤細(xì)胞信號通路的互相作用，用系統(tǒng)生物學(xué)的方法建立了計算模型模擬腫瘤細(xì)胞的增值和凋亡。

(2)用常微分方程去描述信號通路中的反應(yīng)，用petri網(wǎng)描述生物信號通路網(wǎng)絡(luò)結(jié)構(gòu)，實現(xiàn)生物定量數(shù)據(jù)與圖形的結(jié)合與自動轉(zhuǎn)化，更好地刻畫生物系統(tǒng)結(jié)構(gòu)與動態(tài)性質(zhì)之間的關(guān)系。然后將整個信號通路基于規(guī)則進(jìn)行模塊性的劃分，用尋優(yōu)算法獲取最優(yōu)值，降低了計算成本，提高了計算效率。

(3)通過模擬手段提高對腫瘤的生長過程中的分子機(jī)制的認(rèn)識，實現(xiàn)對轉(zhuǎn)移精確地進(jìn)行預(yù)測。

附圖說明

為了更清楚地說明本發(fā)明實施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對實施例或現(xiàn)有技術(shù)描述中所需要使用的附圖作簡單地介紹，顯而易見地，下面描述中的附圖僅僅是本發(fā)明的一些實施例，對于本領(lǐng)域普通技術(shù)人員來講，在不付出創(chuàng)造性勞動的前提下，還可以根據(jù)這些附圖獲得其它的附圖。

圖1為本發(fā)明實施例的方法流程圖；

圖2為本發(fā)明實施例的RPPA粗粒度篩選蛋白質(zhì)計算結(jié)果示例圖；

圖3為本發(fā)明實施例的構(gòu)建與細(xì)胞剛性相關(guān)的信號通路示例圖；

圖4為本發(fā)明實施例的最大限度簡化通路示例圖；

圖5為本發(fā)明實施例的用混合Petri網(wǎng)去描述信號通路示例圖；

圖6為本發(fā)明實施例的將信號通路分解成n個小模塊示例圖；

圖7為本發(fā)明實施例的信號通路分解成的子網(wǎng)通路一示例圖；

圖8為本發(fā)明實施例的信號通路分解成的子網(wǎng)通路二示例圖；

圖9為本發(fā)明實施例的參數(shù)穩(wěn)定性分析結(jié)果示例圖；

圖10為本發(fā)明實施例的參數(shù)敏感性分析結(jié)果示例圖；

圖11為本發(fā)明實施例的不同條件下蛋白質(zhì)對應(yīng)的參數(shù)變化示例圖。

具體實施方式

下面結(jié)合附圖對本發(fā)明實施例進(jìn)行詳細(xì)描述。應(yīng)當(dāng)明確，所描述的實施例僅是本發(fā)明一部分實施例，而不是全部的實施例?；诒景l(fā)明中的實施例，本領(lǐng)域普通技術(shù)人員在沒有作出創(chuàng)造性勞動前提下所獲得的所有其它實施例，都屬于本發(fā)明保護(hù)的范圍。

本領(lǐng)域技術(shù)人員應(yīng)當(dāng)知曉，下述具體實施例或具體實施方式，是本發(fā)明為進(jìn)一步解釋具體的發(fā)明內(nèi)容而列舉的一系列優(yōu)化的設(shè)置方式，而該些設(shè)置方式之間均是可以相互結(jié)合或者相互關(guān)聯(lián)使用的，除非在本發(fā)明明確提出了其中某些或某一具體實施例或?qū)嵤┓绞綗o法與其他的實施例或?qū)嵤┓绞竭M(jìn)行關(guān)聯(lián)設(shè)置或共同使用。同時，下述的具體實施例或?qū)嵤┓绞絻H作為最優(yōu)化的設(shè)置方式，而不作為限定本發(fā)明的保護(hù)范圍的理解。

本發(fā)明是通過找出影響骨髓瘤細(xì)胞生存的信號通路中的關(guān)鍵蛋白質(zhì)，輔助對于腫瘤細(xì)胞增生的模型建立方法研究。根據(jù)骨髓瘤細(xì)胞的剛性環(huán)境(骨髓瘤細(xì)胞壓強(qiáng)高)有助于腫瘤細(xì)胞的增生的結(jié)論，用壓強(qiáng)為100pa和400pa的細(xì)胞分別建立正常和腫瘤細(xì)胞的模型，測得得到兩組反相蛋白質(zhì)陣列數(shù)據(jù)，通過實驗數(shù)據(jù)(100pa VS 400pa)粗粒度篩選出可能影響腫瘤細(xì)胞增長的蛋白質(zhì)，并使用常微分方程組和Petri網(wǎng)描述信號相關(guān)的細(xì)胞信號通路，利用模塊化因子圖算法優(yōu)化模型中的關(guān)鍵參數(shù)，在最后通過合理和準(zhǔn)確性驗證，最終通過該模型得出敏感度較高的蛋白質(zhì)FAK、NFκB、mTOR1。圖1是本發(fā)明的總體流程圖，以下結(jié)合圖1對本發(fā)明的模型建立和計算方法進(jìn)行詳細(xì)闡述。

(1)獲取RPPA(Reverse Phase Protein Arrays)數(shù)據(jù)，建立模型基礎(chǔ)數(shù)據(jù)庫

獲取RPPA數(shù)據(jù)的方式可以是多樣的，例如，可以采用現(xiàn)有的公用數(shù)據(jù)庫或資料庫中提供的已有的骨髓瘤細(xì)胞比對RPPA數(shù)據(jù)，即反向蛋白質(zhì)陣列數(shù)據(jù)，其中對比數(shù)據(jù)模擬生長壓強(qiáng)需設(shè)置為100pa和400pa，以此作為模型建立的對比數(shù)據(jù)，從而建立數(shù)據(jù)庫?，F(xiàn)有技術(shù)中有多種途徑可以獲得上述RPPA公開數(shù)據(jù)，這里不再贅述。

也可以通過做模擬對比實驗，建立細(xì)胞生長的模型，設(shè)置例如壓強(qiáng)為100pa和400pa的細(xì)胞膠體模型模擬正常和腫瘤細(xì)胞，記錄在不同時間點細(xì)胞內(nèi)蛋白質(zhì)的濃度，在一個具體的實施方式中，該不同的時間點例如可以設(shè)置為0min,30min,60min,overnight，共4個時間點，當(dāng)然，也可以根據(jù)具體的算法需要，設(shè)定不同數(shù)量的采樣時間點。利用蛋白質(zhì)芯片得到兩組RPPA數(shù)據(jù)。該數(shù)據(jù)中包括在不同壓強(qiáng)下173種蛋白質(zhì)在不同時間點下的濃度。

當(dāng)然，此處也可以采用其他的已有技術(shù)或途徑，獲取可以作為后續(xù)模型算法所需要的骨髓瘤細(xì)胞的RPPA數(shù)據(jù)集。

(2)數(shù)據(jù)預(yù)處理，粗粒度篩選關(guān)鍵蛋白質(zhì)

對所有蛋白質(zhì)數(shù)據(jù)以t＝0min為標(biāo)準(zhǔn)對其規(guī)范化，其計算公式為：

其中t₀表示t＝0min，表示第i個蛋白質(zhì)在t_j時刻的濃度，表示第i個蛋白質(zhì)在t₀時刻的濃度，為規(guī)范化后的蛋白質(zhì)濃度。

計算兩種細(xì)胞內(nèi)蛋白質(zhì)濃度變化率，計算公式為：

根據(jù)此公式找出其中濃度變化大于50％的蛋白質(zhì)作為有意義的表達(dá)的蛋白質(zhì)。圖2是根據(jù)RPPA粗粒度篩選蛋白質(zhì)計算結(jié)果示例圖。

(3)構(gòu)建多發(fā)性骨髓瘤信號通路

信號通路是指能將細(xì)胞外的分子信號經(jīng)細(xì)胞膜傳入細(xì)胞內(nèi)發(fā)揮效應(yīng)的一系列酶促反應(yīng)通路。這些細(xì)胞外的分子信號(稱為配體，ligand)包括激素、生長因子、細(xì)胞因子、神經(jīng)遞質(zhì)以及其它小分子化合物等。細(xì)胞內(nèi)各種不同的生化反應(yīng)途徑都是由一系列不同的蛋白組成的，執(zhí)行著不同的生理生化功能。各個信號通路中上游蛋白對下游蛋白活性的調(diào)節(jié)(包括激活或抑制作用)主要是通過添加或去除磷酸基團(tuán)，從而改變下游蛋白的立體構(gòu)象完成的。所以，構(gòu)成信號通路的主要成員是蛋白激酶和磷酸酶，它們能夠快速改變和恢復(fù)下游蛋白的構(gòu)象。

在本發(fā)明中，根據(jù)步驟(2)粗粒度篩選的蛋白質(zhì)來構(gòu)造在細(xì)胞剛性環(huán)境下MIC和MSC細(xì)胞內(nèi)蛋白質(zhì)的相互作用通路。首先將篩選出的蛋白質(zhì)通過IPA(Ingenuity Pathway Analysis)數(shù)據(jù)庫尋找相互作用的重要的通路，IPA是基于云計算的一體化應(yīng)用軟件，它可以分析來源于基因組、microRNA、SNP、芯片、代謝組、蛋白組、RNA-Seq實驗以及各類小規(guī)模實驗數(shù)據(jù)。利用IPA，可以搜索基因、蛋白、化學(xué)分子、藥物的各類信息，并且?guī)椭鷺?gòu)建相互作用模型。在高度結(jié)構(gòu)化、集成豐富詳實生物化學(xué)知識的Ingenuity Knowledge Base支持下，IPA的分析和搜索可以幫助所獲得數(shù)據(jù)在生物體系中重要性。然后在搜索出的通路中選擇p≤0.05的通路，其中p值表示某蛋白質(zhì)在該pathway(即通路)中出現(xiàn)的誤差率。由于在找出的通路中有一些其他的蛋白質(zhì)，通過共享蛋白質(zhì)來整合所有細(xì)胞通路并參考MM細(xì)胞的相關(guān)文獻(xiàn)或現(xiàn)有的公開數(shù)據(jù)，來構(gòu)建與細(xì)胞剛性相關(guān)的信號通路，如圖3所示。由于估計由整個信號通路建立的常微分方程組中的所有參數(shù)是不可能的，因此我們需要重新定義信號通路，在保證通路結(jié)構(gòu)完整的情況下最大限度簡化通路，如圖4所示。

(4)建立系統(tǒng)生物學(xué)模型

系統(tǒng)生物學(xué)模型所包含生化反應(yīng)的動力學(xué)參數(shù)能夠明顯影響數(shù)學(xué)模擬結(jié)果，可以通過現(xiàn)有技術(shù)或相關(guān)文獻(xiàn)獲得，也可以通過體外的生理和生化實驗測定。但是，某些生化反應(yīng)的動力學(xué)參數(shù)無法直接獲得，需要通過實驗數(shù)據(jù)分析加以估計幾乎所有描述真實生物系統(tǒng)的模型都過于巨大，其內(nèi)部的動態(tài)變化過程無法被實驗數(shù)據(jù)全部描述。因此，系統(tǒng)生物學(xué)模型在不斷涌現(xiàn)的多實驗數(shù)據(jù)的整合分析中能夠起到關(guān)鍵的作用。

(4.1)在本發(fā)明中，我們使用常微分方程組描述信號通路，并采用RPPA試驗數(shù)據(jù)中高水平表達(dá)的蛋白質(zhì)在各個時間點的采樣數(shù)據(jù)來確定信號通路中的關(guān)鍵參數(shù)，接上例，設(shè)定該采樣時間點位4個。在一個具體的實施方式中，結(jié)合圖5至圖8，常微分方程可以通過以下方式進(jìn)行構(gòu)建：以pCylinD為例，結(jié)合圖5至圖8中可見，GSK3β促進(jìn)CylinD的磷酸化，而pp21抑制CylinD的磷酸化，根據(jù)這些相互作用可以得到：

其中[CyclinD]、[pp21]、[pGSK3β]分別表示蛋白質(zhì)CyclinD、pp21、pGSK3β的濃度；kCyclinD_pp21表示CylinD被pp21抑制的磷酸率、kCyclinD_pGSK3β表示CylinD被pGSK3β激活的磷酸率。

遵循這樣的規(guī)律，根據(jù)圖6,我們得到用來描述其它蛋白質(zhì)的40個常微分方程，該其它蛋白質(zhì)即如附圖6信號通路中涉及到的蛋白質(zhì)，由于其微分方程的表達(dá)式相同，此處不再一一贅述。由于這些常微分方程包含45未知參數(shù)并且蛋白質(zhì)之間有相互作用，使用一般的智能算法來確定參數(shù)值是相當(dāng)復(fù)雜的。本發(fā)明使用調(diào)整后的馬爾科夫鏈模塊分解方法把整個通路分解成多個小模塊，然后用智能算法來優(yōu)化每個小模塊的參數(shù)。該方法不僅可以減少計算機(jī)的壓力，也可以快速獲得更優(yōu)結(jié)果。

(4.2)采用信號通路的微分方程模型可以表征細(xì)胞在不同時間內(nèi)分子濃度的變化,觀察在不同外界條件下細(xì)胞的動態(tài)響應(yīng)過程,深入的揭示信號轉(zhuǎn)導(dǎo)的作用機(jī)制。同時，用混合Petri網(wǎng)去描述信號通路中的每一個生化反應(yīng)及級聯(lián)反應(yīng)，如圖5所示，把整個信號傳導(dǎo)通路基于Petri網(wǎng)和分解規(guī)則分解成n個小模塊，如圖6所示，應(yīng)用PSO算法(即粒子群優(yōu)化算法)對每個小模塊進(jìn)行估參(即參數(shù)估計)，從而在初始范圍搜索空間中獲得一個相對較小的參數(shù)范圍。

在一個具體的實施方式中，具體的分解規(guī)則可以設(shè)置為：a.每個子模塊中蛋白質(zhì)數(shù)據(jù)盡可能少；b.每個子模塊中至少有一個蛋白質(zhì)濃度是有實驗數(shù)據(jù)的，或者是有確定的數(shù)據(jù)支撐的，這樣才能保證計算的高效性。并使用粒子群算法來優(yōu)化各個小模塊的參數(shù)。

得到的蛋白質(zhì)濃度時間序列數(shù)據(jù)用表示，通過ODE方程(即常微分方程)獲得的模擬的蛋白質(zhì)濃度時間序列其中，i表示蛋白質(zhì)索引,t^j表示時間點，Θ表示ODE方程中的參數(shù)，參數(shù)優(yōu)化的目標(biāo)函數(shù)是：

其中，M表示蛋白質(zhì)數(shù)量，N表示時間點數(shù)量，最終得到一個參數(shù)初始范圍。

(4.3)然后把整個信號通路分解成兩個子網(wǎng)通路，如圖7、圖8所示，為每一個子通路構(gòu)造因子圖，并為每一個因子節(jié)點構(gòu)造適應(yīng)函數(shù)，即公式(1)。

在一個具體的實施方式中，可以為分解成子網(wǎng)通路設(shè)定規(guī)則：a.從細(xì)胞表型出發(fā)(細(xì)胞增長、細(xì)胞死亡)，依次找出促進(jìn)細(xì)胞增長或細(xì)胞死亡的蛋白質(zhì)；b.如果按兩條表型找出的兩類蛋白質(zhì)中共享蛋白質(zhì)數(shù)量超過90％，則把兩條子網(wǎng)合并為一個大網(wǎng)；c.如果其中一個大的子網(wǎng)的蛋白質(zhì)數(shù)量是另外一個子網(wǎng)蛋白質(zhì)數(shù)量的2倍或以上，則重新分解這個大的子網(wǎng)。另外，因子圖是指將一個具有多變量的全局函數(shù)因子分解，得到幾個局部函數(shù)的乘積，以此為基礎(chǔ)得到的一個雙向圖，叫做因子圖，因子圖是由變量節(jié)點和因子節(jié)點，以及連接兩個節(jié)點的邊構(gòu)成。

(4.4)應(yīng)用置信度傳播(Belief Propagation,BP)方法調(diào)和兩個子通路中共享的蛋白質(zhì)參數(shù)，解決兩個子網(wǎng)通路相同蛋白質(zhì)對應(yīng)的沖突的參數(shù)，從而得到一更優(yōu)的參數(shù)范圍，以(4.2)步中輸出的相對較小的參數(shù)作為BP算法的輸入?yún)?shù)，然后得到一個更小的參數(shù)范圍，最后把這個范圍作為最后應(yīng)用PSO算法進(jìn)行系統(tǒng)估參的輸入。從而實現(xiàn)整個模型形成一個不斷縮小搜索范圍的過程。

其中表示適應(yīng)度函數(shù)，表示模擬結(jié)果與實驗結(jié)果的誤差，分別表示參數(shù)集合和與因子節(jié)點對應(yīng)的分子濃度水平的集合。表示通過包括參數(shù)集Θ的ODEs方程計算得到的蛋白質(zhì)m在時間點t_j模擬濃度。表示蛋白質(zhì)m在時間點t_j的實驗水平的濃度。

(4.5)估參方法中，本發(fā)明采用粒子群優(yōu)化算法PSO，以步驟(4.4)中得到的參數(shù)范圍作為初始搜索空間繼續(xù)對參數(shù)進(jìn)行優(yōu)化，從而得到最優(yōu)值。

在一具體的實施方式中，可在該步中，重復(fù)運(yùn)行例如5遍或更多遍，保證算法的穩(wěn)定性。

PSO中，每個優(yōu)化問題的解都是搜索空間中的一個粒子。所有的粒子都有一個由被優(yōu)化的函數(shù)決定的適應(yīng)值(fitness value)，每個粒子還有一個速度決定他們飛翔的方向和距離。在本發(fā)明中，PSO初始化為一群隨機(jī)粒子(隨機(jī)解)。然后通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個"極值"來更新自己。第一個就是粒子本身所找到的最優(yōu)解，這個解叫做個體極值。另一個極值是整個種群目前找到的最優(yōu)解，這個極值是全局極值。這個過程一直迭代，直到找到比較滿意的解為止。此解作為整個模型估參的最優(yōu)解。

(5)參數(shù)分析

a.穩(wěn)定性分析

某些生物系統(tǒng)(如代謝網(wǎng)絡(luò)和信號通路等)具有特殊的潛能，能夠聚合于特定的平衡態(tài)(或稱為穩(wěn)定狀態(tài))。當(dāng)生物系統(tǒng)達(dá)到穩(wěn)定狀態(tài)時，系統(tǒng)內(nèi)所有組分濃度改變速率為零。穩(wěn)態(tài)分析可以鑒定并合理描述這些穩(wěn)定狀態(tài)，這有助于闡釋生物系統(tǒng)內(nèi)部的穩(wěn)態(tài)維持機(jī)制，以及細(xì)胞凋亡響應(yīng)等生物學(xué)過程所涉及的狀態(tài)快速轉(zhuǎn)換機(jī)制。在本發(fā)明中，穩(wěn)定性分析指分析參數(shù)是否穩(wěn)定，對于某個參數(shù)，如果其變異系數(shù)大于1，則認(rèn)為不穩(wěn)定，否則就是穩(wěn)定的。變異系數(shù)的計算公式為：C·V＝(標(biāo)準(zhǔn)偏差SD/平均值Mean)×100％

根據(jù)此公式，對本文中的所有參數(shù)進(jìn)行穩(wěn)定性分析。如圖9所示。其中有10個參數(shù)的變異系數(shù)大于1，即所占比例為22％，意味著在這個模型里將近80％的參數(shù)是穩(wěn)定的。

b.敏感性分析

在系統(tǒng)生物學(xué)模型構(gòu)建完成之后，簡單的結(jié)構(gòu)化和參數(shù)化并不能很好地解釋生物系統(tǒng)的內(nèi)在調(diào)控機(jī)制，因此對模型特性的整體分析是非常重要的。敏感性分析(Sensitivity Analysis)是指從定量分析的角度研究在一定范圍內(nèi)擾動系統(tǒng)組分的初始濃度或反應(yīng)動力學(xué)參數(shù)，對模型輸出結(jié)果影響的大小(特定組分的濃度變化或模型整體狀態(tài)的改變)。生物學(xué)實驗上可以理解為敲除或過表達(dá)特定基因產(chǎn)物，之后測定生物體或細(xì)胞特定表型的變化。敏感性分析能夠闡釋系統(tǒng)輸出對特定參數(shù)值的依賴性，這通常是研究者最為感興趣的。值得注意的是，從理論分析的角度進(jìn)行敏感性分析時，初始濃度或反應(yīng)參數(shù)的擾動范圍可以隨意選取。但是從實踐的視角出發(fā)，如果期望獲得的敏感性分析結(jié)果是可靠的，那么擾動的選取范圍必須符合真實的生理生化狀態(tài)。此外，敏感性分析方法對于尋找生物網(wǎng)絡(luò)的關(guān)鍵調(diào)控位點具有重要的指導(dǎo)意義。

在本發(fā)明中，敏感性分析是一種衡量某一特定參數(shù)對輸出的影響，分析參數(shù)的敏感性，即參數(shù)變化是否會引起整個系統(tǒng)輸出值的變化，簡單來說，給每一個參數(shù)增加或減少1％，觀察對系統(tǒng)輸出的影響，從而找出敏感的蛋白質(zhì)。根據(jù)信號通路可知輸出為蛋白質(zhì)Casp3、p90RSK、CyclinD1、p21、p7056k濃度。參數(shù)的敏感性計算公式：

其中，ΔP_i是第i個參數(shù)值的變化連,表示一個很小的變化(增加或減少1％)，[ProteinName]表示系統(tǒng)輸出蛋白質(zhì)Casp3、p90RSK、CyclinD1、p21、p7056k，P_i表示優(yōu)化的參數(shù)，s_i表示敏感性參數(shù)。其中，每個參數(shù)通過增加或減少1％，觀察Casp3、p90RSK、CyclinD1、p21、p7056k變化百分比，示例如圖10，通過分析可得系統(tǒng)輸出變化較大的參數(shù)占所有參數(shù)的比例不超過2％，該結(jié)果表示所有參數(shù)中只有5-7個參數(shù)比較敏感，這說明系統(tǒng)輸出蛋白質(zhì)的濃度變化受這幾個比較敏感的參數(shù)影響。

綜合上述結(jié)果，對比正常細(xì)胞和腫瘤細(xì)胞對應(yīng)的參數(shù)變化值，如圖11所示，由上圖可以看出幾種蛋白質(zhì)濃度變化范圍比較大，則認(rèn)為FAK、NFκB、mTOR1這幾種蛋白質(zhì)可能是通路中對腫瘤細(xì)胞生長或增生影響比較大的，從而實現(xiàn)對蛋白質(zhì)的篩選。

以上所述，僅為本發(fā)明的具體實施方式，但本發(fā)明的保護(hù)范圍并不局限于此，任何熟悉本技術(shù)領(lǐng)域的技術(shù)人員在本發(fā)明揭露的技術(shù)范圍內(nèi)，可輕易想到的變化或替換，都應(yīng)涵蓋在本發(fā)明的保護(hù)范圍之內(nèi)。因此，本發(fā)明的保護(hù)范圍應(yīng)以權(quán)利要求的保護(hù)范圍為準(zhǔn)。