高飛 胡道楠 童恒慶 王傳美

(武漢理工大學(xué)理學(xué)院,武漢 430070)

(2018年2月2日收到;2018年4月16日收到修改稿)

1 引 言

Willis環(huán)腦動(dòng)脈血管瘤是一種高致死率疾病[1],治療過(guò)程中有可能會(huì)因基于植入血流導(dǎo)向裝置[2?4]和支架聯(lián)合彈簧圈[5]等治療方式引起不明原因的遲發(fā)性動(dòng)脈瘤破裂[6].動(dòng)脈瘤的延遲性破裂危害巨大,一旦出現(xiàn)將嚴(yán)重危及病人生命.

生理上動(dòng)脈瘤破裂與否與血流速度密切相關(guān),故以血流速度為主要研究對(duì)象,以血流動(dòng)力學(xué)相關(guān)原理建立的腦動(dòng)脈血管瘤動(dòng)力學(xué)模型系統(tǒng)[7?10]在臨床及理論研究中發(fā)揮了重要作用.動(dòng)脈瘤破裂表現(xiàn)為血流速度的巨大改變,亦即出現(xiàn)“尖峰”,雖然不是所有尖峰都會(huì)引起動(dòng)脈瘤破裂,但是沒有尖峰意味著血流流速穩(wěn)定,病情穩(wěn)定,系統(tǒng)表現(xiàn)為穩(wěn)定狀態(tài);反之,存在尖峰的系統(tǒng)則表示為混沌狀態(tài).

最初的Willis環(huán)腦動(dòng)脈血管瘤系統(tǒng)(Willis aneurysm system,WAS)由Austin[11]利用實(shí)驗(yàn)?zāi)M得到,現(xiàn)有研究多以整數(shù)階阻尼項(xiàng)腦動(dòng)脈瘤系統(tǒng)[12]為基本模板,對(duì)其進(jìn)行混沌理論分析[13?15]和模型改進(jìn)[16,17].在模型改進(jìn)方面,近年來(lái)主要有基于藥物的整數(shù)階WAS模型[16]和基于血液軟物質(zhì)性[18]、黏彈性[19]特性的分?jǐn)?shù)階WAS[17]等.但是,對(duì)于不明原因引起的遲發(fā)性動(dòng)脈瘤破裂,也就是系統(tǒng)中尖峰延遲出現(xiàn)的情況(即“時(shí)滯”),上述模型并不能給出合理的描述和解釋.

近年來(lái),分?jǐn)?shù)階微積分作為一種有用的數(shù)學(xué)工具被廣泛應(yīng)用于生物及醫(yī)學(xué)方面[20,21];而時(shí)滯一直存在于現(xiàn)實(shí)中并影響著系統(tǒng)的動(dòng)態(tài),故分?jǐn)?shù)階時(shí)滯系統(tǒng)引起了學(xué)者們的廣泛研究和關(guān)注[22,23].時(shí)滯加入后,會(huì)破壞原來(lái)系統(tǒng)的穩(wěn)定性并影響著系統(tǒng)動(dòng)力學(xué)行為,使系統(tǒng)變得更加復(fù)雜,而臨床上腦動(dòng)脈瘤里血管內(nèi)情況亦是錯(cuò)綜復(fù)雜.現(xiàn)有WAS的相關(guān)機(jī)理及理論滯后于臨床現(xiàn)實(shí),因此,研究帶有時(shí)滯的相關(guān)模型,將在一定程度上為腦動(dòng)脈血管瘤的臨床診斷給出理論指導(dǎo).

鑒于此,本文構(gòu)造了分?jǐn)?shù)階Willis環(huán)腦遲發(fā)性動(dòng)脈瘤時(shí)滯系統(tǒng)(fractional Willis aneurysm system with time-delay,FWASTD)并對(duì)其進(jìn)行了數(shù)值仿真和理論分析:通過(guò)與非時(shí)滯分?jǐn)?shù)階Willis環(huán)腦動(dòng)脈血管瘤系統(tǒng)(fractional Willis aneurysm system,FWAS)做對(duì)比,驗(yàn)證了其時(shí)滯有效性;用傳統(tǒng)動(dòng)力學(xué)方法驗(yàn)證系統(tǒng)混沌,探究了時(shí)滯給系統(tǒng)帶來(lái)的豐富動(dòng)力學(xué)行為;利用分?jǐn)?shù)階時(shí)滯穩(wěn)定性原理實(shí)現(xiàn)了FWASTD的混沌控制和自同步混沌控制.本文為腦動(dòng)脈瘤系統(tǒng)研究和臨床診斷提供了相應(yīng)的理論基礎(chǔ)和相關(guān)參考.

