SIR-DNA傳染病動(dòng)力學(xué)優(yōu)化算法*

陸秋琴，黃光球

西安建筑科技大學(xué) 管理學(xué)院，西安 710055

1 引言

DNA計(jì)算[1]是利用DNA雙螺旋結(jié)構(gòu)和堿基互補(bǔ)配對規(guī)律進(jìn)行信息編碼，將要運(yùn)算的對象映射成DNA分子鏈，通過生物酶的作用，生成各種數(shù)據(jù)池，再按照一定的規(guī)則將原始問題的數(shù)據(jù)運(yùn)算高度并行地映射成DNA分子鏈的可控的生化反應(yīng)過程，然后利用分子生物技術(shù)，檢測所需要的運(yùn)算結(jié)果[2]。

DNA計(jì)算包括分子內(nèi)DNA計(jì)算[3]、分子間DNA計(jì)算[1]和超分子DNA計(jì)算[2]。分子內(nèi)DNA計(jì)算借助于分子內(nèi)的形態(tài)轉(zhuǎn)移操作，用單DNA分子構(gòu)建可編程的狀態(tài)機(jī)[3]；分子間DNA計(jì)算集中在不同DNA分子間的雜交反應(yīng)，使其作為計(jì)算中的一個(gè)基本步驟[1]；超分子DNA計(jì)算是利用不同序列的原始DNA分子的自裝配過程進(jìn)行計(jì)算[4-5]。

DNA計(jì)算利用DNA反應(yīng)的強(qiáng)大并行計(jì)算能力，已成功地解決了許多NP難題。吳帆等人[4]提出了一種基于DNA自組裝模型來求解N皇后問題的DNA計(jì)算方法，該算法降低了生化解的錯(cuò)誤率。李肯立等人[5]通過引入DNA自組裝模型，提出了一種求解最大團(tuán)問題的DNA計(jì)算方法，該算法給出了DNA分子的編碼方案及結(jié)果檢測的實(shí)驗(yàn)方法。姚慶安等人[6]提出一種基于k-臂分子和粘貼計(jì)算求解最短路徑問題的DNA計(jì)算模型。徐京雷等人[7]提出了一種基于改進(jìn)的閉環(huán)模型求解TSP問題的DNA算法，該算法能在較低的時(shí)間復(fù)雜度內(nèi)解決TSP問題。馬瑩等人[8]將著色問題轉(zhuǎn)化為求最大獨(dú)立集問題，然后通過質(zhì)粒DNA分子生物實(shí)驗(yàn)得到最大獨(dú)立集，證明了該質(zhì)粒DNA算法有效并且是可行的。吳雪等人[9]針對最大匹配問題，應(yīng)用DNA粘貼計(jì)算模型尋求該問題最優(yōu)解的生物操作過程，實(shí)例表明該算法是可行的。

目前，DNA計(jì)算的大量研究還停留在紙面上，很多設(shè)想和方案沒有條件付諸實(shí)驗(yàn)[2]，對維數(shù)較高的優(yōu)化問題的求解還存在一些困難。為解決此問題，本文從一個(gè)新的角度提出了一種新的DNA算法。本文算法依據(jù)一種傳染病在一個(gè)生態(tài)系統(tǒng)中傳播的場景構(gòu)造而得，該生態(tài)系統(tǒng)中動(dòng)物被感染后，傳染病攻擊該動(dòng)物的某些DNA片斷，從而導(dǎo)致該動(dòng)物具有類似遺傳病的缺陷。本文算法不需依賴生物實(shí)驗(yàn)即可進(jìn)行優(yōu)化計(jì)算。

在現(xiàn)實(shí)世界，人群中經(jīng)常爆發(fā)的很多傳染病均可用KM倉室建模方法[10-11]描述。傳染病在個(gè)體之間的傳播，使得個(gè)體之間的相互作用關(guān)系體現(xiàn)得淋漓盡致；一種傳染病攻擊的是人類個(gè)體的少部分器官，將該現(xiàn)象映射到對優(yōu)化問題的求解，就是每次處理的變量個(gè)數(shù)只是全部變量的極少部分。因此，該算法具有求解高維問題的天然優(yōu)勢。為簡單起見，本文利用SIR疾病傳播場景來說明本文的思路，所提出的優(yōu)化算法稱為SIR-DNA算法。

2 SIR-DNA算法設(shè)計(jì)

考慮優(yōu)化問題：

式中，Rn是n維歐氏空間；X=(x1,x2,…,xn)是一個(gè)n維決策向量；H為搜索空間，又稱解空間；f(X)為目標(biāo)函數(shù)；gi(X)≥0為第i個(gè)不等式約束條件，i=1,2,…,I，I為不等式約束條件個(gè)數(shù)。目標(biāo)函數(shù)f(X)和約束條件gi(X)不需要特殊的限制條件。

把優(yōu)化模型式（1）中約束條件所形成的可行解空間設(shè)想成一個(gè)生態(tài)系統(tǒng)ES。假設(shè)ES包含有N個(gè)動(dòng)物，每個(gè)動(dòng)物稱為一個(gè)個(gè)體，這些個(gè)體的編號是1,2,…,N。有一種傳染病Z在該生態(tài)系統(tǒng)中傳播，動(dòng)物被感染后，傳染病Z攻擊該動(dòng)物的某些DNA片斷，從而導(dǎo)致該動(dòng)物具有類似遺傳病的缺陷。在現(xiàn)實(shí)世界，HIV、天花、SARS、寨卡等傳染病毒就是這樣的病毒。

2.1 動(dòng)物個(gè)體DNA分子

