改進(jìn)雙鈣庫(kù)模型的分岔分析*

牛帥，孫小強(qiáng)，祁宏

(1. 山西大學(xué)復(fù)雜系統(tǒng)研究所，山西 太原 030006；2. 中山大學(xué)中山醫(yī)學(xué)院，廣東 廣州 510089)

鈣離子(Ca2+)是細(xì)胞內(nèi)最常使用的第二信使[1]，可以調(diào)控幾乎一切生理過(guò)程[2]，例如，肌肉細(xì)胞收縮、神經(jīng)元細(xì)胞釋放神經(jīng)遞質(zhì)、新陳代謝、學(xué)習(xí)記憶以及細(xì)胞的成熟、分化和凋亡等[3-4]。Ca2+信號(hào)之所以能調(diào)節(jié)如此多的生命活動(dòng)，歸結(jié)于它擁有復(fù)雜的時(shí)空動(dòng)力學(xué)特性[3-4]。當(dāng)細(xì)胞處于靜息狀態(tài)時(shí)，細(xì)胞質(zhì)中Ca2+濃度([Ca2+])很低，大約在0.1 μmol/L左右，細(xì)胞核和線粒體中[Ca2+]與細(xì)胞質(zhì)接近；細(xì)胞內(nèi)的Ca2+大多數(shù)儲(chǔ)藏于其它細(xì)胞器中，這些細(xì)胞器被稱為鈣庫(kù)。當(dāng)細(xì)胞受到刺激時(shí)，Ca2+從鈣庫(kù)中釋放出來(lái)，細(xì)胞質(zhì)中[Ca2+]迅速達(dá)到1 μmol/L以上，隨后被鈣泵重新攝入鈣庫(kù)中，形成鈣振蕩[5-6]。Ca2+正是通過(guò)其振蕩的振幅、周期和持續(xù)時(shí)間等動(dòng)力學(xué)特性來(lái)控制各種生理活動(dòng)[7-8]。

細(xì)胞中主要的鈣庫(kù)為內(nèi)質(zhì)網(wǎng)(在肌肉細(xì)胞中被稱為肌質(zhì)網(wǎng))[6]。當(dāng)內(nèi)質(zhì)網(wǎng)釋放Ca2+時(shí)，線粒體會(huì)將大量Ca2+暫時(shí)吸收[9]，在短時(shí)間內(nèi)發(fā)揮鈣庫(kù)的作用[10]。內(nèi)質(zhì)網(wǎng)膜上有兩種不同的Ca2+釋放通道：IP3R (inositol-1,4,5-triphosphate receptor)和RyR (ryanodine receptor)[11]。IP3R對(duì)IP3分子敏感[12]，IP3分子與IP3R結(jié)合可使其釋放Ca2+，稱為IP3誘導(dǎo)鈣釋放(IP3-induced calcium release)[13]。RyR對(duì)IP3分子不敏感但對(duì)Ca2+敏感，Ca2+可促使RyR釋放Ca2+，稱為鈣致鈣釋放(calcium-induced calcium release)[14]。1990年Goldbeter等[15]根據(jù)這些生物學(xué)事實(shí)提出關(guān)于鈣振蕩的雙鈣庫(kù)模型，該模型中雙鈣庫(kù)分別指對(duì)IP3敏感的鈣庫(kù)(IP3-sensitive Ca2+store, IS)和對(duì)Ca2+敏感的鈣庫(kù)(Ca2+-sensitive Ca2+store, CS)。

有關(guān)該模型的研究可分為兩大類：一部分將此模型擴(kuò)展為隨機(jī)模型，研究噪聲對(duì)鈣振蕩的影響[16-20]；另一部分從數(shù)值模擬方面分析比較了不同刺激下鈣振蕩是否存在及其振蕩頻率的變化[21-23]。我們注意到雙鈣庫(kù)模型[15]產(chǎn)生鈣振蕩的核心機(jī)制是合作性與正反饋，該機(jī)制在數(shù)學(xué)上體現(xiàn)為希爾方程，本研究想對(duì)其希爾系數(shù)和激活常數(shù)進(jìn)行詳細(xì)分析，以了解它們?cè)诋a(chǎn)生鈣振蕩方面的作用。本研究在降低原始雙鈣庫(kù)模型希爾系數(shù)的基礎(chǔ)上提出了改進(jìn)模型，運(yùn)用定性理論分析了平衡點(diǎn)的存在性和穩(wěn)定性，并利用分岔理論[24]進(jìn)一步分析了分岔點(diǎn)的類型，還結(jié)合參數(shù)敏感性分析和雙參數(shù)分岔分析，研究了關(guān)鍵參數(shù)對(duì)鈣振蕩的影響，最后根據(jù)生物學(xué)意義提出了相關(guān)生理學(xué)研究的建議。