2 理論知識(shí)

定義1[24]設(shè)α是一個(gè)正實(shí)數(shù),令n?1 6 α

定理1[25]對(duì)于帶初值的時(shí)滯系統(tǒng)

定理2[26]對(duì)于分?jǐn)?shù)階時(shí)滯非線性系統(tǒng),Dαx(t)=f(x(t),x(t?τ)),當(dāng)分?jǐn)?shù)階微分階次0< α<1,f(x(t),x(t?τ))滿足Lipschiz條件時(shí),若存在正定矩陣P和半正定矩陣Q,對(duì)于任意的狀態(tài)變量分?jǐn)?shù)階時(shí)滯非線性系統(tǒng)仍然滿足

則分?jǐn)?shù)階時(shí)滯非線性系統(tǒng)是Lyapunov穩(wěn)定的.

3 FWASTD及其性質(zhì)

3.1 建立FWASTD

從文獻(xiàn)[17]可知,FWAS如下:

其中,x和y分別為血流變化率和變化率的加速度;F和μ作為重要的生理參量,分別代表脈沖壓和血流阻力;ω,α,β,γ是涉及到心率和血管的生理參量指數(shù),在病理上與動(dòng)脈瘤狀況息息相關(guān)[14];q1和q2作為分?jǐn)?shù)階次可以精細(xì)刻畫腦動(dòng)脈瘤系統(tǒng).

由于臨床中血流速度呈現(xiàn)上下波動(dòng)、不斷改變的狀態(tài),說(shuō)明血流加速度存在于血流運(yùn)動(dòng)變化中并起著推動(dòng)血流速度改變的重要作用.一旦血流有顯著變化,最開始的病變一定是從血流加速度開始,故引言中提到的尖峰延遲的根源應(yīng)在血流加速度上,而表現(xiàn)在血流速度上.取時(shí)滯因子為τ,將其加入到血流加速度中,構(gòu)造時(shí)滯系統(tǒng)FWASTD如下:

3.2 FWASTD的有效性

本節(jié)用FWAS(系統(tǒng)(1))與FWASTD(系統(tǒng)(2))做對(duì)比來(lái)說(shuō)明FWASTD的時(shí)滯有效性.根據(jù)參考文獻(xiàn)[17],取FWAS在混沌狀態(tài)下的系數(shù),即:α=0.9,β=3,γ=2,F=0.1,μ=0.1,ω=1,q1=1,q2=0.95,對(duì)于時(shí)滯τ,嘗試性地取為1,設(shè)FWAS與FWASTD的血流速度分別為x1和x2,利用Matlab 2015b進(jìn)行仿真,用修正的Adams-Bashforth-Moulton方法[27,28]分別求解FWAS與FWASTD的血流速度,得到兩者血流速度的時(shí)間序列圖如圖1所示.

圖1 血流速度x1和x2的時(shí)間序列圖Fig.1.Time course of blood speed x1and x2.

從圖1可以看出,x1和x2在區(qū)間[700,2000]內(nèi)幾乎重合,但從第2000個(gè)時(shí)間點(diǎn)左右開始,相比于x1,x2的尖峰值出現(xiàn)延后,延后的時(shí)間點(diǎn)大約為800個(gè)時(shí)間點(diǎn),從該時(shí)間點(diǎn)后,x1和x2僅有部分重合.從引言部分可知,尖峰延后意味著動(dòng)脈瘤有延遲性破裂的可能.由此可見,FWASTD可以描述不明原因引起的遲發(fā)性動(dòng)脈瘤破裂,也證明了FWASTD的時(shí)滯有效性.

3.3 FWASTD的取值和仿真

3.3.1 FWASTD的分?jǐn)?shù)階次取值

FWASTD的分?jǐn)?shù)階次是整數(shù)階次的推廣.分?jǐn)?shù)階算子本身具有記憶性而優(yōu)于整數(shù)階算子,符合刻畫系統(tǒng)生理病情的需要.為了觀察分?jǐn)?shù)階次取值對(duì)于系統(tǒng)狀態(tài)的影響,以期得到最佳刻畫混沌狀態(tài)的取值,現(xiàn)分別取兩組組合,第一組組合中固定q1取值分別為0.8,0.975和1,得到關(guān)于q2的分岔圖和最大Lyapunov指數(shù)圖;第二組組合中固定q2取值分別為0.75,0.95和1.15,得到關(guān)于q1的分岔圖和最大Lyapunov指數(shù)圖.具體結(jié)果如圖2所示,其中,圖2(a)、圖2(c)、圖2(e)對(duì)應(yīng)第一組組合,圖2(b)、圖2(d)、圖2(f)對(duì)應(yīng)第二組組合,其他參數(shù)參照3.2節(jié)進(jìn)行取值.