動(dòng)物的每條染色體攜帶一個(gè)DNA分子，DNA是由分別帶有4種堿基A、T、C、G的脫氧核苷酸鏈接組成的雙螺旋長鏈分子。在這條雙螺旋的長鏈中，共有約30億個(gè)堿基對，而基因則是DNA長鏈中有遺傳效應(yīng)的一些片段。在組成DNA的數(shù)量浩瀚的堿基對中，有一些特定位置的單個(gè)核苷酸經(jīng)常發(fā)生變異引起DNA的多態(tài)性，人們稱之為位點(diǎn)。染色體、基因和位點(diǎn)的結(jié)構(gòu)關(guān)系如圖1所示。在DNA長鏈中，位點(diǎn)個(gè)數(shù)約為堿基對個(gè)數(shù)的1/1 000。由于位點(diǎn)在DNA長鏈中出現(xiàn)頻繁，多態(tài)性豐富，近年來成為人們研究DNA遺傳信息的重要載體。大量研究表明，動(dòng)物的許多表型性狀差異以及對傳染病毒的易感性等都可能與某些位點(diǎn)相關(guān)聯(lián)，或和包含有多個(gè)位點(diǎn)的基因相關(guān)聯(lián)。

Fig.1 Structure of chromosomes,genes and loci圖1 染色體、基因和位點(diǎn)的結(jié)構(gòu)關(guān)系

對于ES中的每個(gè)動(dòng)物個(gè)體，采用A、T、C、G的編碼方式來獲取每個(gè)位點(diǎn)的信息，因?yàn)槿旧w具有雙螺旋結(jié)構(gòu)，所以用兩個(gè)堿基的組合表示一個(gè)位點(diǎn)的信息，如表1所示。表1中，在位點(diǎn)rs100015位置，不同樣本的編碼都是T和C的組合，有3種不同編碼方式TT、TC和CC。類似地，其他位點(diǎn)雖然堿基的組合不同，但也只有3種不同編碼。

Table 1 Loci name and loci alleles of DNAfragments表1 DNA片段位點(diǎn)名稱和位點(diǎn)等位基因信息

在ES中，傳染病Z攻擊的是動(dòng)物個(gè)體W個(gè)特殊基因中某些位點(diǎn)，這些特殊基因稱為致病基因，而這W個(gè)致病基因均由R個(gè)位點(diǎn)組成，其編號為1,2,…,R；每個(gè)位點(diǎn)均由3個(gè)堿基對組成。對于不同的動(dòng)物個(gè)體，傳染病Z攻擊的是W個(gè)致病基因中的哪些基因完全是隨機(jī)的；若某個(gè)致病基因被攻擊，組成該基因的R個(gè)位點(diǎn)中的哪些位點(diǎn)受到攻擊也完全是隨機(jī)的。染病個(gè)體治愈后會(huì)獲得短期免疫力，與免疫相關(guān)的特殊基因也有W個(gè)，這些特殊基因稱為免疫基因，每個(gè)免疫基因也包含R個(gè)位點(diǎn)；愈后個(gè)體在哪些免疫基因中的哪些位點(diǎn)獲得免疫，也是完全隨機(jī)的；雖然個(gè)體在不同免疫基因的不同位點(diǎn)上獲得免疫力，但是免疫效果都是一樣。為計(jì)算處理方便，將A、T、C、G采用二進(jìn)制編碼分別表示為00、01、10、11。

2.2 SIR模型轉(zhuǎn)換成個(gè)體類別轉(zhuǎn)換方程

將ES中的動(dòng)物群分成易感者類（S）、染病者類（I）和治愈者類（R）。易感個(gè)體一旦與染病者接觸，就會(huì)因感染而變成染病者；染病者經(jīng)過治療后會(huì)成為治愈者；治愈者能夠獲得暫時(shí)的免疫，但免疫喪失后又會(huì)變成易感者。假設(shè)該傳染病不是惡性的，故在一個(gè)時(shí)間段內(nèi)，該傳染病不會(huì)造成動(dòng)物死亡；此外，不考慮該時(shí)間段內(nèi)個(gè)體的出生、死亡、遷入和遷出。故可假設(shè)該ES中的動(dòng)物個(gè)體數(shù)N為常數(shù)。若傳染病的發(fā)生率是飽和接觸率，則基于KM經(jīng)典的假設(shè)[10-11]，該生態(tài)系統(tǒng)可用SIR傳染病模型描述[12]，如圖2所示，其數(shù)學(xué)模型為式（2）所示。

式中，β為動(dòng)物種群飽和接觸傳染率；γ為治愈系數(shù)；α為免疫喪失率；S(t)、I(t)、R(t)分別為時(shí)期t易感者、染病者和治愈者的占比。

Fig.2 Bin structure of SIR epidemic model圖2 SIR傳染病倉室結(jié)構(gòu)

因S(t)+I(t)+R(t)=1，S(t)、I(t)和R(t)實(shí)際上描述了一個(gè)概率分布，故對任意一個(gè)個(gè)體來說，S(t)、I(t)、R(t)相當(dāng)于該個(gè)體屬于易感者類、染病者類和治愈者類的概率，不妨將其分別稱為易感概率、染病概率和治愈概率。因此，式（2）實(shí)際上是計(jì)算任意一個(gè)個(gè)體在時(shí)期t的屬于易感者類、染病者類和治愈者類的概率的方程。為快速計(jì)算起見，將式（2）寫成遞推形式，并注意到式（2）中有一個(gè)方程是多余的，可以去掉，故對于個(gè)體i來說，有：

