劉秋麗，楊 潔

(1．華南師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院，廣東廣州510631;2．北京信息科技大學(xué)理學(xué)院，北京100192)

布爾網(wǎng)絡(luò)和概率布爾網(wǎng)絡(luò)(Probabilistic Boolean Networks，簡稱PBN)能模擬真實世界的基因調(diào)控網(wǎng)絡(luò)［1-2］．PBN本質(zhì)上是傳統(tǒng)布爾網(wǎng)絡(luò)的概率推廣．布爾網(wǎng)絡(luò)B=(V，F(xiàn))是由基因節(jié)點的集合V={x1，x2，…，xn}，xi{0，1}(i=1，…，n)和函數(shù)序列 F={f1，…，fn}，fi:{0，1}n→{0，1}(i=1，…，n)組成的二元組．布爾函數(shù)序列F表示基因之間的調(diào)控準(zhǔn)則，函數(shù)fi是基因i的預(yù)測函數(shù)．假設(shè)基因i在時刻t的狀態(tài)記為xi(t)，xi(t)=0表示基因i沒有表達(dá)，而xi(t)=1表示基因i表達(dá)．在t時所有基因的表達(dá)可用狀態(tài)向量x(t)=(x1(t)，x2(t)，…，xn(t))表示．布爾網(wǎng)絡(luò)和它的動態(tài)都是確定的．但是基因表達(dá)是隨機(jī)的，為了彌補(bǔ)這種模型的不足，有必要將布爾網(wǎng)絡(luò)推廣到PBN．概率布爾網(wǎng)絡(luò)結(jié)合了基因之間的不確定性，因此經(jīng)常被用來模擬基因調(diào)控網(wǎng)的隨機(jī)性質(zhì)．給定一個概率布爾網(wǎng)絡(luò)，從一個狀態(tài)到其他狀態(tài)的轉(zhuǎn)移根據(jù)固定的轉(zhuǎn)移概率．PBN的動態(tài)性質(zhì)能用馬爾可夫鏈來描述，進(jìn)而Markov決策過程理論能用來分析和研究PBN中的最優(yōu)控制問題．

為了避免和疾病相聯(lián)系的狀態(tài)，利用最優(yōu)控制方法設(shè)計干預(yù)策略，從而影響網(wǎng)絡(luò)的動態(tài)演化［2-6］．本文考慮的優(yōu)化問題與文獻(xiàn)［2］-［6］不同，這導(dǎo)致我們的處理方法跟文獻(xiàn)［2］-［6］不同．PBN的最優(yōu)控制已在 Markov決策過程的框架下刻畫［2-3，5-6］，但是半Markov決策過程在PBN的應(yīng)用至今研究尚少．本文將利用半Markov決策過程的首達(dá)目標(biāo)理論研究PBN中的最優(yōu)控制問題，給出PBN的定義，建立PBN的半Markov決策過程模型，給出最優(yōu)問題的算法并用一個數(shù)值例子說明如何應(yīng)用本文所給的方法找到最優(yōu)治療方式．

1 概率布爾網(wǎng)絡(luò)

PBN是布爾網(wǎng)絡(luò)的推廣，本質(zhì)上是由布爾網(wǎng)絡(luò)組成的集類．布爾網(wǎng)絡(luò)中每個基因只有一個預(yù)測函數(shù)，和布爾網(wǎng)絡(luò)不同，在PBN中每個基因i的狀態(tài)由多個預(yù)測函數(shù)f(i)j(j=1，2，…，l(i))決定．預(yù)測函數(shù)f(i)j的選取概率為(i=1，2，…，n)．假設(shè)fj為第j個可能的網(wǎng)絡(luò)，則(1≤ji≤l(i)，i=1，2，…，n)．

這樣的網(wǎng)絡(luò)被選取的概率為

其中N=∏ni=1l(i)是布爾網(wǎng)絡(luò)的個數(shù)．

狀態(tài)(x1，…，xn)與(x'1，…，x'n)之間的轉(zhuǎn)移概率可由下式得到