1 數(shù)學(xué)模型

1.1 原始雙鈣庫(kù)模型

雙鈣庫(kù)模型[15]提出的鈣振蕩機(jī)制可分為四個(gè)階段(圖1)：①IP3形成階段，外界刺激，如激動(dòng)劑與細(xì)胞膜上的G蛋白耦聯(lián)受體結(jié)合，使細(xì)胞內(nèi)產(chǎn)生IP3分子。②IP3誘導(dǎo)鈣釋放階段，IP3分子誘導(dǎo)IS上的Ca2+釋放通道打開(v1)，釋放出的Ca2+在細(xì)胞質(zhì)內(nèi)聚集，構(gòu)成原發(fā)Ca2+；即使在無(wú)刺激狀態(tài)下，也會(huì)有少量Ca2+從細(xì)胞外漏入細(xì)胞質(zhì)中(v0)。③鈣致鈣釋放階段，原發(fā)Ca2+被位于CS上的鈣泵泵入CS中(v2)，當(dāng)CS中Ca2+積累到一定量時(shí)，原發(fā)Ca2+還可引發(fā)CS上的Ca2+釋放通道打開(v3)，釋放出的Ca2+是構(gòu)成鈣振蕩(也稱鈣尖峰)上升相的主要因素；同時(shí)CS中Ca2+可漏出到細(xì)胞質(zhì)中(kf)。④鈣尖峰恢復(fù)階段，位于細(xì)胞膜上的鈣泵將細(xì)胞質(zhì)中的Ca2+泵出細(xì)胞外(k)，構(gòu)成鈣尖峰的下降相。

雙鈣庫(kù)模型有兩點(diǎn)假設(shè)：(i)細(xì)胞外Ca2+的快速補(bǔ)充使得IS中Ca2+的量恒定。(ii) IS釋放Ca2+的量與IP3分子濃度成正比，而后者與刺激強(qiáng)度成正比。模型中β表示外界因素對(duì)細(xì)胞的刺激強(qiáng)度，其取值范圍為0～1。

基于以上機(jī)制，Goldbeter等建立了二維模型，分別用Z和Y表示細(xì)胞質(zhì)和CS中[Ca2+]，其表達(dá)式為：

(1)

在文獻(xiàn)[15]中，當(dāng)刺激強(qiáng)度β為30%時(shí)，鈣振蕩的振幅為1.5 μmol/L，周期為1 s。系統(tǒng)存在兩個(gè)霍普夫分岔點(diǎn)(Hopf bifurcation point，簡(jiǎn)記為HB)：β1=0.291和β2=0.775，當(dāng)0.291<β<0.775時(shí)，Ca2+產(chǎn)生振蕩，振幅變化范圍較小，而周期變化范圍較大。

圖1 雙鈣庫(kù)模型機(jī)制圖，修改自文獻(xiàn)[15]Fig.1 The mechanism of two-pool model, which is modified from reference[15]

1.2 改進(jìn)雙鈣庫(kù)模型

原始雙鈣庫(kù)模型中希爾系數(shù)的取值為m=n=2，p=4，這些系數(shù)的取值并無(wú)嚴(yán)格的生物學(xué)意義，反而不利于數(shù)學(xué)分析。為了方便分析模型平衡點(diǎn)的類型及穩(wěn)定性，將原始模型中的希爾系數(shù)降低為m=n=1，p=2，相應(yīng)地將k的取值由10調(diào)整為20，于是得到改進(jìn)模型：

(2)

1.3 參數(shù)敏感性分析

在原始雙鈣庫(kù)模型參數(shù)值[15]的基礎(chǔ)上對(duì)各參數(shù)值進(jìn)行±10%的變化，分析參數(shù)敏感性[25]。根據(jù)文獻(xiàn)[26]中的參數(shù)敏感性分析理論，局部敏感性系數(shù)(Slocal)的計(jì)算公式為：

p和△p分別指參數(shù)值和參數(shù)值的變化量，output和△output是相應(yīng)的鈣振蕩的振幅或者周期和它們的變化量。若Slocal>0，表示隨著參數(shù)值的增加(減少)，鈣振蕩的振幅或者周期也增加(減少)；若Slocal<0，表示隨著參數(shù)值的增加(減少)，鈣振蕩的振幅或者周期減少(增加)。

2 結(jié) 果

2.1 原始雙鈣庫(kù)模型參數(shù)敏感性分析