在時(shí)期t，采用式（3）計(jì)算出個(gè)體i的Si(t)、Ii(t)和Ri(t)。顯然，在某個(gè)時(shí)期，任何一個(gè)個(gè)體只能處于S類、I類和R類中的某一個(gè)類。于是，個(gè)體i在時(shí)期t處于S類、I類和R類3個(gè)類別中的哪個(gè)類，可由Si(t)、Ii(t)和Ri(t)所決定的概率分布確定。表2中所列情況符合圖2描述的SIR模型的類別轉(zhuǎn)換情形。

Table 2 Legal class transfer of SIR epidemic model表2 SIR傳染病模型的合法類別轉(zhuǎn)換

既然優(yōu)化問題式（1）的可行解空間與ES相對應(yīng)，那么該ES中一個(gè)個(gè)體對應(yīng)于優(yōu)化問題式（1）的一個(gè)可行解，N個(gè)個(gè)體所對應(yīng)的可行解集就是：

個(gè)體i的特征j對應(yīng)于可行解Xi中的一個(gè)變量xij，j=1,2,…,K+n。因此，個(gè)體i與可行解Xi是等價(jià)概念。個(gè)體的體質(zhì)強(qiáng)弱用體質(zhì)指數(shù)IPI來表示，IPI指數(shù)對應(yīng)于優(yōu)化問題式（1）的目標(biāo)函數(shù)值。好的可行解對應(yīng)具有較高IPI值的個(gè)體，即體質(zhì)強(qiáng)壯的個(gè)體，差的可行解對應(yīng)具有較低IPI值的個(gè)體，即體質(zhì)虛弱的個(gè)體。個(gè)體i的IPI指數(shù)計(jì)算方法為：

2.3 基因表達(dá)設(shè)計(jì)

在ES中，因動(dòng)物個(gè)體的種類相同，故對于不同的個(gè)體，傳染病毒Z攻擊的W個(gè)特殊基因是相似的（如圖3（a）所示），但不會(huì)相同，該特殊基因組稱為致病基因組；這W個(gè)致病基因被攻擊后發(fā)生改變，就會(huì)使動(dòng)物患病。類似地，對于同種動(dòng)物個(gè)體，能夠?qū)魅静《綵免疫的W個(gè)免疫基因是類似的，但不會(huì)相同，該基因組稱為免疫基因組，如圖3（b）所示；對于同種動(dòng)物中的不同個(gè)體，傳染病被治愈后，這W個(gè)免疫基因會(huì)發(fā)生變化，會(huì)使動(dòng)物獲得免疫。免疫基因組和致病基因組是不同序列的基因組。

因動(dòng)物個(gè)體之間存在差異，故傳染病毒攻擊致病基因組中的哪些致病基因，以及若一個(gè)致病基因被攻擊，則該致病基因中的哪些位點(diǎn)發(fā)生改變，則完全因人而異。同理，當(dāng)動(dòng)物個(gè)體染病后被治愈并獲得短期免疫力，但該個(gè)體在免疫基因組中的哪些免疫基因，以及若一個(gè)免疫基因獲得免疫力，則該免疫基因中的哪些位點(diǎn)發(fā)生改變，完全因人而異。

傳染病毒Z攻擊動(dòng)物個(gè)體的致病基因組和免疫基因組序列分別為：

式（5）、（6）中的基因Di、Vi，i=1,2,…,W，均由R個(gè)位點(diǎn)組成，即：

因每個(gè)動(dòng)物個(gè)體是不相同的，故組成每個(gè)基因的位點(diǎn)是不相同的，從而對不同的動(dòng)物個(gè)體來說，式（7）中的取值是不相同的，每個(gè)位點(diǎn)的取值均是從A、T、C、G中任取兩個(gè)進(jìn)行組合。例如，假設(shè)X和Y是某個(gè)位點(diǎn)的堿基對，則可能的組合只有XX、XY、YY3個(gè)。若用二進(jìn)制表示，則每個(gè)基因的總長度是12R個(gè)比特，而W個(gè)基因序列的總長度是12RW個(gè)比特。

2.4 基因表達(dá)與變量的對應(yīng)關(guān)系設(shè)計(jì)

模型式（1）中的未知數(shù)變量有n+K個(gè)，其中n=|H|為實(shí)數(shù)變量個(gè)數(shù)，K=|D|為0、1整數(shù)變量個(gè)數(shù)。顯然，取W=n+K。于是，W個(gè)基因被傳染病毒Z攻擊，等價(jià)于n+K個(gè)未知數(shù)變量被處理。特別注意，傳染病毒Z每次只選擇幾個(gè)基因?qū)嵤┕簦艘馕吨看沃恍杼幚韓+K個(gè)未知數(shù)變量中的幾個(gè)變量，于是達(dá)到了天然降維的目的。

W個(gè)致病基因被攻擊后發(fā)生改變，在外顯特征上，動(dòng)物個(gè)體會(huì)表現(xiàn)出一些奇特的特征，如艾滋病人的免疫喪失，麻風(fēng)病人臉上出天花，SARS病人呼吸窘迫等。然而，在SIR-DNA算法中，將動(dòng)物個(gè)體的外顯特征解釋為個(gè)體性狀，而一種個(gè)體性狀對應(yīng)著模型式（1）的一個(gè)可行解。由所有個(gè)體性狀組成的表稱為個(gè)體性狀表，如圖4（a）所示。

Fig.3 Gene expression before and after immunization and individual animals are infected圖3 動(dòng)物個(gè)體被傳染病毒攻擊前后和獲得免疫前后的基因表達(dá)

Fig.4 Corresponding relation between W disease-causing and immune-causing genes and W variables圖4 W個(gè)致病基因和免疫基因與W個(gè)變量的對應(yīng)關(guān)系