為了表達(dá)方便，通過下列方程可以把二進(jìn)制數(shù)轉(zhuǎn)換為十進(jìn)制:

由于x(t)取值于{0，1}n，所以 zt取值從1到2n．由于x(t)到zt是一對一的，所以用十進(jìn)制表示和分量表示是等價的．

本文假設(shè)PBN兩次狀態(tài)轉(zhuǎn)移之間的時間間隔是隨機(jī)的．考慮時刻 tk和tk+1，對所有的tk≤t＜tk+1，PBN的基因狀態(tài)向量z(tk)=x．在tk+1時，網(wǎng)絡(luò)到達(dá)新的基因狀態(tài)向量z(tk+1)=y．設(shè)Xk+1表示轉(zhuǎn)移到狀態(tài)向量y之前停留在狀態(tài)向量x的時間，有Xk+1=tk+1-tk．時間長度Xk+1是隨機(jī)的，它的分布函數(shù)依賴當(dāng)前的狀態(tài)向量和將要到達(dá)的狀態(tài)向量．PBN的動態(tài)性質(zhì)可由半馬爾可夫過程刻畫．

2 概率布爾網(wǎng)絡(luò)的半M arkov決策過程模型

在這個部分把PBN描述為半Markov決策過程首達(dá)目標(biāo)模型．

首先給出PBN的狀態(tài)．在任意時刻t≥0，PBN的狀態(tài)x是基因狀態(tài)向量的十進(jìn)制表示．因此，狀態(tài)空間由所有可能的狀態(tài)向量組成，即S={1，2，…，2n}．給定PBN，一些狀態(tài)表示基因調(diào)控網(wǎng)的某些性質(zhì)，比如某些基因的表達(dá)與否導(dǎo)致網(wǎng)絡(luò)到達(dá)不想要狀態(tài)．根據(jù)想要狀態(tài)和不想要狀態(tài)的定義可以把狀態(tài)空間分為不同的種類．假設(shè)S0?S為想要狀態(tài)的集合，則S1=S-S0對應(yīng)不想要狀態(tài)的集合．

第二，在治療癌癥時諸如放射線和化學(xué)療法等外部變量(控制輸入)的使用可以使病人遠(yuǎn)離癌細(xì)胞擴(kuò)散的狀態(tài)．和PBN中狀態(tài)的二元表示一致，我們假設(shè)外部變量(控制輸入)為0或1．特別地假設(shè)n基因的PBN有m個控制輸入u1，u2，…，um．包含控制輸入向量v=(u1，…，um)能被看作行動．與狀態(tài)向量的表示一致，用控制輸入的十進(jìn)制a表示行動．A(x)表示當(dāng)PBN在狀態(tài)xS時行動的集合．換言之，A(x)表示在狀態(tài)x時對于控制者可選的控制輸入a的集合．本文中行動空間為有限集，即A(x)={1，2，…，2m}．

第三，引入半Markov核來描述PBN的動態(tài)演化．對任意的tRR+，x，yS，給定的控制輸入aA(x)，由文獻(xiàn)［7］，半Markov核Q為

第四，為了完成對PBN的半Markov決策過程首達(dá)目標(biāo)模型的描述，需要介紹一個優(yōu)化準(zhǔn)則．首先給出策略的定義．策略實際上是決策者選取行動的操作準(zhǔn)則．對所有的xS，策略為作用在S× RR上的函數(shù)g且滿足g(x，λ)屬于A(x)．換言之，函數(shù)g(x，λ)是當(dāng)前狀態(tài)為x時所選取的行動．設(shè)G為所有策略的集合．

設(shè)τ表示PBN首次到達(dá)不想要的狀態(tài)集S1的時間，即 τ=inf{t:ZtS1}．

至此，完成了PBN的半Markov決策過程首達(dá)目標(biāo)模型:

定義相應(yīng)的最優(yōu)值函數(shù)為策略g*G稱為最優(yōu)的，如果對所有的(x，λ)S0×R R有D(x，λ，g*)=D*(x，λ)．對于最優(yōu)策略g*，稱{g*(x，λ)}最優(yōu)行動．