為了確定原始雙鈣庫(kù)模型中各參數(shù)在鈣振蕩中發(fā)揮作用的大小，進(jìn)行了參數(shù)敏感性分析，結(jié)果見(jiàn)圖2。本研究規(guī)定：若某參數(shù)對(duì)應(yīng)的Slocal的絕對(duì)值小于1，則該參數(shù)為低敏感參數(shù)；若其絕對(duì)值大于2，則為高敏感參數(shù)；介于兩者之間的為中等敏感參數(shù)。對(duì)于鈣振蕩的振幅來(lái)說(shuō)，CS上鈣泵的激活常數(shù)K2和CS上鈣釋放通道的激活常數(shù)KA為高敏感參數(shù)；對(duì)于鈣振蕩的周期來(lái)說(shuō)，KA為高敏感參數(shù)。

2.2 原始與改進(jìn)雙鈣庫(kù)模型的時(shí)間序列比較

文獻(xiàn)[15]中提到即使模型中的三個(gè)希爾系數(shù)m、n、p均取1，仍會(huì)得到振蕩結(jié)果，但是，本研究發(fā)現(xiàn)這需要同時(shí)改變其它多個(gè)參數(shù)，且其分岔圖只有一個(gè)分岔點(diǎn)(結(jié)果未給出)。當(dāng)m=n=1，p=2時(shí)，僅需改變一個(gè)中等敏感參數(shù)，即將Ca2+流出細(xì)胞的速率常數(shù)k的值由10變?yōu)?0，便可得到與原模型相似的結(jié)果。原始模型與改進(jìn)模型中細(xì)胞質(zhì)鈣振蕩的時(shí)間序列圖分別見(jiàn)圖3(a)和3(b)。為了便于分析，在改進(jìn)模型中取m=n=1，p=2。

2.3 改進(jìn)雙鈣庫(kù)模型的數(shù)學(xué)計(jì)算

根據(jù)定性理論對(duì)系統(tǒng)(2)求解可知Z0=0.365β+0.05，Y0由a1Y2+a2Y+a3=0確定。其中：

a1=389 017β3+1 225 670β2+

2 679 100β+6 825 000，

a2=170 000 429β3+

604 894 790β2+7 146 700β+3 025 000，

a3=-50 572 210β3-

20 783 100β2-310 323 000β-4 250 000

顯然a1>0且a3<0，則系統(tǒng)(2)僅存在一個(gè)正平衡點(diǎn)E0(Z0,Y0+)，記作E0(Z0,Y0)。系統(tǒng)(2)在平衡點(diǎn)E0處的雅可比矩陣為：

(3)

(4)

圖2 參數(shù)改變±10%時(shí)關(guān)于(a)振幅和(b)周期的敏感性分析Fig.2 Sensitivity analysis of (a) amplitude and (b) period when parameters are varied ±10% from standard parameters

圖3 β=0.5時(shí)，(a)原始模型與(b)改進(jìn)模型的鈣振蕩時(shí)間序列圖Fig.3 The time series of [Ca2+]Cyt in (a) original model and (b) improved model at β=0.5

將E0(Z0,Y0)代入(4)式中得到Trace(J)的符號(hào)由-A-BC決定，其中：

A=4707.048 9(β-0.088 8)·

(β-8.122 5)(β2+9.319 8β+108.679 8)·

(β2+7.094 7β+12.589 8)·

(β2+5.138 4β+7.988 9)·

(β2+0.373 4β+0.035 5)·

(β2-0.298 8β+0.067 5)，

B=-0.027 6(β+3.555 9)(β+0.191 6)·

(β-8.124 1)(β2+9.321 8β+108.717 5)·

(β2+5.138 4β+7.988 9)·

(β2-0.220 2β+0.023 3)，

C= [2.897 9·1010(β2+7.077 9β+12.536 3)·

分析可得：當(dāng)0<β<0.176 8或0.899 3<β<1時(shí)，Trace(J)<0，平衡點(diǎn)穩(wěn)定；0.176 8<β<0.899 3時(shí)，Trace(J)>0，平衡點(diǎn)不穩(wěn)定。當(dāng)β=0.176 8或0.899 3時(shí)，Trace(J)=0，系統(tǒng)平衡點(diǎn)為非雙曲平衡點(diǎn)。根據(jù)分岔理論對(duì)這兩個(gè)平衡點(diǎn)進(jìn)行詳細(xì)分析討論。

1)當(dāng)β=0.176 8時(shí)，平衡點(diǎn)為E1(0.114 5, 2.363 6)，作線性變換為Z1=Z-0.114 5，Y1=Y-2.363 6。引入相似變換

(5)

可得到新系統(tǒng)