從圖2整體來(lái)看,在q1和q2的所有組合里,時(shí)滯系統(tǒng)FWASTD均通過(guò)倍周期分岔道路通往混沌,分岔圖形基本一致,只是取值范圍不同,故可從圖2得出描述混沌的最佳取值區(qū)間,從而解決分?jǐn)?shù)階次取值問(wèn)題.圖2(c)和圖2(d)取值區(qū)間類似,其最大值均逼近1,在分岔圖中剛好互為固定參數(shù),本文選其組合(q1=0.975,q2=0.95)為下文所用.

3.3.2 FWASTD的時(shí)滯取值

為了研究時(shí)滯對(duì)于系統(tǒng)(2)的影響,給出關(guān)于時(shí)滯的分岔圖,目的是研究時(shí)滯可能的取值范圍等問(wèn)題,為后續(xù)研究做鋪墊.取α=0.9,β=3,γ=2,F=0.1,μ=0.1,ω=1,q1=0.975,q2=0.95,變動(dòng)時(shí)滯參數(shù)τ,可以得到hτ(h=0.05為步長(zhǎng))的分岔圖,如圖3所示.

從圖3(a)可知,隨著時(shí)滯的增加,FWASTD系統(tǒng)(2)從混沌到穩(wěn)定.為了更細(xì)致地研究時(shí)滯的區(qū)間段,圖3(a)里混沌部分被放大得到圖3(b),從圖3(b)可以看出,混沌主要集中在hτ=0.05和hτ=0.1這兩個(gè)點(diǎn)上,其他地方并未出現(xiàn)散點(diǎn),這說(shuō)明FWASTD仿真的時(shí)滯點(diǎn)是離散點(diǎn),由此可以確定混沌狀態(tài)下的時(shí)滯τ的具體數(shù)值分別為τ=1和τ=2,即這兩個(gè)值是產(chǎn)生系統(tǒng)混沌的關(guān)鍵.之后進(jìn)行的關(guān)于FWASTD仿真的取值點(diǎn)將從這兩個(gè)點(diǎn)中選取.

3.3.3 FWASTD的仿真

FWASTD作為整數(shù)階WAS的推廣,亦能刻畫遲發(fā)性腦動(dòng)脈瘤系統(tǒng)的混沌狀態(tài),下面針對(duì)FWASTD(系統(tǒng)(2))進(jìn)行仿真驗(yàn)證.取τ=2,α=0.9,β=3,γ=2,F=0.1,μ=0.1,ω=1,q1=0.975,q2=0.95作為參數(shù),可以把FWASTD(系統(tǒng) (2))寫為

圖4為系統(tǒng)(3)在時(shí)滯狀態(tài)的時(shí)間序列圖、相圖和Poincaré截面.

圖2 系統(tǒng)在給定初值下,q1=(a)0.8,(c)0.975,(e)1時(shí)q2的分岔圖和最大Lyapunov指數(shù)圖以及q2=(b)0.75,(d)0.95,(f)1.15時(shí)q1的分岔圖和最大Lyapunov指數(shù)圖Fig.2.With a given initial value,bifurcation and largest Lyapunov exponent(LLE)diagram of system versus q2 when q1=(a)0.8,(c)0.975,(e)1;bifurcation and largest Lyapunov exponent diagram of system versus q1when q2=(b)0.75,(d)0.95,(f)1.15.

圖3 時(shí)滯影響下系統(tǒng)的分岔圖 (a)hτ∈[0,3];(b)hτ∈[0.03,0.12]Fig.3.Bifurcation diagram of the system with time-delay:(a)hτ∈[0,3];(b)hτ∈[0.03,0.12].

從圖4(a)和圖4(b)可以看出,血流速度紊亂,頻頻出現(xiàn)尖峰,說(shuō)明系統(tǒng)呈現(xiàn)出混沌狀態(tài);從圖4(c)可以看出,其軌跡無(wú)規(guī)律;再結(jié)合圖4(d)中系統(tǒng)(3)的Poincaré截面里具有層次結(jié)構(gòu)且成片密集的點(diǎn),亦證明了系統(tǒng)(3)處于混沌狀態(tài),說(shuō)明階次為分?jǐn)?shù)的FWASTD也可刻畫系統(tǒng)的混沌狀態(tài).