從外顯特征看，若一個(gè)動(dòng)物個(gè)體染上傳染病Z，也獲得對傳染病Z的免疫，該個(gè)體屬于S類；若一個(gè)動(dòng)物個(gè)體染上傳染病Z，即其致病基因發(fā)生改變，該個(gè)體屬于I類；若一個(gè)動(dòng)物個(gè)體染上傳染病Z后被治愈，即其致病基因復(fù)原，其免疫基因發(fā)生改變，該個(gè)體屬于R類。

圖4描述了W個(gè)致病基因（如圖4（b）所示）、W個(gè)免疫基因（如圖4（c）所示）與W個(gè)變量的對應(yīng)關(guān)系。圖4（a）中，每個(gè)圓點(diǎn)代表一個(gè)位點(diǎn)；前n個(gè)為基因?qū)?yīng)n個(gè)實(shí)數(shù)型變量；后K個(gè)為基因?qū)?yīng)K個(gè)0-1整數(shù)型變量。每個(gè)基因上選擇一個(gè)圓點(diǎn)，并組成一種個(gè)體性狀，它對應(yīng)著動(dòng)物個(gè)體的一種外顯特征。在本算法中，一種個(gè)體性狀對應(yīng)著一個(gè)可行解。例如，圖中的所有黑點(diǎn)組成一個(gè)可行解。于是，對模型式（1）的求解最終轉(zhuǎn)化為在圖4（a）所示的各種離散點(diǎn)組合中，找出使模型式（1）的目標(biāo)函數(shù)值達(dá)到最小的可行解，此為要尋找的全局最優(yōu)可行解。

當(dāng)傳染病毒Z發(fā)起對動(dòng)物個(gè)體的攻擊時(shí)，因?yàn)樵搨魅静《久看沃贿x擇幾個(gè)基因?qū)嵤┢茐?，例如圖4中傳染病毒Z在一次攻擊中只選擇基因1、基因2和基因n+1實(shí)施破壞，故只需對基因1、基因2和基因n+1所對應(yīng)的變量進(jìn)行處理即可。

由于一個(gè)位點(diǎn)只有12比特寬，要想獲得非常高精度的可行解，一個(gè)基因所包含的位點(diǎn)數(shù)為必需滿足212R≥Mr，即R≥(lbMr)/12。一般情況下，當(dāng)R=1,2,3時(shí)，可以獲得10-4、10-7、10-10級別的精度。

模型式（1）中的變量既包含實(shí)數(shù)變量，又包含整數(shù)變量，其各值的編碼方法如下所述。

（1）實(shí)變量編碼方法

式中，ω為權(quán)重；δj稱為變量偏差，該值是從下列集合中選擇：

（2）整數(shù)變量編碼方法

K個(gè)0-1整數(shù)變量，可以直接采用的編碼方法。

2.5 演化算子設(shè)計(jì)

對圖2進(jìn)行分解，存在下列兩種情況：

情況1在時(shí)期t，個(gè)體從類別A轉(zhuǎn)移到類別C，如圖5（a）所示，其中A,C∈{S,I,R}，但A≠C。大量的個(gè)體類別轉(zhuǎn)移都屬于這種情況。例如S→I、I→R、R→S等。

Fig.5 Two classes transfer of an individual at time t圖5 時(shí)期t個(gè)體的兩種類別轉(zhuǎn)移情形

為了能使某個(gè)體從類別A轉(zhuǎn)移到類別C，將已處于類別C的若干其他個(gè)體的某些基因上的位點(diǎn)傳給該個(gè)體的對應(yīng)位點(diǎn)，也即使該個(gè)體的對應(yīng)位點(diǎn)具有類別C的個(gè)體的基因上的位點(diǎn)。此舉實(shí)現(xiàn)了該個(gè)體從類別A轉(zhuǎn)移到類別C。例如，對于S→I轉(zhuǎn)移，將已處于I類的若干個(gè)體的某些基因上的位點(diǎn)傳給處于S類的某個(gè)體，即可使其感染上傳染病，即實(shí)現(xiàn)了S→I轉(zhuǎn)移。

情況2在時(shí)期t，當(dāng)個(gè)體處于某個(gè)類別A時(shí)，A∈{S,I,R}，并沒有發(fā)生類別轉(zhuǎn)移，即相當(dāng)于A→A，如圖5（b）所示。圖2中的每個(gè)節(jié)點(diǎn)隱含了圖5（b）所示的情形，例如S→S、I→I和R→R。

當(dāng)某個(gè)體處于類別A時(shí)，為了能使該個(gè)體向好的方向發(fā)展，但其類別又保持不變，將已處于同樣類別A，但其IPI指數(shù)要高于該個(gè)體IPI指數(shù)的若干個(gè)強(qiáng)壯個(gè)體的某些基因上的位點(diǎn)傳給該個(gè)體的對應(yīng)基因上的位點(diǎn)，也就是將IPI指數(shù)高的強(qiáng)壯個(gè)體向IPI指數(shù)低的虛弱個(gè)體傳遞強(qiáng)壯基因的信息，使得這些虛弱個(gè)體能向好的方向發(fā)展。

從易感、發(fā)病和治愈的個(gè)體中分別隨機(jī)挑選出L個(gè)個(gè)體，L≥1，u∈{S,I,R}，這些個(gè)體的IPI指數(shù)要高于當(dāng)前個(gè)體i的IPI指數(shù)，分別形成強(qiáng)壯易感者集合、染病者集合和治愈者集合CPu={Xj1(t),Xj2(t),…,XjL(t)}。

（1）S-I算子設(shè)計(jì)。設(shè)當(dāng)前易感個(gè)體i被感染，變?yōu)槿静≌?，其類別從S→I，操作步驟如下：