(6)

其中：

f1(u,v)=-1.033 1u2-11.867 6uv-25.905 0v2-

0.798 9u3-4.553 1u2v-8.061 7uv2-4.272 2v3,

g1(u,v)=-3.409 1u2-39.159 4uv-85.478 8v2-

2.636 0u3-15.023 7u2v-26.601 4uv2-14.096 9v3

由Hopf分岔定理可知：第一李雅普諾夫系數(shù)

fuuguu+fvvgvv))|(0,0,β)=-34.467 3<0

其中ω=-6.061 2。

當(dāng)β=0.176 8時(shí)，因?yàn)锳<0，所以，系統(tǒng)(2)在平衡點(diǎn)E1處發(fā)生超臨界分岔。隨著參數(shù)β的增加，當(dāng)β>0.176 8時(shí)，平衡點(diǎn)由穩(wěn)定變?yōu)椴环€(wěn)定，同時(shí)在平衡點(diǎn)附近的小鄰域內(nèi)產(chǎn)生一個(gè)穩(wěn)定的周期軌，即系統(tǒng)產(chǎn)生振蕩。

2)當(dāng)β=0.899 3時(shí)，平衡點(diǎn)為E2(0.378 2, 0.596 5)。根據(jù)上述方法變換得到新系統(tǒng)：

(7)

其中：f2(u,v)=34.972 3u2+11.797 2uv-13.181 2v2-38.89u3-67.994 7u2v-42.914 2uv2-10.622 4v3,

g2(u,v)=32.425 3u2+10.938uv-12.221 2v2-36.057 7u3-63.042 7u2v-39.788 8uv2-9.848 8v3同樣由Hopf分岔定理可知：第一李雅普諾夫系數(shù)A=-20.455 2<0，其中ω=-21.571 0。

當(dāng)β=0.899 3時(shí)，因?yàn)锳<0，所以，系統(tǒng)(2)在平衡點(diǎn)E2處發(fā)生超臨界分岔。隨著參數(shù)β的增加，當(dāng)β>0.899 3時(shí)，平衡點(diǎn)由不穩(wěn)定變?yōu)榉€(wěn)定，同時(shí)在平衡點(diǎn)附近的小鄰域內(nèi)有一個(gè)穩(wěn)定周期軌的消失，即系統(tǒng)振蕩消失。

2.4 原始模型和改進(jìn)模型的比較

通過(guò)數(shù)學(xué)分析發(fā)現(xiàn)，與原始模型相比，改進(jìn)模型的HB分岔點(diǎn)發(fā)生改變，即系統(tǒng)產(chǎn)生振蕩或者振蕩消失的點(diǎn)發(fā)生改變。兩個(gè)模型的動(dòng)力學(xué)行為的數(shù)值模擬見(jiàn)圖4。

圖4(a)顯示的是β為分岔參數(shù)時(shí)兩個(gè)模型分岔圖的比較，其中黑色實(shí)線表示穩(wěn)定的平衡點(diǎn)，紅色虛線表示不穩(wěn)定的平衡點(diǎn)，綠色和橙色實(shí)心點(diǎn)分別表示原始模型和改進(jìn)模型的分岔圖。HB1和HB2為原始模型的左右分岔點(diǎn)，HB3和HB4為改進(jìn)模型的左右分岔點(diǎn)。從圖4(a)中可以看出改進(jìn)模型與原始模型有3點(diǎn)區(qū)別，即改進(jìn)模型的兩個(gè)分岔點(diǎn)分別向兩側(cè)移動(dòng)，因此，鈣振蕩區(qū)間更大；隨著β的增加鈣振蕩時(shí)的振幅減小更明顯；[Ca2+]穩(wěn)態(tài)值更低。圖4(b)為原始模型(綠線)和改進(jìn)模型(橙線)在振蕩區(qū)間的周期變化，在同等刺激強(qiáng)度下，改進(jìn)模型的振蕩周期稍微變短。總之，改進(jìn)模型與原始模型的動(dòng)力學(xué)行為雖稍有差別，但整體趨勢(shì)相同，可以較好地體現(xiàn)原始模型的性質(zhì)。

圖4 原始模型和改進(jìn)模型的(a)振幅和(b)周期比較Fig.4 Comparison of (a) amplitude and (b) period between the original model and the improved model

2.5 改進(jìn)模型中β和k的雙參數(shù)分岔分析