圖4 系統(tǒng)在給定初值下,(a)x-t,y-t的時(shí)間歷程圖,(b)x-t的時(shí)間歷程圖,(c)相圖,(d)Poincaré截面Fig.4.The system with a given initial value:(a)Time course of x-t,y-t;(b)time course of x-t;(c)phase diagram;(d)Poincaré section.

3.4 FWASTD中脈沖壓和血流阻力系數(shù)對(duì)系統(tǒng)的影響

在系統(tǒng)(2)中分別把脈沖壓F和血流阻力系數(shù)μ作為變量,其他參數(shù)參照系統(tǒng)(3)保留不變,得到關(guān)于脈沖壓和血流阻力系數(shù)的分岔圖和最大Lyapunov指數(shù)圖,如圖5所示.

多數(shù)文獻(xiàn)[14,16,17]把研究的重點(diǎn)放在生理參量脈沖壓F上,控制也是從脈沖壓這一生理參量入手,但從圖5(a)分岔圖結(jié)合最大Lyapunov指數(shù)可以看出,其最大Lyapunov指數(shù)一直處于0以上,系統(tǒng)持續(xù)處于混沌狀態(tài).

圖5 (a)系統(tǒng)隨脈沖壓F變化的分岔圖和最大Lyapunov指數(shù)圖;(b)系統(tǒng)隨血流阻力系數(shù)μ變化的分岔圖和最大Lyapunov指數(shù)圖Fig.5.(a)Bifurcation and largest Lyapunov exponent diagram of versus pulse pressure F;(b)bifurcation and largest Lyapunov exponent diagram of versus coefficient of blood f l ow dampingμ.

從圖5(b)分岔圖可以看出,隨著血流阻力μ的增大,系統(tǒng)從開始處于混沌狀態(tài)逐漸變?yōu)榉植?直至穩(wěn)定.結(jié)合μ的最大Lyapunov指數(shù)圖可證實(shí)系統(tǒng)會(huì)隨著血流阻力的增加而穩(wěn)定,臨床上亦有促進(jìn)血栓形成來(lái)輔助治療的記載[14],以上種種均說(shuō)明時(shí)滯狀態(tài)下研究血流阻力系數(shù)對(duì)于臨床診斷具有重大意義.

4 FWASTD的控制

由于FWASTD的混沌狀態(tài)表現(xiàn)為遲發(fā)性腦動(dòng)脈瘤破裂,第一要義是避免動(dòng)脈瘤破裂,控制混沌.故本節(jié)將根據(jù)分?jǐn)?shù)階時(shí)滯系統(tǒng)穩(wěn)定性定理對(duì)FWASTD進(jìn)行控制并實(shí)現(xiàn)同步控制.

4.1 時(shí)滯系統(tǒng)解的惟一性分析

下面將證明具有初值條件的FWASTD存在惟一解且滿足Lipschitz條件.

定理3 構(gòu)造具有初值條件的系統(tǒng)如下:

其中,

則含有初值條件的FWASTD(4)存在惟一解.

證 明 取|·|和∥·∥分 別 為 向量范 數(shù) 和 矩陣范數(shù). 令F2(X(t)).對(duì)于?δ>0,區(qū)間連續(xù)并有界.取有

則因區(qū)間[X0?δ,X0+δ]連續(xù)有界,α,β,γ是有范圍的生理參數(shù),因此L1是個(gè)常數(shù),于是對(duì)于有

故G(t,X(t))除t外所有變?cè)獫M足Lipschitz條件,根據(jù)已知條件G(t,X(t))在給定初值的鄰域內(nèi)連續(xù),所以G(t,X(t))滿足定理1,即具有初值條件的FWASTD存在惟一解.定理3證畢.

4.2 FWASTD的混沌控制

取和系統(tǒng)(3)相同的參數(shù),得到有初值條件的FWASTD如下:

從圖4的仿真結(jié)果可以看出相圖軌線無(wú)規(guī)律,Poincaré截面具有一定層次結(jié)構(gòu)和形狀,因此具有初值條件的系統(tǒng)(5)不處于穩(wěn)定狀態(tài).

對(duì)系統(tǒng)(5)設(shè)計(jì)線性控制器:

定理4 當(dāng)k1>0.5,0 6 k26 0.9,k3>1.95時(shí),系統(tǒng)(5)穩(wěn)定.

根據(jù)定理2,系統(tǒng)(6)是Lyapunov穩(wěn)定的.