（1.1）從W個(gè)致病基因中以概率p(p＜E0)隨機(jī)挑選出若干個(gè)致病基因，E0表示個(gè)體的致病基因被選中的最大概率，即基因被病毒攻擊的最大概率。

（1.2）若致病基因j被選中，則：

（1.2.1）從致病基因j的R個(gè)位點(diǎn)中隨機(jī)選出一個(gè)位點(diǎn)k與Y進(jìn)行異或運(yùn)算，記錄Inf(i,1)=j，Inf(i,2)=k。其中，若j≤n，則Y=111111；若n＜j≤n+K，則Y=1。此操作表示傳染病Z攻擊個(gè)體i的DNA。此步驟稱為基因位點(diǎn)操作。

（1.2.2）從強(qiáng)壯染病者集合CPI中隨機(jī)挑選出一個(gè)個(gè)體s，將其對應(yīng)性狀表中的第j列的狀態(tài)值傳給個(gè)體i。如圖6所示，個(gè)體i的第j個(gè)致病基因原來的狀態(tài)值vi轉(zhuǎn)變?yōu)閭€(gè)體s的第j個(gè)致病基因原來的狀態(tài)值vs。該步驟相當(dāng)于染病者傳播病毒信息給易感者。此步驟稱為個(gè)體性狀操作。

Fig.6 State value vitransferring into vsof the j th disease gene of individual i圖6 個(gè)體i的第j個(gè)致病基因狀態(tài)值vi變?yōu)関s

（2）I-R算子設(shè)計(jì)。設(shè)當(dāng)前染病個(gè)體i被治愈，變?yōu)橹斡?，其類別從I→R，操作步驟如下：

（2.1）將個(gè)體i的所有發(fā)生變異的致病基因位點(diǎn)進(jìn)行復(fù)原，即查找個(gè)體i的Inf(i,1)和Inf(i,2)，找到個(gè)體i的發(fā)生變異的致病基因和位點(diǎn)，將其與Y再次進(jìn)行異或運(yùn)算，即可實(shí)現(xiàn)復(fù)原。

（2.2）從W個(gè)免疫基因中以概率p(p＜F0)隨機(jī)挑選出若干個(gè)免疫基因，F(xiàn)0表示個(gè)體的免疫基因被選中的最大概率，F(xiàn)0=E0。若免疫基因j被選中，則：

（2.2.1）從免疫基因j的R個(gè)位點(diǎn)中隨機(jī)選出一個(gè)位點(diǎn)k與Y進(jìn)行異或運(yùn)算，記錄Rec(i,1)=j，Rec(i,2)=k。其中，若j≤n，則Y=111111；若n＜j≤n+K，則Y=1。此操作表示治愈者獲得免疫。

（2.2.2）從強(qiáng)壯染病者集合CPR中隨機(jī)挑選出一個(gè)個(gè)體s，將其對應(yīng)性狀表中的第j列的狀態(tài)值傳給個(gè)體i。如圖6所示。該步驟相當(dāng)于強(qiáng)壯治愈者將其免疫信息傳給染病者，使其治愈并獲得免疫。

（3）R-S算子設(shè)計(jì)。設(shè)當(dāng)前治愈者個(gè)體i失去免疫力，變?yōu)橐赘姓?，其類別從R→S，操作步驟如下：

（3.1）將個(gè)體i的所有使其獲得免疫力的免疫基因位點(diǎn)進(jìn)行復(fù)原，即查找個(gè)體i的Rec(i,1)和Rec(i,2)，找到個(gè)體i的使其獲得免疫力的免疫基因和位點(diǎn)，將這些免疫基因位點(diǎn)與Y再次進(jìn)行異或運(yùn)算，即可實(shí)現(xiàn)復(fù)原，從而使其喪失免疫力。

（3.2）從W個(gè)免疫基因中以概率p(p＜F0)隨機(jī)挑選出若干個(gè)免疫基因，若免疫基因j被選中，則從強(qiáng)壯易感者集合CPS中隨機(jī)挑選出一個(gè)個(gè)體s，將其對應(yīng)性狀表中的第j列的狀態(tài)值傳給個(gè)體i。如圖6所示。該步驟相當(dāng)于強(qiáng)壯易感者將其易感信息傳給治愈者，使其治愈并喪失免疫力。

（4）S-S算子、I-I算子、R-R算子設(shè)計(jì)。因這3個(gè)算子不存在類別轉(zhuǎn)換，故不進(jìn)行基因位點(diǎn)操作，但要進(jìn)行個(gè)體性狀操作。這3個(gè)算子的設(shè)計(jì)方法類似，其設(shè)計(jì)方法為：設(shè)當(dāng)前個(gè)體i的類別為u，u∈{S,I,R}，從個(gè)體i的W個(gè)基因中以概率p(p＜E0)隨機(jī)挑選出若干個(gè)基因，若基因j被選中，則從該類別的強(qiáng)壯者集合CPu中隨機(jī)挑選出一個(gè)個(gè)體s，將其對應(yīng)個(gè)體性狀表中的第j列的性狀值傳給個(gè)體i。如圖6所示。該步驟相當(dāng)于強(qiáng)壯者將其特有的信息傳給當(dāng)前個(gè)體i，使該個(gè)體變?yōu)楦鼜?qiáng)狀態(tài)。此操作相當(dāng)于弱者向強(qiáng)者學(xué)習(xí)，使自身變強(qiáng)壯。W個(gè)基因選擇方法是：若u=S，則選W個(gè)未染病的致病基因；若u=I，則選W個(gè)已染病的致病基因；若u=R，則選W個(gè)已免疫的免疫基因。