原始模型和改進(jìn)模型的區(qū)別不僅體現(xiàn)在希爾系數(shù)上，而且Ca2+流出細(xì)胞的速率常數(shù)k也不同。這說(shuō)明鈣振蕩不僅受刺激強(qiáng)度β的影響，還和k有一定的關(guān)系。于是我們通過(guò)對(duì)β和k的雙參數(shù)分岔分析，得到了系統(tǒng)產(chǎn)生極限環(huán)振蕩的參數(shù)區(qū)域(圖5陰影部分)。該區(qū)域之外為非振蕩區(qū)域，其穩(wěn)定吸引子為相應(yīng)參數(shù)組合下系統(tǒng)的穩(wěn)定平衡點(diǎn)。

結(jié)果表明：當(dāng)k<2.085 8或k>40.198 7時(shí)系統(tǒng)并無(wú)振蕩行為；當(dāng)2.085 8

圖5 β和k雙參數(shù)分岔圖Fig.5 Two-parameter bifurcation diagram for β and k

2.6 改進(jìn)模型中β和KA、β和K2的雙參數(shù)分岔分析

參數(shù)敏感性分析結(jié)果(圖2)表明，CS上鈣釋放通道的激活常數(shù)KA和CS上鈣泵的激活常數(shù)K2為高敏感參數(shù)，它們對(duì)鈣振蕩的振幅和周期的影響較大。因此，我們分別以β和KA、β和K2為雙分岔參數(shù)，分析了不同刺激強(qiáng)度β下KA和K2對(duì)Ca2+是否振蕩以及振蕩范圍的影響(圖6)。圖6(a)顯示，隨著KA的增大，鈣振蕩行為、Hopf分岔點(diǎn)的個(gè)數(shù)、振蕩所需刺激強(qiáng)度的變化方式與k相似。圖6(b)顯示，K2不斷增大導(dǎo)致鈣振蕩區(qū)域逐漸減小直至消失。因此，KA和K2對(duì)鈣振蕩均有影響，但與k相比，它們的數(shù)值變化范圍更小，對(duì)鈣振蕩的影響更大，這與參數(shù)敏感性分析的結(jié)果是一致的。

圖6 (a) β和KA雙參數(shù)分岔圖；(b) β和K2雙參數(shù)分岔圖Fig.6 Two-parameter bifurcation diagram for (a) β and KA,(b) β and K2

3 結(jié) 論

Ca2+作為細(xì)胞內(nèi)最重要的信使分子，在將外界信號(hào)傳遞到細(xì)胞內(nèi)的過(guò)程中發(fā)揮著重要作用[5]。雙鈣庫(kù)模型是關(guān)于鈣振蕩的經(jīng)典模型，它描述了外界刺激如何引發(fā)鈣振蕩的詳細(xì)機(jī)制[15]。但是，其核心機(jī)制鈣致鈣釋放的數(shù)學(xué)表達(dá)式中希爾方程的希爾系數(shù)太高，不便于數(shù)學(xué)分析。我們?cè)诮档推湎栂禂?shù)的基礎(chǔ)上提出了改進(jìn)模型，并首次進(jìn)行了數(shù)學(xué)分析。結(jié)合數(shù)值模擬結(jié)果，我們發(fā)現(xiàn)改進(jìn)雙鈣庫(kù)模型不僅能夠較好地反映原始雙鈣庫(kù)模型的動(dòng)力學(xué)性質(zhì)，其鈣振蕩范圍更廣，振幅變化更明顯，更好地反映了Ca2+編碼細(xì)胞信號(hào)的方式[27]。

通過(guò)參數(shù)敏感性分析和雙參數(shù)分岔分析，我們發(fā)現(xiàn)Ca2+流出細(xì)胞的速率常數(shù)k、CS上鈣釋放通道的激活常數(shù)KA、CS上鈣泵的激活常數(shù)K2均可對(duì)鈣振蕩產(chǎn)生影響，后兩者的影響尤為明顯。這些參數(shù)都有一定的生理學(xué)意義，KA和K2值的改變可能是由于RyR和鈣泵蛋白質(zhì)表達(dá)水平的改變所導(dǎo)致，這會(huì)改變Ca2+的動(dòng)力學(xué)從而引發(fā)相關(guān)疾病[28-29]。此外，這些參數(shù)之間相互也有關(guān)系，如k和KA對(duì)鈣振蕩的影響相似，因此作用于它們相應(yīng)蛋白的藥物便有異曲同工之處，這可為從系統(tǒng)水平調(diào)節(jié)鈣信號(hào)提供理論依據(jù)。

綜上所述，我們的工作不僅可以使人們更好地從數(shù)學(xué)層面研究雙鈣庫(kù)模型，對(duì)生物醫(yī)藥學(xué)者研究鈣信號(hào)及其相關(guān)疾病也有一定的指導(dǎo)意義。