3 優(yōu)化問題的算法

在上節(jié)建立的PBN的半Markov決策過程首達(dá)目標(biāo)模型的基礎(chǔ)上，本節(jié)主要給出最優(yōu)策略的計算方法．具體的算法如下:

由文獻(xiàn)［8］可知，最優(yōu)策略 g*存在．另外，假設(shè){Dn+1，n≥-1}滿足上述的值迭代算法，最優(yōu)值函數(shù)D*等于Dn+1，即對每對(x，λ)S0× RR，有D*(x，λ)=Dn+1(x，λ)．詳見文獻(xiàn)［8］．

4 數(shù)值模擬

考慮三基因x1，x2，x3的概率布爾網(wǎng)絡(luò)．基因x1由2個布爾函數(shù)f(1)1、f(1)2預(yù)測，基因x2由布爾函數(shù)f1(2)、f2(2)預(yù)測，基因x3由布爾函數(shù)f1(3)、f2(3)預(yù)測．假設(shè)x1取值為0或1的控制輸入．對于新的PBN，變量x1、x2和x3重新命名如下:變量x1為控制輸入a而x2和x3分別變?yōu)閤1和x2．由表1可知有4個布爾網(wǎng)絡(luò)(f1(1)，f1(2))、(f1(1)，f2(2))、(f2(1)，f1(2))以及(f(1)2，f(2)2)被選的概率P1、P2、P3和P4分別為P1=0．1、P2=0．4、P3=0．1 以及 P4=0．4．因此 c(1)1=0．15，c(1)2=0．85，c(2)1=0．2，c(2)2=0．8．

表1 真值表Table 1 Truth table

由于 n=2．因此變量 z可能取值于 1，2，3，4 中的任一個，狀態(tài)空間為S={1，2，3，4}．假設(shè)最想要的狀態(tài)是1、2和3，而不想要的狀態(tài)是4．所以S0={1，2，3}，S1={4}．因為 m=1，所以控制輸入 a 可能取1或者2，行動集A(x)={1，2}．報酬函數(shù)取值如下:r(1，1)=1，r(1，2)=2，r(2，1)=0．5，r(2，2)=0．8，r(3，1)=0．2，r(3，2)=0．6．

另外取ε=10-12．由上節(jié)給出的算法，一些數(shù)值結(jié)果見圖1．

圖1 函數(shù)Ta D*(x，λ)Figure 1 The function of Ta D*(x，λ)

由圖1，最優(yōu)策略g*為

［1］SHMULEVICH I，DOUGHERTY E R，KIM S，et al．Probabilistic Boolean networks:A rule-based uncertainty model for gene regulatory networks［J］．Bioinformatics，2002，18:261-274．

［2］DATTA A，BITTNER M L，DOUGHERTY E R．External control inmarkovian genetic regulatory networks［J］．Machine Learning，2003，52:169-191．

［3］LIU Q L，GUO X P，ZHOU T S，Optimal control for probabilistic Boolean networks［J］．IET Systems Biology，2010，4(2):99-107．

［4］CHING W，ZHANGhang SQ，JIAO Y，et al．Optimal control policy for probabilistic Boolean networkswith hard constraints［J］．IET Systems Biology，2009，3:90-99．

［5］PAL R，DATTA A，DOUGHERTY ER．Optimal infinite horizon control for probabilistic Boolean networks［J］．IEEE Transactions On Signal Processing，2006，52:169-191．

［6］CHEN PC Y，CHEN JW．A markovian approach to the control of genetic regulatory networks［J］．BioSystems，2007，90:535-545．

［7］GUO X P，HERNANDEZ-LERMA O．Continuous-time markov decision processes:Theory and applications［M］．New York:Springer，2009．

［8］HUANG Y H，GUO X P．Minimizing risk models in denumerable semi-Markov decision processes with a target set［C］∥Proceedings of the 29th Chinese control conference，Beijing，China，2010:1576-1581．