（5）生長算子設(shè)計(jì)。個(gè)體i的生長算子可以描述為：

式中，Vi(t)是Xi(t)的一個(gè)副本，用于臨時(shí)保存Xi(t)的值；Xi(t)=(xi1(t),xi2(t),…,xin(t))；函數(shù)IPI(Vi(t))和IPI(Xi(t))按式（4）計(jì)算。

2.6 SIR-DNA算法構(gòu)造方法

算法1 SIR-DNA算法

（1）初始化：①令時(shí)期t=0，最大迭代次數(shù)G=108，N=100～200，精度ε=10-8，R=3～5，L=3～ 5，F(xiàn)0=E0=0.000 1～0.01；②初始化N個(gè)正常個(gè)體X1(0),X2(0),…,XN(0)；③在N個(gè)正常個(gè)體中隨機(jī)選擇30%的個(gè)體使其致??；④在剩下的正常個(gè)體中隨機(jī)選擇30%的個(gè)體使其獲得免疫力。

（2）計(jì)算ai=Rand(0,1)，bi=Rand(0,1)，ci=Rand(0,1)；Si(0)=ai/(ai+bi+ci)，Ii(0)=bi/(ai+bi+ci)，Ri(0)=1-Si(0)-Ii(0),i=1,2,…,N。

（3）計(jì)算個(gè)體i的SIR類別：SIRi(0)=GetSIR(Si(0),Ii(0),Ri(0))，i=1,2,…,N；函數(shù)GetSIR()用于確定個(gè)體i將處于何種類別。

（4）執(zhí)行下列操作：

（5）結(jié)束。

函數(shù)GetSIR(pS,pI,pR)的定義如下：

2.7 SIR-DNA算法的降維特征的理論分析

將式（1）的目標(biāo)函數(shù)展開成如下形式：

SIR-DNA算法的求解過程最終歸結(jié)為不斷修改變量x1,x2,…,xn的取值，使得目標(biāo)函數(shù)值A(chǔ)不斷向理論最優(yōu)解靠近。在演化過程中，n個(gè)變量中的每個(gè)變量取值的修改過程是相互獨(dú)立的；當(dāng)一個(gè)變量的取值被修改后，只能使目標(biāo)函數(shù)值A(chǔ)向3個(gè)方向變化：比當(dāng)前值更好，比當(dāng)前值更差，與當(dāng)前值一樣。不妨假定這3種情形出現(xiàn)的概率相同，均為1/3。假設(shè)SIR-DNA算法每次進(jìn)化有m個(gè)變量的取值同時(shí)改變，則能使目標(biāo)函數(shù)值A(chǔ)向好的方向改變的概率為p=1/3m。

顯然，當(dāng)m=n時(shí)，即每次進(jìn)化時(shí)n個(gè)變量的取值同時(shí)改變，此時(shí)能使目標(biāo)函數(shù)值A(chǔ)向好的方向發(fā)展的概率最低，即p=3-n。因此，每次進(jìn)化時(shí)同時(shí)對n個(gè)變量的取值進(jìn)行更新是最不可取的策略，該策略既耗時(shí)，又使目標(biāo)函數(shù)值A(chǔ)向好的方向發(fā)展的概率達(dá)到最低。

另一方面，當(dāng)m=1時(shí)，即每次進(jìn)化時(shí)只修改1個(gè)變量的取值，此時(shí)能使目標(biāo)函數(shù)值A(chǔ)向好的方向發(fā)展的概率最大，即p=1/3。因此，每次進(jìn)化時(shí)只對1個(gè)變量的取值進(jìn)行更新是可取的策略，該策略能使A向好的方向發(fā)展的概率達(dá)到最大。

上面的分析給出了SIR-DNA算法在求解過程中每次處理變量數(shù)的一般規(guī)律，這是一個(gè)普適規(guī)律。然而，由于求解過程的復(fù)雜性，每次只修改1個(gè)變量的取值不一定性能最佳，但一般規(guī)律是：無論n多大，每次參與運(yùn)算的變量數(shù)都不要超過10個(gè)，即1≤nE0＜10。

上述分析表明，SIR-DNA算法具有天然降維特征。

2.8 SIR-DNA算法時(shí)間復(fù)雜度

SIR-DNA算法的時(shí)間復(fù)雜度如表3所示。

Table 3 Time complexity of SIR-DNA表3 SIR-DNA算法的時(shí)間復(fù)雜度計(jì)算表

2.9 SIR-DNA算法收斂性和穩(wěn)定性分析

從S-S、S-I、I-I、I-R、R-R、R-S算子的定義可知，任何一個(gè)試探解的新一代的生成只與該試探解的當(dāng)前狀態(tài)有關(guān)，而與該試探解以前是如何演變到當(dāng)前狀態(tài)的歷程無關(guān)，因而SIR-DNA算法的演化過程具有Markov特性。因老的試探解無需保留，故可以使算法的空間復(fù)雜度降到最低。

從生長算子的定義可知，若從當(dāng)前位置出發(fā)，下一步搜索方向只有兩個(gè)，即要么向比當(dāng)前位置更好的方向搜索，要么留在當(dāng)前位置不動(dòng)。因此，SIRDNA算法具有“步步不差”的搜索特征。

由于SIR-DNA算法的演化過程具有Markov特性和“步步不差”特性，根據(jù)文獻(xiàn)[13]可得如下定理。

定理1 SIR-DNA算法具有全局收斂性。