取k1=0.5,k2=0.9,k3=1.95,系統(tǒng)(6)仿真結(jié)果如圖6所示.

從圖6(a)和圖6(b)可以看出,x,y隨著時(shí)間增長(zhǎng)在[?0.2,0.1]波動(dòng),仿真結(jié)果證明,該控制器可以把FWASTD控制在一個(gè)比較穩(wěn)定的范圍之內(nèi),雖然有周期性波動(dòng),但是波動(dòng)振幅小(上下不超過(guò)0.1),滿足混沌控制需要.圖6(c)和圖6(d)也說(shuō)明了具有線性控制器的FWASTD處于比較穩(wěn)定的狀態(tài),體現(xiàn)在實(shí)際中表現(xiàn)為血管中血流速度和加速度基本保持穩(wěn)定,病情處于可控狀態(tài).近年來(lái),在血管內(nèi)使用血流導(dǎo)向裝置成為治療腦動(dòng)脈瘤的一種主流辦法.血流導(dǎo)向裝置可以通過(guò)改變血流速度和加速度等血流動(dòng)力學(xué)因子來(lái)改善動(dòng)脈瘤內(nèi)部的血流情況從而達(dá)到治療疾病的目的,但現(xiàn)有的血流導(dǎo)向裝置在設(shè)計(jì)和材料使用上仍有缺陷[29],本文為其提供了相關(guān)的設(shè)計(jì)參考.

4.3 FWASTD的同步控制

下面考慮FWASTD同步控制,把FWAS作為驅(qū)動(dòng)系統(tǒng):

圖6 (a)x-t曲線,y-t曲線;(b)時(shí)間歷程;(c)相圖;(d)Poincaré截面Fig.6.(a)Curves of x-t and y-t;(b)time course;(c)phase diagram;(d)Poincaré section.

把帶控制器的FWASTD作為響應(yīng)系統(tǒng):

誤差系統(tǒng)為

數(shù)值仿真驗(yàn)證時(shí),取τ=2,α=0.9,β=3,γ=2,F=0.1,ω=1,q1=0.975,q2=0.95,μ=0.1,k1=?0.5,k2=?2.1,k3=1.1,t表示時(shí)間,誤差系統(tǒng)仿真結(jié)果如圖7所示.

圖7 (a)e1隨時(shí)間的變化;(b)e2隨時(shí)間的變化Fig.7.(a)Curve of e1with time change;(b)curve of e2with time change.

從圖7(a)和圖7(b)可以看出,其誤差逐漸趨于0,由此可見,FWASTD和FWAS可以實(shí)現(xiàn)同步.實(shí)際治療中,由于許多生理參數(shù)無(wú)法確定等原因,時(shí)滯系統(tǒng)FWASTD無(wú)法得到像4.2節(jié)那樣精準(zhǔn)的混沌控制,然而從自同步控制過(guò)程可以看出,其控制參數(shù)主要依賴于血流阻力參數(shù),控制變量少且易得,而FWAS可以通過(guò)脈沖藥物等進(jìn)行控制,操作簡(jiǎn)單可控,更符合實(shí)際.

5 結(jié) 論

鑒于治療過(guò)程中有可能會(huì)因治療方式等不明原因引起遲發(fā)性動(dòng)脈瘤破裂的實(shí)際情況,本文提出FWASTD,經(jīng)過(guò)與非時(shí)滯系統(tǒng)的對(duì)比,說(shuō)明了時(shí)滯的有效性,通過(guò)理論證明和數(shù)值仿真,論證了其豐富的混沌性質(zhì)并分析了時(shí)滯對(duì)系統(tǒng)造成的影響.其中尤為重要的是,研究表明FWASTD系統(tǒng)在一定條件下隨血流阻力的增加而穩(wěn)定,與臨床上的促進(jìn)血栓形成來(lái)輔助治療[10]形成較為明確的對(duì)應(yīng)關(guān)系,說(shuō)明血流阻力系數(shù)研究對(duì)于臨床診斷具有一定意義.

同時(shí),利用分?jǐn)?shù)階時(shí)滯系統(tǒng)的穩(wěn)定性理論設(shè)計(jì)了合適的線性控制器以對(duì)FWASTD進(jìn)行有效控制及同步控制.本文提出的FWASTD對(duì)腦動(dòng)脈瘤里的時(shí)滯研究提供了理論基礎(chǔ),在一定程度上為腦動(dòng)脈血管瘤的臨床診斷和治療提出了理論指導(dǎo).