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

背景技術(shù)：

多發(fā)性骨髓瘤(multiple myeloma，MM)是一種以惡性漿細(xì)胞克隆性增殖為特點(diǎn)的惡性腫瘤，是目前血液系統(tǒng)第二大惡性腫瘤。目前MM的治療主要包括傳統(tǒng)化療、新藥靶向治療及免疫治療等。雖然新的靶向治療明顯提高了MM的療效,但患者中位生存時(shí)間仍在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)境對(duì)MM細(xì)胞的存活、生長以及耐藥等重要環(huán)節(jié)有著息息相關(guān)的作用。骨髓基質(zhì)細(xì)胞(bone marrow stromal cells,BMSCs)作為骨髓微環(huán)境的主要成員,與MM的發(fā)生、發(fā)展有著密切的關(guān)聯(lián)。

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

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

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

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

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

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

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

上述方案主要存在的缺陷有以下幾點(diǎn)：

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

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

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

以下對(duì)本發(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ù)實(shí)現(xiàn)要素：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟2.2、計(jì)算正常細(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的通路，作為信號(hào)通路，其中p表示某蛋白質(zhì)在該通路中出現(xiàn)的誤差率。

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

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

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

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

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

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

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

優(yōu)選地，所述步驟4.2中，整個(gè)信號(hào)通路分解成多個(gè)小模塊，可依據(jù)如下分解原則：a.每個(gè)子模塊中蛋白質(zhì)數(shù)據(jù)盡可能少；b.每個(gè)子模塊中至少有一個(gè)蛋白質(zhì)濃度是有實(shí)際數(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)合并為一個(gè)大網(wǎng)；c.如果其中一個(gè)大的子網(wǎng)的蛋白質(zhì)數(shù)量是另外一個(gè)子網(wǎng)蛋白質(zhì)數(shù)量的2倍或以上，則重新分解這個(gè)大的子網(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)定性分析通過計(jì)算變異系數(shù)對(duì)參數(shù)進(jìn)行分析，所述變異系數(shù)的計(jì)算方法如下：

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

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

其中，ΔP_i是第i個(gè)參數(shù)值的變化量,表示一個(gè)很小的變化，例如可以是增加或減少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ì)胞信號(hào)通路的互相作用，用系統(tǒng)生物學(xué)的方法建立了計(jì)算模型模擬腫瘤細(xì)胞的增值和凋亡。

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

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

附圖說明

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

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

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

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

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

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

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

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

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

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

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

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

具體實(shí)施方式

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

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

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

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

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

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

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

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

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

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

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

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

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

信號(hào)通路是指能將細(xì)胞外的分子信號(hào)經(jīng)細(xì)胞膜傳入細(xì)胞內(nèi)發(fā)揮效應(yīng)的一系列酶促反應(yīng)通路。這些細(xì)胞外的分子信號(hào)(稱為配體，ligand)包括激素、生長因子、細(xì)胞因子、神經(jīng)遞質(zhì)以及其它小分子化合物等。細(xì)胞內(nèi)各種不同的生化反應(yīng)途徑都是由一系列不同的蛋白組成的，執(zhí)行著不同的生理生化功能。各個(gè)信號(hào)通路中上游蛋白對(duì)下游蛋白活性的調(diào)節(jié)(包括激活或抑制作用)主要是通過添加或去除磷酸基團(tuán)，從而改變下游蛋白的立體構(gòu)象完成的。所以，構(gòu)成信號(hào)通路的主要成員是蛋白激酶和磷酸酶，它們能夠快速改變和恢復(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是基于云計(jì)算的一體化應(yīng)用軟件，它可以分析來源于基因組、microRNA、SNP、芯片、代謝組、蛋白組、RNA-Seq實(shí)驗(yàn)以及各類小規(guī)模實(shí)驗(yàn)數(shù)據(jù)。利用IPA，可以搜索基因、蛋白、化學(xué)分子、藥物的各類信息，并且?guī)椭鷺?gòu)建相互作用模型。在高度結(jié)構(gòu)化、集成豐富詳實(shí)生物化學(xué)知識(shí)的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)的信號(hào)通路，如圖3所示。由于估計(jì)由整個(gè)信號(hào)通路建立的常微分方程組中的所有參數(shù)是不可能的，因此我們需要重新定義信號(hào)通路，在保證通路結(jié)構(gòu)完整的情況下最大限度簡化通路，如圖4所示。

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

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

(4.1)在本發(fā)明中，我們使用常微分方程組描述信號(hào)通路，并采用RPPA試驗(yàn)數(shù)據(jù)中高水平表達(dá)的蛋白質(zhì)在各個(gè)時(shí)間點(diǎn)的采樣數(shù)據(jù)來確定信號(hào)通路中的關(guān)鍵參數(shù)，接上例，設(shè)定該采樣時(shí)間點(diǎn)位4個(gè)。在一個(gè)具體的實(shí)施方式中，結(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個(gè)常微分方程，該其它蛋白質(zhì)即如附圖6信號(hào)通路中涉及到的蛋白質(zhì)，由于其微分方程的表達(dá)式相同，此處不再一一贅述。由于這些常微分方程包含45未知參數(shù)并且蛋白質(zhì)之間有相互作用，使用一般的智能算法來確定參數(shù)值是相當(dāng)復(fù)雜的。本發(fā)明使用調(diào)整后的馬爾科夫鏈模塊分解方法把整個(gè)通路分解成多個(gè)小模塊，然后用智能算法來優(yōu)化每個(gè)小模塊的參數(shù)。該方法不僅可以減少計(jì)算機(jī)的壓力，也可以快速獲得更優(yōu)結(jié)果。

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

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

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

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

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

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

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

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

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

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

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

(5)參數(shù)分析

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

某些生物系統(tǒng)(如代謝網(wǎng)絡(luò)和信號(hào)通路等)具有特殊的潛能，能夠聚合于特定的平衡態(tài)(或稱為穩(wěn)定狀態(tài))。當(dāng)生物系統(tǒng)達(dá)到穩(wěn)定狀態(tài)時(shí)，系統(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)定，對(duì)于某個(gè)參數(shù)，如果其變異系數(shù)大于1，則認(rèn)為不穩(wěn)定，否則就是穩(wěn)定的。變異系數(shù)的計(jì)算公式為：C·V＝(標(biāo)準(zhǔn)偏差SD/平均值Mean)×100％

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

b.敏感性分析

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

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

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

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

以上所述，僅為本發(fā)明的具體實(shí)施方式，但本發(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)。