SIR-DNA算法的穩(wěn)定性依賴于SIR傳染病動(dòng)力學(xué)模型的解的穩(wěn)定性，即微分方程組（2）的解的穩(wěn)定性。文獻(xiàn)[12]給出了該模型的解保持全局穩(wěn)定性的條件，本文不再贅述。于是，SIR傳染病模型可以幫助SIR-DNA算法選擇最合理的參數(shù)實(shí)現(xiàn)穩(wěn)定收斂。

3 SIR-DNA算法與其他算法比較

本文使用CEC2013[14]所提供的國際上通用的基準(zhǔn)函數(shù)來測試SIR-DNA算法的性能，本文選擇的6個(gè)基準(zhǔn)函數(shù)，如表4所示。在表4中，O是一個(gè)n維決策向量。這里用SIR-DNA算法去求解表4所示的基準(zhǔn)函數(shù)，其參數(shù)是n=50，N=200，O的值隨機(jī)產(chǎn)生。選擇7種優(yōu)化算法與SIR-DNA算法進(jìn)行比較，這些算法如表5所示。

用這些算法獨(dú)立求解每個(gè)基準(zhǔn)函數(shù)51次，計(jì)算結(jié)果如表6所示。表6的排名1是按平均最佳目標(biāo)函數(shù)值排序，排名2是按平均最佳目標(biāo)函數(shù)值和CPU時(shí)間排序。表7是8種算法所得最優(yōu)解的非參數(shù)Wilcoxon秩和檢驗(yàn)。表中，h-val=1表明SIR-DNA算法能夠以99%的概率優(yōu)于其他算法，h-val=-1表示SIR-DNA算法明顯劣于其他算法，而h-val=0表示SIR-DNA算法與其他算法的結(jié)果差異不顯著。

Table 4 Benchmark functions表4 基準(zhǔn)函數(shù)

從表6中可以看出，這8個(gè)算法按精度排序如下：

按精度和CPU耗時(shí)的排序如下：

從表7中可以知道，SIR-DNA算法的性能明顯優(yōu)于7種被比較的算法。

圖7（a）～（c）說明8個(gè)算法求解F2、F20和F22時(shí)的樣本收斂曲線，其中的水平和垂直軸采用對數(shù)刻度。從表6中可以看出，當(dāng)SIR-DNA算法求解F2、F20時(shí)，其發(fā)現(xiàn)最優(yōu)解所需時(shí)間要比被比較算法少；當(dāng)SIR-DNA算法求解F22時(shí)，其發(fā)現(xiàn)最優(yōu)解的精度和所需時(shí)間稍遜于DSDA算法，但比其他算法要好。綜合看來，SIR-DNA算法的綜合性能要優(yōu)于被比較算法，表明其計(jì)算速度快。從圖7（a）～（b）可以看出，SIR-DNA算法的收斂曲線大部分在被比較算法的左側(cè)，表明其收斂速度很快；從圖7（c）可以看出，SIRDNA算法的收斂曲線在被比較算法的右側(cè)，表明其收斂速度慢于某些算法，但是只有SIR-DNA算法獲得的最優(yōu)解最好。

Table 5 Parameters of 7 compared algorithms表5 7種參與比較的算法的參數(shù)

Table 7 Results comparison of Wilcoxon rank sum test(α=0.01)表7 Wilcoxon秩和檢驗(yàn)結(jié)果比較(α=0.01)

4 結(jié)論

本文采用SIR傳染病動(dòng)力學(xué)理論與DNA分子結(jié)構(gòu)理論相結(jié)合的方法構(gòu)建出SIR-DNA算法。在該算法中，SIR傳染病模型的優(yōu)勢和特性如下：

（1）S-S、S-I、I-I、I-R、R-R、R-S算子不與優(yōu)化問題相關(guān)，從而使SIR-DNA算法具有通用性。

（2）利用傳染病傳播時(shí)所誘發(fā)的6個(gè)狀態(tài)轉(zhuǎn)移，從多種角度天然地實(shí)現(xiàn)個(gè)體之間信息的充分交換，降低了個(gè)體陷入局部最優(yōu)的概率。

（3）因病毒每次攻擊的是個(gè)體的很少部分特征，每次運(yùn)算只涉及到很少一部分特征參與運(yùn)算，從而實(shí)現(xiàn)天然降維。

（4）SIR-DNA算法利用SIR傳染病動(dòng)力學(xué)理論，使搜索過程達(dá)到生態(tài)穩(wěn)定時(shí)出現(xiàn)收斂，而SIR傳染病模型可以幫助算法選擇合理的參數(shù)實(shí)現(xiàn)收斂。

在SIR-DNA算法中，由于采用SIR傳染病模型作為基礎(chǔ)，從一個(gè)新的角度構(gòu)造了能反映DNA特征的優(yōu)化算法，其特征在于：

（1）SIR-DNA算法不是直接與DNA算法進(jìn)行結(jié)合，從而避免了DNA算法的缺陷，即DNA計(jì)算時(shí)需要借助生物實(shí)驗(yàn)，無法求解高維優(yōu)化問題等。

Fig.7 Convergence curves when 8 algorithms solve F2、F20、F22圖7 8個(gè)算法求解F2、F20、F22時(shí)的收斂曲線

（2）SIR-DNA算法納入DNA特征的方法是，將DNA中的基因和位點(diǎn)進(jìn)行編碼，采用隨機(jī)修改關(guān)鍵基因及其位點(diǎn)的方法來模擬病毒對DNA的攻擊或獲得免疫力，從而無需生物實(shí)驗(yàn)即可進(jìn)行計(jì)算。

迄今為止，已經(jīng)發(fā)現(xiàn)了數(shù)以千計(jì)的傳染病。本文的SIR-DNA算法為將這些傳染病轉(zhuǎn)化為能求解的非常復(fù)雜的優(yōu)化問題的群智能優(yōu)化算法提供了參考。

