高飛 胡道楠 童恒慶 王傳美

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

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

1 引 言

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

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

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

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

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

2 理論知識

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

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

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

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

3 FWASTD及其性質(zhì)

3.1 建立FWASTD

從文獻[17]可知,FWAS如下:

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

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

3.2 FWASTD的有效性

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

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

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

3.3 FWASTD的取值和仿真

3.3.1 FWASTD的分數(shù)階次取值

FWASTD的分數(shù)階次是整數(shù)階次的推廣.分數(shù)階算子本身具有記憶性而優(yōu)于整數(shù)階算子,符合刻畫系統(tǒng)生理病情的需要.為了觀察分數(shù)階次取值對于系統(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)對應第一組組合,圖2(b)、圖2(d)、圖2(f)對應第二組組合,其他參數(shù)參照3.2節(jié)進行取值.

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

3.3.2 FWASTD的時滯取值

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

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

3.3.3 FWASTD的仿真

FWASTD作為整數(shù)階WAS的推廣,亦能刻畫遲發(fā)性腦動脈瘤系統(tǒng)的混沌狀態(tài),下面針對FWASTD(系統(tǒng)(2))進行仿真驗證.取τ=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)在時滯狀態(tài)的時間序列圖、相圖和Poincaré截面.

圖2 系統(tǒng)在給定初值下,q1=(a)0.8,(c)0.975,(e)1時q2的分岔圖和最大Lyapunov指數(shù)圖以及q2=(b)0.75,(d)0.95,(f)1.15時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 時滯影響下系統(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)尖峰,說明系統(tǒng)呈現(xiàn)出混沌狀態(tài);從圖4(c)可以看出,其軌跡無規(guī)律;再結(jié)合圖4(d)中系統(tǒng)(3)的Poincaré截面里具有層次結(jié)構(gòu)且成片密集的點,亦證明了系統(tǒng)(3)處于混沌狀態(tài),說明階次為分數(shù)的FWASTD也可刻畫系統(tǒng)的混沌狀態(tài).

圖4 系統(tǒng)在給定初值下,(a)x-t,y-t的時間歷程圖,(b)x-t的時間歷程圖,(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ù)對系統(tǒng)的影響

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

多數(shù)文獻[14,16,17]把研究的重點放在生理參量脈沖壓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ù)圖可證實系統(tǒng)會隨著血流阻力的增加而穩(wěn)定,臨床上亦有促進血栓形成來輔助治療的記載[14],以上種種均說明時滯狀態(tài)下研究血流阻力系數(shù)對于臨床診斷具有重大意義.

4 FWASTD的控制

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

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

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

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

其中,

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

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

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

故G(t,X(t))除t外所有變元滿足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é)果可以看出相圖軌線無規(guī)律,Poincaré截面具有一定層次結(jié)構(gòu)和形狀,因此具有初值條件的系統(tǒng)(5)不處于穩(wěn)定狀態(tài).

對系統(tǒng)(5)設計線性控制器:

定理4 當k1>0.5,0 6 k26 0.9,k3>1.95時,系統(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隨著時間增長在[?0.2,0.1]波動,仿真結(jié)果證明,該控制器可以把FWASTD控制在一個比較穩(wěn)定的范圍之內(nèi),雖然有周期性波動,但是波動振幅小(上下不超過0.1),滿足混沌控制需要.圖6(c)和圖6(d)也說明了具有線性控制器的FWASTD處于比較穩(wěn)定的狀態(tài),體現(xiàn)在實際中表現(xiàn)為血管中血流速度和加速度基本保持穩(wěn)定,病情處于可控狀態(tài).近年來,在血管內(nèi)使用血流導向裝置成為治療腦動脈瘤的一種主流辦法.血流導向裝置可以通過改變血流速度和加速度等血流動力學因子來改善動脈瘤內(nèi)部的血流情況從而達到治療疾病的目的,但現(xiàn)有的血流導向裝置在設計和材料使用上仍有缺陷[29],本文為其提供了相關(guān)的設計參考.

4.3 FWASTD的同步控制

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

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

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

誤差系統(tǒng)為

數(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表示時間,誤差系統(tǒng)仿真結(jié)果如圖7所示.

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

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

5 結(jié) 論

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

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