[1]Adleman L M.Molecular computation of solution to combinatorial problems[J].Science,1994,266(5187):1021-1024.

[2]Zhang Xuncai,Zhao Hailan,Cui Guangzhao,et al.Research advances and prospect of DNAcomputing[J].Computer Engineering andApplications,2007,43(10):44-47.

[3]Takahashi K,Yaegashi S,Asanuma H,et al.Photo-and thermoregulation of DNA nanomachines[C]//LNCS 3892:Proceedings of the 11th International Workshop on DNA Computing,London,Jun 6-9,2005.Berlin,Heidelberg:Springer,2005:336-346.

[4]Wu Fan,Li Kenli.An algorithm in tile assembly model forNqueen problem[J].Acta Electronica Sinica,2013,41(11):2174-2180.

[5]Li Kenli,Luo Xing,Wu Fan,et al.An algorithm in tile assembly model for maximum clique problem[J].Journal of Computer Research and Development,2013,50(3):666-675.

[6]Yao Qing'an,Zheng Hong,Wang Hongmei.DNA computing model for shortest path problem based onk-armed molecule and sticker operation[J].Journal of Jilin University:Information Science Edition,2014,32(6):653-656.

[7]Xu Jinglei,Zhao Hongchao,Liu Xiyu.Closed circle DNA algorithm of traveling salesman problem[J].Computer Engineering&Science,2014,36(1):111-114.

[8]Ma Ying,Yin Zhixiang.DNA computing model for the graph vertex coloring problem by plasmids[J].Journal of Anhui University of Science and Technology:Natural Science,2015,35(2):64-67.

[9]Wu Xue,Song Chenyang,Zhang Nan,et al.DNA algorithm for maximum matching problem based on sticker computation model[J].Computer Science,2013,40(12):127-132.

[10]Kermack W O,Mckendrick A G.Contributions to the mathematical theory of epidemics[J].Proceedings of the Royal Society of London:SeriesA,1927,115(772):700-721.

[11]Kermack W O,Mckendrick A G.Contributions to the mathematical theory of epidemics(II):the problem of endemicity[J].Proceedings of the Royal Society of London,1932,138(795):55-83.

[12]Wang Qian,Wang Tingting.Stability analysis of a SIRS epidemic model[J].Journal of Capital Normal University:Natural Science Edition,2016,37(2):5-11.

[13]Huang Guangqiu.SIS epidemic model-based optimization[J].Journal of Computation Science,2014,5(1):32-50.

[14]Liang J J,Qu B Y,Suganthan P N,et al.Problem definitions and evaluation criteria for the CEC 2013 special session on real-parameter optimization,201212[R].Singapore:Nanyang Technological University,2013.

[15]Chuang Yaochen,Chen C T,Hwang C.A simple and efficient real-coded genetic algorithm for constrained optimization[J].Applied Soft Computing,2016,38:87-105.

[16]Koro?ec P,?ilc J,Filipi? B.The differential ant-stigmergy algorithm[J].Information Sciences,2012,192:82-97.

[17]Beheshti Z,Shamsuddin S M.Non-parametric particle swarm optimization for global optimization[J].Applied Soft Computing,2015,28:345-359.

[18]Al-Roomi A R,El-Hawary M E.Metropolis biogeographybased optimization[J].Information Sciences,2016,360:73-95.[19]Mukherjee R,Debchoudhury S,Das S.Modified differential evolution with locality induced genetic operators for dynamic optimization[J].European Journal of Operational Research,2016,253(2):337-355.

[20]Zhao Zhiwei,Yang Jingming,Hu Ziyu,et al.A differential evolution algorithm with self-adaptive strategy and control parameters based on symmetric Latin hypercube design for unconstrained optimization problems[J].European Journal of Operational Research,2016,250(1):30-45.

[21]Li Genghui,Cui Laizhong,Fu Xianghua,et al.Artificial bee colony algorithm with gene recombination for numerical function optimization[J].Applied Soft Computing,2017,52:146-159.

附中文參考文獻(xiàn):

[2]張勛才,趙海蘭,崔光照.DNA計(jì)算的研究進(jìn)展及展望[J].計(jì)算機(jī)工程與應(yīng)用,2007,43(10):44-47.

[4]吳帆,李肯立.基于自組裝的N皇后問題DNA計(jì)算算法[J].電子學(xué)報(bào),2013,41(11):2174-2180.

[5]李肯立,羅興,吳帆,等.基于自組裝模型的最大團(tuán)問題DNA計(jì)算算法[J].計(jì)算機(jī)研究與發(fā)展,2013,50(3):666-675.

[6]姚慶安,鄭虹,王紅梅.基于k-臂分子求解最短路徑的DNA計(jì)算模型[J].吉林大學(xué)學(xué)報(bào):信息科學(xué)版,2014,32(6):653-656.

[7]徐京雷,趙洪超,劉希玉.旅行商問題的閉環(huán)DNA算法[J].計(jì)算機(jī)工程與科學(xué),2014,36(1):111-114.

[8]馬瑩,殷志祥.圖頂點(diǎn)著色問題的質(zhì)粒DNA計(jì)算[J].安徽理工大學(xué)學(xué)報(bào):自然科學(xué)版,2015,35(2):64-67.

[9]吳雪,宋晨陽,張楠,等.最大匹配問題的粘貼DNA算法[J].計(jì)算機(jī)科學(xué),2013,40(12):127-132.

[12]王茜,王婷婷.SIRS傳染病模型的穩(wěn)定性分析[J].首都師范大學(xué)學(xué)報(bào):自然科學(xué)版,2016,37(2):5-11.