跳转至

第 7 章:钙动力学(Calcium Dynamics)

7.1 Calcium Oscillations and Waves

Ca²⁺ 动力学研究的核心动因来自实验观察到的丰富时空行为:许多细胞类型在激动剂刺激下呈现 [Ca²⁺] 的振荡。Keener 把这些振荡分成两大类——一类依赖膜电位的周期性波动和随之而来的经电压门控通道的 Ca²⁺ 内流(例如心肌细胞),另一类在电压钳下发生。本章聚焦后者,其中又可按振荡所依赖的释放通道(RyR 还是 IPR)进一步区分。本章从 IPR 依赖的模型开始讨论。作者特别强调 IP₃ 依赖振荡的一个核心事实:即便在无胞外 Ca²⁺ 条件下振荡仍能持续一段时间(虽然形状和周期会变化),因此必然涉及 Ca²⁺ 在胞内储库之间的往返转运。

Berridge 等人(2003)的综述把 Ca²⁺ 振荡与大量生理过程联系起来,包括受精时卵母细胞激活、轴突生长、细胞迁移、基因表达、植物根毛结节的形成、肌肉发育、上皮细胞因子释放等。多数情形下振荡充当频率编码信号,使细胞既能以 Ca²⁺ 作为第二信使,又避免长时间高 [Ca²⁺] 的毒性效应,例如促性腺激素释放激素刺激下的胞吐依赖振荡频率(Tse 等 1993),Dolmetsch 等 1998、Li 等 1998 也证明基因表达可被 Ca²⁺ 尖峰频率调制。然而仍有不少情形下振荡所携带的信号尚未被明确解码。

IP₃ 依赖振荡的周期从几秒到几分钟不等。某些细胞类型振荡在恒定 [IP₃] 下发生,另一些则由 [IP₃] 自身的振荡驱动。磷酸酶、激酶对 IPR 的修饰也影响周期,胞外 Ca²⁺ 内流是另一重要调控机制。机制多样性意味着不可能指望单一模型涵盖所有细胞类型,但大多数模型有共同结构,作者因此主张通过研究少数模型学习总体方法论。Berridge 等提出的 "Ca²⁺ 信号工具包"(receptors, G 蛋白, channels, buffers, pumps, exchangers 等)构成模块化建模框架——每种工具包组件都有对应的建模组件,组合它们即可构造全细胞模型。

Ca²⁺ 振荡在细胞内常常组织成重复的胞内波(Rooney and Thomas 1993; Thomas 等 1996; Røttingen and Iversen 2000; Falcke 2004)。最著名的例子来自 Xenopus 卵母细胞:Lechleiter 与 Clapham 等 1991-1992 年用钙敏染料负载卵母细胞并光解释放 IP₃,通过共聚焦显微镜观察到了高度空间组织化的胞内 Ca²⁺ 波——同心圆、平面波、螺旋波等典型图样(图 7.3A)。该实验得以实现的关键是卵母细胞巨大(直径 >600 μm,比一般细胞大一个数量级);典型 Ca²⁺ 波宽约 100 μm,在小细胞中无法观察全貌,也无法形成螺旋。螺旋波的形成使卵母细胞成为 Ca²⁺ 波研究的重要系统,但外推到其他细胞类型必须谨慎。

另一经典胞内 Ca²⁺ 波是受精后沿卵皮层传播的波(Ridgway 等 1977; Nuccitelli 等 1993),它推动了最早的 Ca²⁺ 波传播模型(Gilkey 等 1978; Cheer 等 1987; Lane 等 1987),但因这些早期模型未考虑底层生理学,已被 Wagner 等 1998 和 Bugrim 等 2003 的工作取代。

Ca²⁺ 波不仅在单细胞内传播,还可经由细胞间隙连接跨越多个细胞形成胞间波(Sanderson 等 1990, 1994; Charles 等 1991, 1992; Cornell-Bell 等 1990; Kim 等 1994; Robb-Gaspers and Thomas 1995),距离可达数个细胞长度。并非所有胞间协调都如此长程——胰腺或腮腺腺泡细胞(Yule 等 1996)或肝细胞多倍体(Tordjmann 等 1997, 1998)中常见小群细胞的同步振荡。

胞内 Ca²⁺ 波的精确传播机制存在争议(且因细胞类型而异),但在许多细胞类型中公认胞内 Ca²⁺ 波由 Ca²⁺ 在释放位点间的扩散所驱动:某一群释放位点(通常 IPR 或 RyR)释放的 Ca²⁺ 扩散到邻近位点并触发进一步释放,重复此过程形成推进的高 Ca²⁺ 浓度波前——即主动传播的 Ca²⁺ 波。其理论已在第 6 章给出。但若底层 Ca²⁺ 动力学本身是振荡的(存在稳定极限环),波也可通过运动学或相位波机制传播,即波源于局部振荡器的空间有序触发;这种波不依赖 Ca²⁺ 扩散而存在,仅依赖细胞一端与另一端振荡的相位差,Ca²⁺ 扩散的作用是同步局部振荡器,相位波在无扩散耦合时也能维持。

研究 Ca²⁺ 振荡和波最常用的方法是假设底层机制是确定性的,但实验数据表明 Ca²⁺ 动力学本质上是随机的。最突出的随机事件是小尺度局部释放事件——puff 或 spark——它们由单个或少数 Ca²⁺ 释放通道的开放引起,是全局事件的基本构件。Ca²⁺ puff(IPR 介导)在 Xenopus 卵母细胞与 HeLa 细胞中被广泛研究(Marchant 等 1999; Sun 等 1998; Callamaras 等 1998; Marchant and Parker 2001; Thomas 等 2000; Bootman 等 1997a,b);Ca²⁺ spark(RyR 介导)则在 Cheng 等 1993 首次发现后被大量作者研究(Smith 等 1998; Izu 等 2001; Sobie 等 2002; Soeller and Cannell 1997, 2002)。研究 puff 和 spark 的随机特性需要随机模型——一般基于 IPR 和 RyR 的 Markov 模型的随机仿真——本章第 7.6 节对此简述。Falcke 2004 年约 200 页的综述是迄今最全面的 Ca²⁺ 动力学模型综述。

7.2 Well-Mixed Cell Models: Calcium Oscillations

假定细胞充分混合,每种物种的浓度在空间上均匀。记胞质自由 Ca²⁺ 浓度为 \(c = c(t)\)(无空间依赖),内质网中 Ca²⁺ 浓度为 \(c_e\)(均匀)。\(c\)\(c_e\) 的微分方程由 Ca²⁺ 守恒给出。文字表述:"总钙的变化率 = 进入该区室的净通量",数学形式为

\[ \frac{\mathrm{d}(vc)}{\mathrm{d}t} = \tilde{J}_{\text{net}}, \]

其中 \(v\) 为区室体积,\(\tilde{J}_{\text{net}}\) 为进入该区室的净 Ca²⁺ 通量,单位为 mol/s。通常假定细胞及胞内区室体积恒定,守恒方程简化为

\[ \frac{\mathrm{d}c}{\mathrm{d}t} = J_{\text{net}} = \tilde{J}_{\text{net}}/v, \]

其中 \(J_{\text{net}}\) 为单位胞质体积单位时间的净通入通量,单位为浓度/秒。细胞体积变化情形(如细胞体积控制模型、Ca²⁺ 振荡调控的流体分泌模型)必须用 \(\tilde{J}_{\text{net}}\) 而非 \(J_{\text{net}}\)

具体地,根据图 7.1 所示通量:

\[ \frac{\mathrm{d}c}{\mathrm{d}t} = J_{\text{IPR}} + J_{\text{RyR}} + J_{\text{in}} - J_{\text{pm}} - J_{\text{serca}} - J_{\text{on}} + J_{\text{off}} + J_{\text{uni}} - J_{\text{mito}}. \]

ER 中 Ca²⁺ 方程必须考虑胞质与 ER 体积不同:

\[ \frac{\mathrm{d}c_e}{\mathrm{d}t} = \gamma\,(J_{\text{serca}} - J_{\text{IPR}} - J_{\text{RyR}}) + J_{\text{off,e}} - J_{\text{on,e}}, \]

其中 \(\gamma = v_{\text{cyt}}/v_{\text{ER}}\) 为胞质与 ER 体积比。方程中各通量对应 Ca²⁺ 信号工具包的不同组件,Keener 仅讨论建模中最常用的若干;这些方程通常与描述 IPR/RyR 门控、泵/交换器动力学的附加微分方程耦合,第 7.2.6 节给出一个例子。讨论各工具包组件之前先逐一介绍它们的建模。

7.2.1 Influx

胞外 Ca²⁺ 内流通常是电压依赖的,但在非兴奋细胞中 Ca²⁺ 振荡常在恒定电压下发生,电压依赖性可忽略。此内流受多种其他因素调节,包括 Ca²⁺、IP₃、花生四烯酸、ER 中 Ca²⁺ 浓度。例如有证据表明某些细胞类型 ER 耗竭会通过 store-operated channels(SOCs)增加 Ca²⁺ 内流(Clapham 1995);也有证据表明 SOCs 仅在高激动剂浓度(ER 高度耗竭)下发挥作用,而低激动剂浓度下 Ca²⁺ 内流由花生四烯酸控制(Shuttleworth 1999)。然而精确机制未知。已知 \(J_{\text{in}}\) 随激动剂浓度增加而上升;若非如此稳态 [Ca²⁺] 将与 [IP₃] 无关(参见练习 1),而事实并非如此。一种常见简化是假设 \(J_{\text{leak}}\) 是 IP₃ 浓度 \(p\) 的线性递增函数:

\[ J_{\text{in}} = \alpha_1 + \alpha_2 p, \]

其中 \(\alpha_1\)\(\alpha_2\) 为常数。该极简模型未考虑 ER 装载对内流的影响,但对较低激动剂浓度情形可接受。当前对 Ca²⁺ 内流调控机制的认识不足以支撑更精细的模型。

7.2.2 Mitochondria

线粒体 Ca²⁺ 处理是高度复杂的过程,已构建了若干详细模型。Keener 为简洁起见不予讨论,读者可参阅 Colegrove 等 2000、Friel 2000、Falcke 等 2000、Grubelnik 等 2001、Marhl 等 2000、Schuster 等 2002、Selivanov 等 1998。线粒体的功能之一似乎是较慢地摄取和释放大量 Ca²⁺,从而倾向于调节波的尾缘、降低波幅、改变长期振荡行为;但这当然是过度简化。

7.2.3 Calcium Buffers

胞质中至少 99% 的 Ca²⁺ 与大型蛋白(称 Ca²⁺ buffer)结合。典型 buffer 包括 calsequestrin、calbindin、荧光染料及质膜本身。Ca²⁺ 缓冲及其对振荡与波的影响将在第 7.4 节详细讨论。

7.2.4 Calcium Pumps and Exchangers

Calcium ATPases. 早期 Ca²⁺ ATPase 泵模型采用 Hill 方程形式(第 1.4.4 节)。Lytton 等 1992 数据表明 ATPase 通量近似 \(c\) 的 sigmoid 函数,Hill 系数约 2,故常见模型为

\[ J_{\text{serca}} = \frac{V_p c^2}{K_p^2 + c^2}. \]

此类简单模型有若干严重缺陷:该通量不依赖 ER 中 Ca²⁺ 浓度且符号恒定,但已知 ER 中 Ca²⁺ 浓度足够高时泵可逆向运行并产生 ATP。

MacLennan 等 1997 构建了更详细的模型(图 7.5)。泵有两种基本构象 E1 与 E2。E1 构象从胞质结合两个 Ca²⁺,随即暴露磷酸化位点;磷酸化后切换为 E2 构象,Ca²⁺ 结合位点暴露于 ER 内腔且亲和力大大降低,Ca²⁺ 释放入 ER 后泵去磷酸化并回到 E1 构象完成循环。每将一个 Ca²⁺ 从胞质运至 ER,一个质子被同向运回胞质——这与 Na⁺–K⁺ ATPase 的相似性更加明显。转运循环的限速步骤是 E1-P-C2 → E2-P-C2 跃迁。该模型下通量计算较直接,与第 2 章详尽讨论的例子类似,留作习题。所得

\[ J_{\text{serca}} = \frac{c^2 - K_1 K_2 K_3 K_4 K_5 K_6 c_e^2}{\alpha_1 c^2 + \alpha_2 c_e^2 + \alpha_3 c^2 c_e^2 + \alpha_4}, \]

其中 \(\alpha_i\) 是速率常数的函数,写出全部过于冗长;\(K_i = k_{-i}/k_i\)。若 E1 构象下 Ca²⁺ 结合位点高亲和、E2 构象下低亲和,则 \(k_{-4}/k_4 \ll k_1/k_{-1}\),故 \(K_1 K_4 \ll 1\)。合理假设 \(k_6/k_{-6} = k_{-3}/k_3\)\(k_2/k_{-2} = k_5/k_{-5}\),则 \(K_2 K_3 K_5 K_6 = 1\)。于是 \(K_1 K_2 K_3 K_4 K_5 K_6 \ll 1\),泵可在 \(c_e \gg c\) 时仍支持正向通量。若这是封闭系统模型则细致平衡律会要求 \(K_1 K_2 K_3 K_4 K_5 K_6 = 1\),但此处不成立,因为循环由泵的磷酸化驱动,反应速率依赖 ATP、ADP、P 浓度,能量不断消耗(\(J_{\text{serca}} > 0\))或生成(\(J_{\text{serca}} < 0\))。详见第 2.5.1 节。

MacLennan 等原始描述中假设 Ca²⁺ 的结合与释放快速发生,可得简化模型并保留稳态通量的类似表达式。这种基于快平衡假设简化模型的方法很重要,第 2 章与练习 7 中有详尽探讨。为方便起见 Keener 在此简述推导。在图 7.5 中状态 S1 与 S2 用虚线框分组在一起,T2 与 T3 同样分组。快平衡假设给出 \(s_1 = K_1 c^2 s_2\),对 \(t_2\)\(t_3\) 也有类似关系。定义新变量 \(\bar{s}_1 = s_1 + s_2\)\(\bar{t}_2 = t_2 + t_3\),由上述关系

\[ \bar{s}_1 = s_1\left(1 + \frac{c^2}{K_1}\right) = s_2\left(1 + \frac{K_1}{c^2}\right). \]

\(\bar{S}_1\) 转换为 S3 的速率为

\[ k_2 s_2 = \frac{c^2 k_2 \bar{s}_1}{c^2 + K_1} \]

\(\bar{S}_1\) 转换为 T1 的速率为

\[ k_{-6} s_1 = \frac{K_1 k_{-6} \bar{s}_1}{K_1 + c^2} \]

。对各跃迁重复该过程得到图 7.6 所示简化模型,跃迁率函数 \(\phi\) 的形式可直观理解:\(c\) 增加时 S1 与 S2 间的平衡向 S2 偏移,S3 生成速率增加而 \(\bar{S}_1 \to T_1\) 速率下降;因此 \(\bar{S}_1 \to S_3\)\(c\) 的递增函数、\(\bar{S}_1 \to T_1\)\(c\) 的递减函数。从该图按相同方法计算稳态通量,结果仍取相同形式

\[ J_{\text{serca}} = \frac{c^2 - K_1 K_2 K_3 K_4 K_5 K_6 c_e^2}{\beta_1 c^2 + \beta_2 c_e^2 + \beta_3 c^2 c_e^2 + \beta_4}, \]

其中 \(\beta_i\) 为 (7.7) 中 \(\alpha_i\) 的极限值。质膜上的 Ca²⁺ ATPase 类似 SERCA ATPase,建模方式相同。

Calcium Exchangers. 胞质 Ca²⁺ 移除的另一重要途径是 Na⁺–Ca²⁺ 交换器——它以一个 Ca²⁺ 移出胞质换取三个 Na⁺ 内流。这对心肌细胞 Ca²⁺ 控制特别重要,将在第 12 章兴奋-收缩耦合模型语境下进一步讨论。

7.2.5 IP3 Receptors

IPR 的一项基本性质是其对 Ca²⁺ 或 IP₃ 阶跃变化呈时间依赖响应——对 IP₃ 或 Ca²⁺ 的阶跃增加,受体开放概率先升至峰值再下降至较低平台(图 7.7),这种下降称为受体的适应,因为受体对维持的 Ca²⁺ 或 IP₃ 浓度发生了适应。若在第一个阶跃之上叠加新的阶跃,受体再次出现峰值随后降至平台。因此 IPR 响应的是 [Ca²⁺] 或 [IP₃] 的变化而非绝对浓度。适应是许多生理系统的共同特征(第 16 章、19 章),常表现为快速激活后跟缓慢失活,类似于 Hodgkin–Huxley 方程中 Na⁺ 通道的行为(第 5 章)。

IPR 的适应目前被认为至少部分源自 Ca²⁺ 不仅刺激自身释放、还在较慢时间尺度上抑制释放。这种 Ca²⁺ 对 IPR 的序贯激活与失活被认为是 IP₃ 依赖 Ca²⁺ 振荡与波的机制之一,已出现多个纳入此假说的模型(Sneyd 等 1995b; Tang 等 1996; Schuster 等 2002; Falcke 2004 综述)。但如下所述几乎肯定还有其他重要机制在运作;某些近期模型提示 IPR 动力学重要性不如先前所想,ER 耗竭可能也起重要作用。Sneyd and Falcke 2005 给出 IPR 模型的详细综述,Sneyd 等 2004a 比较若干模型与实验数据。

An Eight-State IP3 Receptor Model. De Young and Keizer 1992 是最早纳入 Ca²⁺ 序贯激活/失活的 IPR 模型之一。该模型假设 IPR 由三个等价且独立的亚基组成,三亚基都处于导通态时才有 Ca²⁺ 通流。每个亚基有一个 IP₃ 结合位点、一个激活 Ca²⁺ 位点、一个失活 Ca²⁺ 位点,每位点可被占据或未占据,故亚基可有 8 种状态。状态标记为 \(S_{ijk}\)\(i,j,k\) 为 0 或 1,0 表示位点未占据、1 表示已占据;第一个下标指 IP₃ 位点、第二个指 Ca²⁺ 激活位点、第三个指 Ca²⁺ 失活位点(图 7.8)。模型共有 24 个速率常数,但因细致平衡要求并非全部独立;加两条简化假设进一步减少独立常数——第一,速率常数与激活 Ca²⁺ 是否结合无关;第二,Ca²⁺ 激活动力学独立于 IP₃ 结合与 Ca²⁺ 失活。剩下 10 个独立常数 \(k_1, \ldots, k_5\)\(k_{-1}, \ldots, k_{-5}\)。在这些假设下细致平衡成立。

记亚基处于状态 \(S_{ijk}\) 的比例为 \(x_{ijk}\)。微分方程基于质量作用,例如

\[ \frac{\mathrm{d}x_{000}}{\mathrm{d}t} = -(V_1 + V_2 + V_3), \]

其中

\[ V_1 = k_1 p\, x_{000} - k_{-1} x_{100}, \quad V_2 = k_4 c\, x_{000} - k_{-4} x_{001}, \quad V_3 = k_5 c\, x_{000} - k_{-5} x_{010}, \]

\(p\) 为 [IP₃],\(c\) 为 [Ca²⁺]。\(V_1\) 描述 IP₃ 与 IP₃ 位点的结合/解离速率,\(V_2\) 描述 Ca²⁺ 与失活位点的结合/解离速率,\(V_3\) 同理。因实验数据显示受体亚基以协同方式作用,模型假设 IPR 仅在三亚基都处于 \(S_{110}\) 状态(即各结合一个 IP₃ 和一个激活 Ca²⁺)时才有 Ca²⁺ 电流,因此受体开放概率为 \(x_{110}^3\)

图 7.9 显示 IPR 开放概率作为 [Ca²⁺] 的函数(这是模型的实验基础)。Bezprozvanny 等 1991 表明该开放概率是 [Ca²⁺] 的钟形函数——低 [Ca²⁺] 时 [Ca²⁺] 增加提高开放概率,高 [Ca²⁺] 时 [Ca²⁺] 增加降低开放概率。模型参数选取使与这些稳态数据吻合。受体动力学性质同等重要:Ca²⁺ 快速激活 IPR,但 Ca²⁺ 在较慢时间尺度上使其失活;模型中这通过速率常数大小体现(\(k_5 > k_2\)\(k_5 > k_4\))。

Reduction of the Eight-State IP3 Receptor Model. 八态受体模型(七条微分方程、众多参数)的复杂性激励人们寻找保留其本质性质的简化模型。因 IP₃ 与 IP₃ 位点结合快速、Ca²⁺ 与激活位点结合快速,可舍弃这些结合过程的瞬态细节,假设受体相对 IP₃ 结合与 Ca²⁺ 激活处于准稳态(De Young and Keizer 1992; Keizer and De Young 1994; Li and Rinzel 1994; Tang 等 1996)。表 7.1 给出的参数值已暗示这一点——\(k_1, k_3, k_5\)\(k_2, k_4\) 大很多,\(k_{-1}, k_{-3}, k_{-5}\) 也比 \(k_{-2}, k_{-4}\) 大。

如图 7.8,将受体状态分为两组:未结合失活位点 Ca²⁺ 的(\(S_{000}, S_{010}, S_{100}, S_{110}\),图 7.8 上排,组 I 状态)与结合了失活位点 Ca²⁺ 的(\(S_{001}, S_{011}, S_{101}, S_{111}\),下排,组 II 状态)。因 IP₃ 结合与 Ca²⁺ 激活结合被假设为快速过程,组内受体状态处于准稳态;但组 I 与组 II 之间的跃迁(失活位点的结合/解离)缓慢,故组 I 状态与组 II 状态不处于平衡。

按标准的简化方法(第 2.4.1 节或第 7.2.4 节),令

\[ y = x_{001} + x_{011} + x_{101} + x_{111}, \]

可得

\[ \frac{\mathrm{d}y}{\mathrm{d}t} = \frac{(k_{-4} K_1 K_2 + k_{-2} p K_4) c}{K_4 K_2 (p + K_1)}(1 - y) - \frac{k_{-2} p + k_{-4} K_3}{p + K_3} y. \]

推导留作习题(练习 4)。可写成

\[ \tau_y(c, p) \frac{\mathrm{d}y}{\mathrm{d}t} = y_\infty(c, p) - y, \]

该形式便于与其他模型(如 Hodgkin–Huxley 方程)比较。开放概率由

\[ x_{110} = \frac{p c (1 - y)}{(p + K_1)(c + K_5)} \]

给出。注意 \(1 - y\)(未被 Ca²⁺ 失活的受体比例)起失活变量作用,类似于 Hodgkin–Huxley 方程中的 \(h\)(第 5 章)。简化模型可写为

\[ x_{110} = \frac{p c}{(p + K_1)(c + K_5)} h, \quad \tau_h(c, p) \frac{\mathrm{d}h}{\mathrm{d}t} = h_\infty(c, p) - h, \]

其中 \(h = 1 - y\)\(\tau_h\)\(h_\infty\) 可从 \(y\) 的方程直接得到。

A Model with Saturating Binding Rates. 最近明确显示八态模型存在严重缺陷。实验上 IPR 开放速率仅在一个数量级内变化,而 Ca²⁺ 或 IP₃ 浓度跨好几个数量级变化;这意味着 IPR 开放不能由简单的 Ca²⁺ 或 IP₃ 结合质量作用动力学支配,而必须采用允许结合速率饱和的动力学方案。

应用第 3.5.4 节学到的内容,有一种显然的方法:将 Ca²⁺ 或 IP₃ 结合步骤与开放步骤分离,使用激动剂控制通道模型中的亲和与效能概念(第 3 章),得到结合速率随激动剂浓度升高而饱和的反应——这本质上等同于第 1.4 节讨论的酶动力学模型。由于反应速率饱和,简单质量作用动力学不再适用,故发展出 Michaelis–Menten 型模型。考虑反应方案

\[ \tilde{A} \underset{k_{-1}}{\overset{k_1 c}{\rightleftharpoons}} \bar{A} \underset{k_{-2}}{\overset{k_2}{\rightleftharpoons}} I \]

。若 \(\tilde{A}\)\(\bar{A}\) 之间的跃迁快于其他反应,二者瞬时平衡,则 \(c \tilde{A} = K_1 \bar{A}\)\(K_1 = k_{-1}/k_1\)。按常用的求慢动力学方法,得

\[ \frac{\mathrm{d}A}{\mathrm{d}t} = k_{-2} I - \phi(c) A, \]

其中 \(A = \bar{A} + \tilde{A}\)\(\phi(c) = \frac{k_2 c}{c + K_1}\)\(c\) 中饱和的速率。因此快平衡假设可将 (7.21) 简化为 \(A \underset{k_{-2}}{\overset{\phi(c)}{\rightleftharpoons}} I\) 。这样就能以简单方式将饱和结合动力学纳入模型。

Sneyd and Dufour 2002 采用这类饱和结合方案构建了 IPR 模型,基于 Taylor 等的定性模型,并与 Hajnóczky and Thomas 1997 的方案一致。除饱和结合速率外,模型的主要特征为:(1) IPR 可在无 Ca²⁺ 时被 IP₃ 开放,但电导较低;(2) IPR 可在无 IP₃ 时被 Ca²⁺ 失活;(3) 一旦 IP₃ 结合,IPR 可自发失活(至关闭态 S),与 Ca²⁺ 无关;(4) 一旦 IP₃ 结合,IPR 还可结合 Ca²⁺ 激活受体,故存在 Ca²⁺ 介导激活与自发失活的内禀竞争;(5) IPR 一旦被 Ca²⁺ 激活,可通过结合更多 Ca²⁺ 被失活(这是与 Taylor 1998 定性模型的主要分歧点);(6) IP₃ 与 Ca²⁺ 结合按序进行。

按这些假设模型如图 7.10,相应方程为

\[ \frac{\mathrm{d}R}{\mathrm{d}t} = \phi_{-2} O - \phi_2 p R + k_{-1} I_1 - \phi_1 R, \]
\[ \frac{\mathrm{d}O}{\mathrm{d}t} = \phi_2 p R - (\phi_{-2} + \phi_4 + \phi_3) O + \phi_{-4} A + k_{-3} S, \]
\[ \frac{\mathrm{d}A}{\mathrm{d}t} = \phi_4 O - \phi_{-4} A - \phi_5 A + k_{-1} I_2, \]
\[ \frac{\mathrm{d}I_1}{\mathrm{d}t} = \phi_1 R - k_{-1} I_1, \quad \frac{\mathrm{d}I_2}{\mathrm{d}t} = \phi_5 A - k_{-1} I_2, \]

其中 \(R + O + A + S + I_1 + I_2 = 1\),且

\[ \phi_1(c) = \frac{\alpha_1 c}{\beta_1 + c}, \quad \phi_3(c) = \frac{\alpha_3}{\beta_3 + c}, \quad \phi_5(c) = \frac{\alpha_5 c}{\beta_5 + c}, \]
\[ \phi_2(c) = \frac{\alpha_2 + \beta_2 c}{\beta_1 + c}, \quad \phi_{-2}(c) = \frac{\alpha_{-2} + \beta_{-2} c}{\beta_3 + c}, \quad \phi_4(c) = \frac{\alpha_4 c}{\beta_3 + c}, \quad \phi_{-4}(c) = \frac{\alpha_{-4}}{\beta_5 + c}. \]

参数由拟合实验数据确定(表 7.2)。假设 IPR 由四个相同且独立的亚基组成,所有四个亚基都处于 O 或 A 态时允许 Ca²⁺ 电流;亚基在 A 态越多电导越大。开放概率的简单表达为

\[ P_o = (a_1 O + a_2 A)^4, \]

原始模型中 \(a_1 = 0.1\)\(a_2 = 0.9\),但这些值可改而不显著影响拟合或模型行为。

Keener 顺带提醒:"参数由拟合实验数据确定"这一简单表述掩盖了许多棘手问题。第一,模型中通常并非所有参数都能由现有实验数据无歧义确定,对上述 IPR 模型确实只有少数几个参数能被较有把握地钉住——除了简化模型或收集更多适当类型的实验数据外无他法,而确定需要哪些额外实验数据本身不是平凡问题。第二,从数据确定参数的过程本身复杂——简单最小二乘拟合有严重缺陷,贝叶斯推断和马尔可夫链蒙特卡洛等更精细方法实现起来更困难。Ball 等 1999 有详细讨论。

7.2.6 Simple Models of Calcium Dynamics

至此已构建了大多数重要 Ca²⁺ 通量模型。将它们组合成 Ca²⁺ 动力学模型很直接:选好 IPR(或 RyR)的最爱模型、ATPase 与内流的最爱模型,代入 (7.3) 与 (7.4)。以简单模型为例,假设只有两种通量——IPR 与 SERCA:

\[ \frac{\mathrm{d}c}{\mathrm{d}t} = (k_f P_o + J_{er})(c_e - c) - J_{\text{serca}}. \]

由于唯一通量在 ER 与胞质之间,故

\[ \frac{\mathrm{d}c_e}{\mathrm{d}t} = -\gamma \frac{\mathrm{d}c}{\mathrm{d}t}, \]

于是 \(c + c_e/\gamma = c_t\) 不变。此类模型称为封闭细胞模型(第 7.2.7 节)。下一步用 Hill 函数(第 1 章)建模 SERCA 泵

\[ J_{\text{serca}} = \frac{V_p c^2}{K_p^2 + c^2}, \]

并以简化的 DeYoung–Keizer 模型 (7.18) 建模 IPR

\[ P_o = \left(\frac{p c (1 - y)}{(p + K_1)(c + K_5)}\right)^3, \]

其中 \(y\) 满足 (7.16)。注意 IPR 通量假设与 ER 与胞质之间的浓度差成正比;这只在 ER 膜两侧无电位差时成立——尽管 Ca²⁺ 带电,Ca²⁺ 流不引发电位差,因为有反离子自由流动维持电荷平衡。

由于 \(c_t\) 固定,该简单模型可简化为二变量模型,可构造相平面。\(p = 0.5\) 的数值解如图 7.11 所示。相平面在结构上与若干二变量激发介质模型(FitzHugh–Nagumo 方程或简化的 Hodgkin–Huxley 方程,第 5 章)相同。失活变量 \(y\) 的零线是 \(c\) 的单调递增函数;\(c\) 的零线是 \(c\) 的 N 形函数(该尺度下 N 形很浅)。两条零线有唯一交点。交点稳定性(粗略但非精确)由其所在 \(c\)-零线的分支决定——最左支上固定点稳定,中间支上不稳定。若固定点不稳定则存在稳定周期极限环,对应自发 Ca²⁺ 振荡。

A More Complex Example. 稍现实的模型如下:忽略 RyR 与线粒体通量;采用 IPR 饱和结合模型;采用四态 Markov 模型描述 SERCA 泵(每循环转运一个 Ca²⁺);用 Hill 函数描述质膜泵;假设 Ca²⁺ 缓冲快速且线性(第 7.4 节)。则得

\[ \frac{\mathrm{d}c}{\mathrm{d}t} = (k_f P_o + J_{er})(c_e - c) - J_{\text{serca}} + J_{\text{in}} - J_{\text{pm}}, \]
\[ \frac{\mathrm{d}c_e}{\mathrm{d}t} = \gamma\,[J_{\text{serca}} - (k_f P_o + J_{er})(c_e - c)], \]

其中

\[ J_{\text{serca}} = \frac{c - \alpha_1 c_e}{\alpha_2 + \alpha_3 c + \alpha_4 c_e + \alpha_5 c c_e}, \quad P_o = (0.1 O + 0.9 A)^4, \]
\[ J_{\text{pm}} = \frac{V_p c^2}{K_p^2 + c^2}, \quad J_{er} = \text{const}, \quad J_{\text{in}} = a_1 + a_2 p. \]

常数 \(J_{er}\) 表示 ER 的恒定泄漏,与稳态下 ATPase 通量平衡所需。(7.36) 与 (7.37) 与 (7.26)–(7.30) 耦合描述 IPR。作者在此采用两种泵模型以强调可混合搭配各通量的个体模型。SERCA ATPase 有五个参数、质膜 ATPase 仅两个,很可能 SERCA 泵被过度参数化——若拟合数据时这会是顾虑,因数据可能不足以无歧义确定所有参数。但作者满足于使用可由机理模型证明的函数,且参数可调至与实验数据合理吻合。参数值见表 7.4。

\(p\) 对应 IP₃ 浓度(间接对应激动剂浓度),作者描述模型行为随 \(p\) 变化。模型的分岔图如图 7.12A。\(p\) 增大时稳态 Ca²⁺ 浓度也增大(因 \(J_{\text{in}}\)\(p\) 增加),振荡出现在中间一段 \(p\) 值范围。\(p\) 取两个不同值时的典型振荡如图 7.12B、C 所示;\(p\) 增大时振荡频率上升。

鉴于 IPR、RyR、ATPase 与其他通量模型的多样性,必须思考 Ca²⁺ 动力学模型能告诉我们什么。例如通过适当选取 IPR 与 ATPase 模型,可构造出呈现各式各样 Ca²⁺ 振荡形状与大小的全细胞模型。图 7.12 中的振荡合理地描述了胰腺腺泡细胞(该模型最初为之设计的细胞类型)的 Ca²⁺ 振荡;然而参数与通量模型可(几乎无限地)调整以模拟其他细胞类型的振荡行为。因此模型仅呈现出性质合理的 Ca²⁺ 振荡这一事实对底层机制几乎无说明力。要对实际细胞中的特定机制得出任何结论还需大量工作——这超出本书范围,Falcke 2004 综述是最易入手的入口。

7.2.7 Open- and Closed-Cell Models

Ca²⁺ 振荡周期由多个因素决定——IPR/RyR 的动力学、IP₃ 生成与降解速率(第 7.2.8 节)、Ca²⁺ 尖峰后 ER 重新装载速率、跨质膜的 Ca²⁺ 转运速率。研究后两个效应的一个好办法是将模型略作改写,使膜 Ca²⁺ 通量与 ER 通量易于分离。为此引入新变量

\[ c_t = c + c_e/\gamma. \]

新变量 \(c_t\) 是细胞内(胞质加 ER)Ca²⁺ 总摩尔数除以胞质体积,故是细胞总 Ca²⁺ 含量的度量。将 (7.36)、(7.37) 用 \(c_t\) 重写:

\[ \frac{\mathrm{d}c_t}{\mathrm{d}t} = \delta (J_{\text{in}} - J_{\text{pm}}), \]
\[ \frac{\mathrm{d}c_e}{\mathrm{d}t} = \gamma\,[J_{\text{serca}} - (k_f P_o + J_{er})(c_e - (c_t - \gamma c_e))], \]

其中为方便起见在 \(c_t\) 方程中引入比例因子 \(\delta\),便于修改膜转运与 ER 转运的比例而不影响模型其他部分。当然也可改用 \(c_t\)\(c\) 表达而非 \(c_t\)\(c_e\)

若膜通量远小于 ER 通量则 \(\delta\) 小——这在许多细胞类型中成立。\(\delta = 0\) 时 Ca²⁺ 不进入也不离开细胞,\(c_t\) 保持不变。此模型称为封闭细胞模型,第 7.2.6 节已讨论过一个简单例子。\(\delta \neq 0\) 时为开放细胞模型,细胞内总 Ca²⁺ 量随时间变化。

某些细胞类型中 \(c_t\) 变化足够慢,可通过先理解封闭细胞模型行为再理解 Ca²⁺ 动力学。\(\delta\) 足够小时开放细胞模型行为与封闭细胞模型相似(因短时间尺度上 \(c_t\) 几近不变),但 \(c_t\) 缓慢漂向稳态时会受到缓慢调制。该分析依赖两个时间尺度的存在——\(c\)\(c_e\) 的快速变化受 \(c_t\) 缓慢变化调制。该方法由 Rinzel 1985 首次提出用以研究胰腺 β-细胞阵发振荡模型,第 9 章有更详细讨论。

该公式的一个特别有趣的用法来自 \(\delta\)\(c_t\) 可在实验上控制的事实(Sneyd 等 2004b)。细胞外施加高浓度 La³⁺ 同时阻断 Ca²⁺ 内流与质膜 ATPase,从而有效使 \(\delta\) 降为零。此外用光解释放性 Ca²⁺ 预先装载细胞可通过闪光调节 \(c_t\)。一旦施加 La³⁺,\(c_t\) 可视为控制参数随意增加(虽然不幸地难以减小)并观察细胞行为——这允许对模型预测进行细致检验。

7.2.8 IP3 Dynamics

上述模型中 Ca²⁺ 振荡在恒定 IP₃ 浓度下发生。IP₃ 的作用是激活 IPR;受体一旦被激活,Ca²⁺ 反馈接管,振荡周期由 Ca²⁺ 对 IPR 的反馈动力学以及膜与 ER 通量的相互作用控制。

这当然是过度简化。已知 IP₃ 生成速率依赖 Ca²⁺,且某些细胞类型中 [Ca²⁺] 振荡伴随 [IP₃] 振荡(Hirose 等 1999; Nash 等 2001; Young 等 2003)。目前尚不清楚 [IP₃] 振荡是否为 Ca²⁺ 振荡所必需——前者可能只是被动跟随后者。实验证据既不统一也不结论性。

已有若干 Ca²⁺ 动力学模型纳入 IP₃ 生成与降解对 Ca²⁺ 的依赖。早期模型有 Meyer and Stryer 1988, 1991; Swillens and Mercan 1990; De Young and Keizer 1992; Cuthbertson and Chay 1991;较近期模型有 Shen and Larter 1995; Borghans 等 1997; Dupont and Erneux 1997; Houart 等 1999; Politi 等 2006。

Ca²⁺ 影响 IP₃ 浓度有两种主要方式(图 7.13)。第一,PLC 被 Ca²⁺ 激活发生于多种细胞类型,构成 Meyer and Stryer 1988、De Young and Keizer 1992 模型的基础。第二,IP₃ 以两种主要方式降解,每种都受反馈调控——被 IP₃ 5-phosphatase 去磷酸化为 inositol 1,4-bisphosphate (IP₂),或被 IP₃ 3-kinase 磷酸化为 inositol 1,3,4,5-tetrakisphosphate (IP₄)。因 IP₄ 也是 5-phosphatase 的底物,IP₄ 充当 IP₃ 去磷酸化的竞争性抑制剂。此外 Ca²⁺(以 Ca²⁺-calmodulin 复合物形式)可增强 3-kinase 活性。

直觉上很清楚这样的复杂反馈阵列容易产生 [IP₃] 与 Ca²⁺ 的相互依赖振荡,上述模型已证实这一直觉期望。但这些模型中只有少数用于做出可被实验检验的预测。两个显著例外是 Dupont 等 2003 与 Politi 等 2006 的模型。

肝细胞中 Ca²⁺ 与 IP₃ 一起振荡。每个 Ca²⁺ 尖峰引起 IP₃ 磷酸化速率增加从而 [IP₃] 下降,随之 [Ca²⁺] 下降。然而虽然两种物种都振荡,目前并不清楚 [IP₃] 振荡是否必需。Dupont 等首先使用 Dupont and Erneux 1997 模型,其中 Ca²⁺ 影响 3-kinase 速率。该模型预测若 5-phosphatase 速率增加 25 倍、IP₃ 生成速率也增加,则具有相同定性性质的振荡仍持续,仅频率略降。在此条件下几乎无 IP₃ 经 3-kinase 降解,故 Ca²⁺ 对 3-kinase 的反馈不能起主要作用——这些振荡源自 Ca²⁺ 对 IPR 的反馈。Dupont 等随后设计了对此模型预测的实验检验:他们直接显微注射 5-phosphatase 以增加其速率,并通过提高激动剂浓度增加 IP₃ 生成速率;观察到振荡基本不变,频率按预测下降。因此当 IP₃ 降解主要经 5-phosphatase 途径(此时不受 [Ca²⁺] 直接影响)时,振荡性质几乎不变。结论是肝细胞中振荡并非源自 [Ca²⁺] 与 IP₃ 动力学的相互作用,而是源自 Ca²⁺ 对 IPR 的反馈。这一工作是建模与实验相互作用的范例。

类似研究也在 CHO 细胞(中国仓鼠卵巢细胞,一种细胞系)中完成。Politi 等构建了包含 Ca²⁺ 对 IP₃ 正负反馈的模型,然后检验外源 IP₃ 缓冲加入的后果。其分析复杂微妙,此处不便详述。结论是 CHO 细胞中 Ca²⁺ 振荡源自 Ca²⁺ 对 IP₃ 生成的正反馈。

外源 IP₃ 脉冲的响应可被用于确定任一细胞类型中 IP₃ 振荡是否为 Ca²⁺ 振荡所必需(Sneyd 等 2006)。若 IP₃ 与 Ca²⁺ 都是动态变量,必然一起振荡,则 IP₃ 脉冲应使解偏离稳定极限环,振荡恢复前有延迟——相位响应曲线应呈现相位滞后(相位提前虽理论可能但极少观察到)。若 Ca²⁺ 振荡发生在恒定 [IP₃] 下,则 IP₃ 脉冲暂时升高振荡频率。这些预测已在两种细胞类型中检验——胰腺腺泡细胞中 IP₃ 脉冲引起下一个振荡峰前的延迟;气道平滑肌细胞中 IP₃ 脉冲暂时升高振荡频率。模型据此预测胰腺腺泡细胞中 Ca²⁺ 振荡必伴随 IP₃ 振荡,而气道平滑肌细胞中则在恒定 [IP₃] 下发生。但该方法在肝细胞中给出不确定的结果(Harootunian 等 1988)且尚未广泛用于其他细胞类型,其效用尚未明确。

7.2.9 Ryanodine Receptors

细胞内储库 Ca²⁺ 释放的另一主要途径是通过 ryanodine 受体,存在于多种细胞——心肌细胞、平滑肌、骨骼肌、嗜铬细胞、垂体细胞、神经元、海胆卵等。RyR 与 IPR 在结构与功能(特别是对 Ca²⁺ 的敏感性)上有许多相似性。正如 Ca²⁺ 能激活 IPR 并增加 Ca²⁺ 通量,Ca²⁺ 也能触发 Ca²⁺ 诱导的 Ca²⁺ 释放(CICR),通过 RyR 从肌浆网或内质网释放(Endo 等 1970; Fabiato 1983)。Ca²⁺ 也可在体外使 RyR 失活,但该失活是否在体内起任何重要生理作用、甚至是否显著发生,尚未明确。

RyR 因对 ryanodine 的敏感性而得名——ryanodine 降低通道开放概率;咖啡因则升高 RyR 开放概率。已有许多 RyR 模型,但因大多用于兴奋-收缩耦合模型,将在第 12.2.4 节讨论。本节讨论少数为其他细胞类型开发的模型之一。

Calcium Oscillations in Bullfrog Sympathetic Neurons. 交感神经元对咖啡因或轻度去极化呈现稳健且可重复的 Ca²⁺ 振荡。这些振荡依赖胞外 Ca²⁺,但在固定膜电位下发生,涉 RyR 介导的 ER Ca²⁺ 释放(由 ryanodine 消除振荡可证明)。典型振荡如图 7.14。

Friel 1995 的简单 CICR 模型对牛蛙交感神经元中这些振荡的行为提供了优秀的定量描述。尽管模型简单(或许正因简单),它是理论补充实验的典范——对实验结果给出解释并提供可被后续检验的定量预测。

模型示意如图 7.15A。单一胞内 Ca²⁺ 储库与胞质交换 Ca²⁺(通量 \(J_{L2}\)\(J_{P2}\)),胞质再与胞外介质交换(\(J_{L1}\)\(J_{P1}\))。故

\[ \frac{\mathrm{d}c}{\mathrm{d}t} = J_{L1} - J_{P1} + J_{L2} - J_{P2}, \quad \frac{\mathrm{d}c_e}{\mathrm{d}t} = \gamma(-J_{L2} + J_{P2}), \]

其中 \(c\) 为胞质 [Ca²⁺],\(c_e\) 为胞内储库 [Ca²⁺],\(\gamma\) 为胞质与 ER 体积比。通量以浓度的简单线性函数给出:

\[ J_{L1} = k_1 (c_o - c), \quad J_{P1} = k_2 c, \quad J_{L2} = k_3 (c_e - c), \quad J_{P2} = k_4 c, \]

其中 \(c_o\) 为固定胞外 [Ca²⁺](图 7.14 实验中 \(c_o\) 固定为 1、0.5、0.7 mM)。高胞外 K⁺ 引起的去极化被建模为 \(k_1\)(Ca²⁺ 从胞外进入的速率)的增加,咖啡因施加(增加 Ca²⁺ 从胞内储库释放的速率)被建模为 \(k_3\) 的增加。

CICR 以简单方式建模:使 \(k_3\) 成为 \(c\) 的递增 sigmoid 函数

\[ k_3 = \kappa_1 + \frac{\kappa_2 c^n}{K_d^n + c^n}, \]

非线性模型的参数由拟合振荡时程确定(表 7.5)。典型结果如图 7.15B。该模型不仅提供了 Ca²⁺ 振荡的优秀定量描述(使用与表 7.5 不同的参数集),还预测了在振荡周期中应观察到的通量;后续对这些通量的测量确认了模型预测。因此 CICR(至少在牛蛙交感神经元中)可由相对简单的模型良好描述——RyR 仅需被 Ca²⁺ 激活即可产生生理振荡,Ca²⁺ 的失活并非必需。

7.3 Calcium Waves

某些细胞类型中 Ca²⁺ 振荡在整个细胞上几乎均匀发生——任何位置的 [Ca²⁺] 时间进程相同,此时适用充分混合模型。但更常见的情形是每个振荡呈现为穿越细胞的波;这些胞内"振荡"实际上是周期性的胞内波。建模并理解这种空间分布行为必须包含 Ca²⁺ 扩散;研究 Ca²⁺ 螺旋波(图 7.3)等对象再次需要偏微分方程模型。

建模空间分布 Ca²⁺ 动力学时,最诱人的办法是简单给 (7.3) 加扩散项。然而稍加思考即知这不正确,原因有几条:大多数通量项表示 Ca²⁺ 跨边界(质膜、ER 或线粒体)的运动,唯一例外是与 buffer 结合/解离的通量;胞质空间高度不均匀,该空间内 Ca²⁺ 运动是否服从标准扩散过程并不清楚。更恰当的建模方式是将细胞视为由两个相互连通区域(胞质 \(\Omega_c\) 与 ER \(\Omega_e\))组成的三维结构(为简化忽略线粒体、高尔基体等其他子域)。在胞质与 ER 中 Ca²⁺ 经正常扩散并与 buffer 反应,故

\[ \frac{\partial c}{\partial t} = \nabla \cdot (D_c \nabla c) - J_{\text{on}} + J_{\text{off}}, \quad \text{在 } \Omega_c \text{ 内}, \]
\[ \frac{\partial c_e}{\partial t} = \nabla \cdot (D_e \nabla c) - J_{\text{on,e}} + J_{\text{off,e}}, \quad \text{在 } \Omega_e \text{ 内}. \]

跨质膜进入胞质的通量产生边界条件

\[ D_c \nabla c \cdot \mathbf{n} = J_{\text{in}} - J_{\text{pm}}, \quad \text{在 } \partial \Omega_{c,m} \text{ 上}, \]

跨 ER 进入胞质的通量产生

\[ D_c \nabla c \cdot \mathbf{n} = -D_e \nabla c_e \cdot \mathbf{n} = J_{\text{IPR}} + J_{\text{RyR}} - J_{\text{serca}}, \quad \text{在 } \partial \Omega_e \text{ 上}, \]

其中 \(\partial \Omega_{c,m}\) 为质膜,\(\partial \Omega_e\) 为 ER 边界,\(\mathbf{n}\) 为所讨论区域的单位外法向量。注意此处假设 ER 与胞外空间无直接通信(此假设与生物学中几乎所有假设一样并非无争议)。由于 ER 的复杂几何,这些模型显然太复杂而难以取得任何进展,必须作简化假设。这是均化方法发挥作用之处——其基本思想是扩散在短距离上快速,局部变化被迅速平滑,只需知道平均或均场行为。即可推导出描述较大空间尺度行为的均场方程。在这些描述中假设 \(c\)\(c_e\) 在空间每一点共存(参见附录 7.8):

\[ \frac{\partial c}{\partial t} = \nabla \cdot (D_c^{\text{eff}} \nabla c) + \chi_c (J_{\text{IPR}} + J_{\text{RyR}} - J_{\text{serca}}) - J_{\text{on}} + J_{\text{off}}, \]
\[ \frac{\partial c_e}{\partial t} = \nabla \cdot (D_e^{\text{eff}} \nabla c_e) - \chi_e (J_{\text{IPR}} + J_{\text{RyR}} - J_{\text{serca}}) - J_{\text{on,e}} + J_{\text{off,e}}, \]

其中 \(D_c^{\text{eff}}\)\(D_e^{\text{eff}}\) 分别为胞质空间与 ER 的有效扩散系数,\(\chi_c\)\(\chi_e\) 为这两个共在空间的表面积与体积比。通常假设胞质各向同性且均匀;对均质性的例外并不少见,但对各向同性的例外罕见。Ca²⁺ 在 ER 中如何扩散、ER 的曲折度在多大程度上决定 ER Ca²⁺ 的有效扩散系数尚不明确。典型假设为 Ca²⁺ 在 ER 中不扩散,或以受限的扩散系数扩散,\(D_e^{\text{eff}} \ll D_c^{\text{eff}}\)。下文省略 eff 上标。

需要表面积与体积比以使单位正确。事实上应注意通量项(如 \(J_{\text{IPR}}\))在全细胞模型中与在空间分布模型中含义不同——全细胞模型中通量单位必须是浓度/单位时间,而空间分布模型中通量单位必须是 mol/单位时间/单位表面积。乘以表面积与体积比可将单位转换为 mol/体积/单位时间。但通常将所有通量缩放为 mol/单位时间/单位胞质体积,此时 (7.57) 与 (7.57) 在缩放后变为

\[ \frac{\partial c}{\partial t} = \nabla \cdot (D_c \nabla c) + J_{\text{IPR}} + J_{\text{RyR}} - J_{\text{serca}} - J_{\text{on}} + J_{\text{off}}, \]
\[ \frac{\partial c_e}{\partial t} = \nabla \cdot (D_e \nabla c_e) - \gamma (J_{\text{IPR}} + J_{\text{RyR}} - J_{\text{serca}}) - J_{\text{on,e}} + J_{\text{off,e}}, \]

其中 \(\gamma\) 为胞质与 ER 体积比。

当细胞某一维度远小于 Ca²⁺ 扩散长度尺度时可作进一步近似。例如许多培养细胞相对较薄,细胞厚度方向的空间变化可忽略;但推导正确的模型方程需稍加注意(详见附录 7.8)。以一长细圆柱形细胞为例,其圆柱壁上有边界通量 \(J_{\text{in}}\)\(J_{\text{pm}}\)(单位 mol/表面积/时间)。若细胞半径小于 Ca²⁺ 扩散长度尺度,则 Ca²⁺ 在每个横截面均匀分布,圆柱壁上的通量作为源项纳入控制偏微分方程。此时可将细胞建模为一维对象。例如长度为 \(L\)、端点无 Ca²⁺ 通量的"一维"细胞,边界条件为 \(\partial c/\partial x = 0\)\(x = L\)\(x = 0\) 处,内部点满足

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} + J_{\text{IPR}} + \frac{\rho}{A}(J_{\text{in}} - J_{\text{pm}}) + J_{\text{RyR}} - J_{\text{serca}} - J_{\text{on}} + J_{\text{off}}, \]

其中 \(\rho\) 为细胞周长、\(A\) 为横截面积。注意原始边界通量(单位 mol/表面积/时间)已转换为适合一维模型的 mol/长度/时间通量。

研究此类反应-扩散模型的方法有多种,但(至少在 Ca²⁺ 波研究中)两种最常用:数值仿真或行波方程的分岔分析。

7.3.1 Simulation of Spiral Waves in Xenopus

Xenopus 卵母细胞中引发波的常用实验程序是在细胞内光解释放 IP₃ 团并观察后续 Ca²⁺ 活动(Lechleiter and Clapham 1992)。充足时间后 Ca²⁺ 波活动随 IP₃ 降解而消失,但短期内观察到的 Ca²⁺ 活动是 Ca²⁺ 扩散与 IP₃ 扩散的结果。另一技术是释放 IP₃S₃——IP₃ 的非水解类似物,对 IP₃ 受体有相似作用但不被细胞降解。充分时间后 IP₃S₃ 在细胞所有部分达到恒定浓度。

在区域中部释放 IP₃S₃ 团会引起释放点处大量 Ca²⁺ 释放。IP₃S₃ 随后在细胞内扩散并在此过程中释放 Ca²⁺。被释放的 Ca²⁺ 激活 IP₃ 受体可引起储库周期性 Ca²⁺ 释放;IP₃ 受体之间 Ca²⁺ 扩散稳定波,产生规则的周期行波。这些周期波是时间模型中振荡的空间对应物,源自底层振荡动力学。若稳态 [IP₃S₃] 处于适当范围(例如图 7.12 显示 [IP₃] 在某中间范围存在极限环)则整个细胞的胞质部分都处于振荡状态。由带振荡动力学的反应-扩散系统标准理论(Kopell and Howard 1973; Duffy 等 1980; Neu 1979; Murray 2002),这些 [IP₃S₃] 值下可存在周期波与螺旋波。当释放 IP₃ 而非 IP₃S₃ 时波活动只持续很短时间,与理论结果一致。当波前被打破时常形成 Ca²⁺ 螺旋波(图 7.3A)。这些结果都已被数值仿真复现;Atri 等 1993 是首批完成此项工作者,但大多数 Ca²⁺ 波模型呈现相同的定性行为。视初始条件而定,这些螺旋波可稳定或不稳定。不稳定情形下螺旋的分支可自交并引起螺旋破碎,此时出现无明显空间结构的复杂模式区(McKenzie and Sneyd 1998)。

Falcke 等 1999, 2000 对螺旋 Ca²⁺ 波稳定性给出了更详细理解——他们显示线粒体 Ca²⁺ 释放速率增加可显著改变观察到的波类型,极端情况下可使胞质变为双稳,增加一个具有高静息 [Ca²⁺] 的稳态。此时螺旋波变得不稳定,导致出现更复杂的时空模式。

7.3.2 Traveling Wave Equations and Bifurcation Analysis

研究 Ca²⁺ 波的第二种主要方法是通过行波方程,最适用于研究一维空间波。引入行波变量 \(\xi = x + s t\)\(s\) 为波速,参见第 6 章 6.2 节),可将 (7.61) 写为方程对

\[ c' = d, \quad D_c d' = s d - \mathcal{J}, \]

其中 \(\mathcal{J}\) 表示 (7.61) 右端所有通量,撇号表示对 \(\xi\) 求导。单个反应-扩散方程转化为两条常微分方程。一般这两条方程与 \(c_e\)\(p\) 及各受体状态的方程耦合。行波脉冲、行波前缘与周期波分别对应行波方程的同宿轨道、异宿轨道与极限环。通过研究行波方程的分岔,可深入了解模型中存在哪些类型的波及其参数值范围。但该方法不给出原始反应-扩散方程行波解的稳定性信息——稳定性要难确定得多。反应-扩散方程的数值仿真可着手解决稳定性问题。

以第 7.2.6 节模型为例简述方法。给该模型加 \(c\) 的扩散(但不扩散 \(c_e\)):

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} + (k_f P_o + J_{er})(c_e - c) - J_{\text{serca}} + J_{\text{in}} - J_{\text{pm}}, \]
\[ \frac{\partial c_e}{\partial t} = \gamma [J_{\text{serca}} - (k_f P_o + J_{er})(c_e - c)], \]

这两方程如前与六态 IPR 模型的五条方程耦合(第 7.2.5 节)。(\(J_{\text{in}}\)\(J_{\text{pm}}\) 缩放为单位细胞长度的通量。)用行波变量重写:

\[ c' = d, \quad D_c d' = s d - (k_f P_o + J_{er})(c_e - c) + J_{\text{serca}} - J_{\text{in}} + J_{\text{pm}}, \]
\[ s c_e' = \gamma [J_{\text{serca}} - (k_f P_o + J_{er})(c_e - c)], \]

与五条受体方程耦合。注意按表 7.4 中参数单位,速度 \(s\) 的自然单位为 \(\mu\text{m}\,\text{s}^{-1}\)

这些方程的双参数数值分岔分析(用 AUTO)结果如图 7.16,\(s\)\(p\) 为分岔参数。\(s \to \infty\) 时行为对应无扩散模型,由一般理论(Maginu 1985)可预期。因此对大 \(s\) 值行为反映图 7.12 所示——有两个 Hopf 分岔,中间段 \(p\) 值有一条周期解支。若在 \(s, p\) 平面上追踪这些 Hopf 分岔形成 U 形曲线(图 7.16A)。为找对应行波(脉冲)的同宿分岔,取该图在恒定 \(s = 10\,\mu\text{m}\,\text{s}^{-1}\) 处的截面,即对该 \(s\) 值画出 \(c\)\(p\) 的分岔图(同图 7.12)。

结果如图 7.16B。两个 Hopf 分岔点对应图 A 中 U 形曲线与虚线(\(s = 10\,\mu\text{m}\,\text{s}^{-1}\))的交点。从最右 Hopf 分岔出发的周期轨道支以一条同宿轨道结束。在 \(s, p\) 平面上追踪该同宿轨道得到图 A 插图所示的同宿轨 C 形曲线。更精确地说,只有大周期轨道(此例 10000 秒)能被追踪,因此图 A 中标记 HC 的分支实际上是大周期轨道的分支。但追踪周期 1000 秒或 30000 秒的轨道支给出几乎相同的数值结果,可较有把握地说在追踪一条同宿轨道。

该同宿轨道分支给出原始反应-扩散方程(即偏微分方程)的稳定行波解。直接数值仿真易验证:固定 \(p = 0.5\),在一维域一端施加大 Ca²⁺ 刺激,所形成的波以约 \(15\,\mu\text{m}\,\text{s}^{-1}\) 的速度传播(与图 7.16A 一致),形状与该 C 形分支上的同宿轨道相同。

同宿轨道分支预测速度约 \(10\)\(17\,\mu\text{m}\,\text{s}^{-1}\) 的行波,处于生理范围内。此外随 \(p\) 增加波速也增加,与实验观察一致。C 形同宿分支的两端在生理上重要——上支对应于 \(p\) 缓慢增加时观察到的波。在某点存在从单个行波脉冲到周期波的转变,但该转变如何发生尚未完全阐明(Champneys 等 2007)。图 7.16 所示分岔图基本结构——同宿轨 C 形分支加 Hopf 分岔 U 形分支——似乎是激发系统模型的通有结构。FitzHugh–Nagumo 方程具有相同基本结构,Hodgkin–Huxley 方程亦如此。所有 Ca²⁺ 波传播模型(至少在已研究此问题的模型中)在 \(s, p\) 平面也具有相同的 C-U 结构。

7.4 Calcium Buffering

所有细胞中 Ca²⁺ 都被重度缓冲,至少 99%(通常更多)的可利用 Ca²⁺ 与大型 Ca²⁺ 结合蛋白结合,人类基因组编码约 200 种此类蛋白(Carafoli 等 2001)。例如 calsequestrin 与 calreticulin 是 ER/SR 中的主要 Ca²⁺ buffer;胞质中 Ca²⁺ 与 calbindin、calretinin、parvalbumin 等结合。Ca²⁺ 泵与交换器以及质膜本身也是重要的 Ca²⁺ buffer。实质上,胞质中处于游离状态的 Ca²⁺ 离子在被某种物质结合前能做的事或能走的距离都很有限。

Ca²⁺ 缓冲的基本化学反应可表示为

\[ P + \text{Ca}^{2+} \underset{k_-}{\overset{k_+}{\rightleftharpoons}} B, \]

其中 \(P\) 为缓冲蛋白、\(B\) 为已缓冲 Ca²⁺。记已结合 Ca²⁺ 的 buffer 浓度为 \(b\),自由 Ca²⁺ 浓度为 \(c\),Ca²⁺ 缓冲的简单模型为

\[ \frac{\partial c}{\partial t} = D_c \nabla^2 c + f(c) + k_- b - k_+ c (b_t - b), \]
\[ \frac{\partial b}{\partial t} = D_b \nabla^2 b - k_- b + k_+ c (b_t - b), \]

其中 \(k_-\) 为 Ca²⁺ 从 buffer 释放的速率,\(k_+\) 为 Ca²⁺ 被 buffer 摄取的速率,\(b_t\) 为 buffer 总浓度,\(f(c)\) 表示所有涉及自由 Ca²⁺ 的其他反应(IP₃ 受体释放、泵重摄等)。

7.4.1 Fast Buffers or Excess Buffers

若 buffer 动力学快速,其对胞内 Ca²⁺ 动力学的影响可简单分析(第 2.2.5 节)。若 \(k_-\)\(k_+\) 远大于 Ca²⁺ 反应的时间常数,则令 \(b\) 处于准稳态

\[ k_- b - k_+ c (b_t - b) = 0, \]

\[ b = \frac{b_t c}{K + c}, \]

其中 \(K = k_-/k_+\)。将 (7.70) 与 (7.71) 相加得"慢"方程

\[ \frac{\partial}{\partial t}(c + b) = D_c \nabla^2 c + D_b \nabla^2 b + f(c), \]

用 (7.73) 消去 \(b\) 后变为

\[ \frac{\partial c}{\partial t} = \frac{1}{1 + \theta(c)} \left[ \nabla^2\left(D_c c + D_b \frac{b_t c}{K + c}\right) + f(c) \right] = \frac{D_c + D_b \theta(c)}{1 + \theta(c)} \nabla^2 c - \frac{2 D_b \theta(c)}{(K + c)(1 + \theta(c))} |\nabla c|^2 + \frac{f(c)}{1 + \theta(c)}, \]

其中

\[ \theta(c) = \frac{b_t K}{(K + c)^2}. \]

假设 \(b_t\) 为常数,不随空间或时间变化。非线性缓冲显著改变模型。具体而言 Ca²⁺ 服从非线性扩散-平流方程,平流由可移动 buffer 输运 Ca²⁺ 引起(Wagner and Keizer 1994)。有效扩散系数

\[ D_{\text{eff}} = \frac{D_c + D_b \theta(c)}{1 + \theta(c)} \]

\(D_c\)\(D_b\) 的凸线性组合,故位于二者之间。因 buffer 是大分子,\(D_{\text{eff}} < D_c\)。若 buffer 不可移动(\(D_b = 0\))则 (7.76) 回到反应-扩散方程。当 Ca²⁺ 梯度较小时非线性平流项可忽略(Irving 等 1990)。最后缓冲也影响非线性反应项 \(f(c)\) 的定性性质——\(f(c)\)\(1 + \theta(c)\) 除。这可能改变模型的许多性质,包括振荡行为与波传播的性质。

若 buffer 不仅快速而且低亲和,使 \(K \gg c\),则

\[ b = \frac{b_t c}{K}, \quad \theta = \frac{b_t}{K}, \]

为常数。故 \(D_{\text{eff}}\) 也为常数。通常假设 buffer 快速、不可移动且低亲和——这三个假设下可得最简单的 Ca²⁺ buffer 模型(不包括它们的最简单情形):

\[ \frac{\partial c}{\partial t} = \frac{K}{K + b_t}\,(D_c \nabla^2 c + f(c)), \]

其中扩散系数与通量都被常数因子 \(K/(K + b_t)\) 缩放;模型中的每个通量可被解释为有效通量,即对自由 Ca²⁺ 浓度变化有贡献的那部分通量。

1986 年 Neher 观察到若 buffer 大大过量则 \(b_t - b \approx b_t\),缓冲反应变为线性,即所谓过量缓冲近似:

\[ \frac{\partial c}{\partial t} = D_c \nabla^2 c + f(c) + k_- b - k_+ c b_t, \quad \frac{\partial b}{\partial t} = D_b \nabla^2 b - k_- b + k_+ c b_t. \]

现再假设 buffer 快速则恢复 (7.79) 与 (7.81)。换言之,(7.81) 给出的 buffer 简单处理可由两种方式得到——要么假设低亲和 buffer,要么假设 buffer 过量。直观上显然这两种近似为何得到相同结果:两种情形下 Ca²⁺ 的结合都几乎不改变未结合 buffer 的比例。

表 7.6 给出三种不同 buffer 的典型参数值。BAPTA 是快速高亲和 buffer,EGTA 是慢速高亲和 buffer,都常用作实验中的外源 buffer;Endog 是典型内源 buffer(典型 \(b_t = 100\,\mu\text{M}\))。

许多 Ca²⁺ 缓冲研究关注单开放通道 Ca²⁺ 通量引起的稳态 Ca²⁺ 分布(Naraghi and Neher 1997; Stern 1992; Smith 1996; Smith 等 1996, 2001; Falcke 2003b)。这些研究的动机是:通道打开时,靠近通道口的稳态 Ca²⁺ 分布在微秒内达到。Smith 等 2001 在该语境下推导了快速缓冲近似与过量缓冲近似作为原方程的渐近解。尽管 (7.76) 复杂,仍保持单方程的优势。但若 buffer 动力学相对 Ca²⁺ 动力学不快速,则唯一办法是对完整系统做数值仿真(Backx 等 1989; Sala and Hernández-Cruz 1990; Nowycky and Pinter 1993; Falcke 2003a)。

当前特别关注的是近期观察——动力学不同的 buffer 对观察到的 Ca²⁺ 响应有显著不同影响。例如 Dargan and Parker 2003 显示在卵母细胞中施加 EGTA(慢速高亲和 buffer)"巴尔干化"了 Ca²⁺ 响应,使局部区域倾向于独立于邻居响应。但施加 BAPTA(快速 buffer)有相反效果——给出较慢的全局性 Ca²⁺ 响应而非局部响应。这些结果可能源自这两种 buffer 对单个释放位点间 Ca²⁺ 扩散的不同影响(第 7.6 节)。胰腺腺泡细胞中类似研究(Kidd 等 1999)显示施加 EGTA 增加频率但降低 Ca²⁺ 尖峰幅度,而施加 BAPTA 有相反效果——降低频率但维持较大尖峰幅度。此外 EGTA 将响应分解为多个独立的、不同空间的释放事件(与卵母细胞中一样)。Kidd 等结论是 EGTA 破坏释放位点之间的远程 Ca²⁺ 扩散,从而降低 Ca²⁺ 释放的全局协调性;BAPTA 快得足以破坏 Ca²⁺ 对 IPR 的失活,从而导致大 Ca²⁺ 尖峰而频率较低。然而尚未有详细建模研究证实这些定性解释。

7.4.2 The Existence of Buffered Waves

因快速 Ca²⁺ buffer 改变了 Ca²⁺ 输运方程的性质,确定 Ca²⁺ 缓冲如何影响波的性质有重大意义。例如施加 buffer 是否能消除波活动?buffer 显著影响行波速度吗?施加外源 buffer(如荧光 Ca²⁺ 染料)是否影响 Ca²⁺ 波的存在或速度?Tsai 及其同事(Tsai and Sneyd 2005, 2007a,b; Guo and Tsai 2006)在这些问题上做了大量工作,结果可简单总结:不可移动 buffer 对缓冲双稳方程中行波的存在性、稳定性与唯一性无影响;可移动 buffer 在数量足够时可消除行波。但当波存在于可移动 buffer 存在下时,它们仍保持唯一与稳定。

这些结果的证明过于技术性,此处不予讨论。Keener 改为讨论 Sneyd 等 1998 早期工作中较简单的若干结果。首先讨论 buffer 是否能消除波活动。

(7.75) 的形式提示变量替换

\[ w = D_c c + D_b \frac{b_t c}{K + c}, \]

由于

\[ \frac{\mathrm{d}w}{\mathrm{d}c} = D_c + D_b \theta(c) > 0, \]

\(w\)\(c\) 的单调递增函数。其唯一逆函数记为 \(c = \phi(w)\)。关于 \(w\),(7.75) 变为

\[ \frac{\partial w}{\partial t} = \left(D_c + \frac{D_b}{1 + \theta}\right) \nabla^2 w + f(\phi(w)), \]

其中 \(\theta = b_t K / (K + \phi(w))^2\)。现假设 \(f(c)\) 为双稳形式,有三个零点 \(C_1 < C_2 < C_3\),其中 \(C_1\)\(C_3\) 稳定。则 \(f(\phi(w))\) 立即有对应三个零点 \(W_1 < W_2 < W_3\)\(W_1\)\(W_3\) 稳定。(7.87) 一维行波解的存在性证明使用与第 6 章中双稳方程完全相同的论证(Sneyd 等 1998)。故提供 \(W_1\)\(W_3\) 转移的行波解存在当且仅当

\[ \int_{W_1}^{W_3} f(\phi(w))\, \mathrm{d}w > 0. \]

(若该不等式反向,行波解仍存在但反向移动,提供 \(W_3\)\(W_1\) 的转移。)用 (7.85) 将该条件改写为关于 \(c\) 的形式:

\[ \int_{C_1}^{C_3} f(c) (D_c + D_b \theta(c))\, \mathrm{d}c > 0. \]

一般而言该积分不能显式计算。但对简单情形双稳动力学 \(f(c) = c(1-c)(c-a)\)\(0 < a < 1/2\),显式计算 (7.89) 积分表明行波存在当且仅当

\[ a < a_c = \frac{1}{2} \cdot \frac{D_c - 12 D_b b_t K \left[\left(3K^2 + 2K\right) \ln\left(\frac{K+1}{K}\right) - \left(3K + \frac{1}{2}\right)\right]}{D_c + 12 D_b b_t K \left[\left(K + \frac{1}{2}\right) \ln\left(\frac{K+1}{K}\right) - 1\right]}. \]

从 (7.89) 可立即得到一个结论:静止 buffer(即 \(D_b = 0\))对双稳方程中行波的存在性无影响——因 \(D_b = 0\) 时行波存在条件 (7.89) 化为

\[ \int_0^1 f(c)\, \mathrm{d}c > 0, \]

这恰是无 buffer 时波存在条件。注意 \(a_c\)\(D_b b_t / D_c\) 的单调递减函数,且

\[ a_c \to a_{c,\min}(K) = -\frac{1}{2} \cdot \frac{\left[(3K^2 + 2K) \ln\left(\frac{K+1}{K}\right) - (3K + \frac{1}{2})\right]}{\left[(K + \frac{1}{2}) \ln\left(\frac{K+1}{K}\right) - 1\right]} \]

\(D_b \to \infty\)。图 7.17 给出 \(a_{c,\min}\)\(K\) 的曲线。\(K\) 大时 \(a_c\) 最小值接近 0.5,故波存在性对 \(D_b\) 不敏感;但 \(K\) 小时 \(a_c\)\(D_b\) 增加也变小,故此时可移动 buffer 易使波停止。

7.5 Discrete Calcium Sources

迄今所有模型均假设 ER 的 Ca²⁺ 释放在空间上均匀。事实并非如此。例如 Xenopus 卵母细胞中 IPR 以约 1/30 \(\mu\text{m}^2\) 的密度成簇分布,每簇约含 25 个 IPR。此外 Ca²⁺ 波传播是跳跃式的,释放在簇间跳跃。为研究这种跳跃式波的性质及其离散结构效应,假设 Ca²⁺ 从离散释放位点释放但被连续移除:

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} - k_s c + \sum_n \delta(x - nL) f(c), \]

其中 \(f\) 是任何描述 Ca²⁺ 释放的函数,\(L\) 是释放位点之间的空间间距。两个有用的释放函数示例

\[ f(c) = A c^2 (c_e - c), \quad f(c) = A H(c - c^*) (c_e - c), \]

其中 \(H\) 为通常的 Heaviside 函数。两者均有 Ca²⁺ 依赖的释放,表示 CICR。

对应连续空间模型可由均化方法得到(练习 16)。不难证明(且直觉上合理)当 \(L^2 k_s / D_c \ll 1\) 时,空间非均匀问题

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} - k_s c + g(x) f(c), \]

其中 \(g(x)\)\(L\) 为周期,可替换为其平均

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} - k_s c + G f(c), \]

其中

\[ G = \frac{1}{L} \int_0^L g(x)\, \mathrm{d}x \]

。因有效释放函数 \(F(c) = G f(c) - k_s c\) 为双稳型,有三个零点 \(0 < c_1 < c_2\),可预期存在 Ca²⁺ 释放行波,只要介质足够激发,即

\[ \int_0^{c_2} F(c)\, \mathrm{d}c > 0. \]

但若 \(L^2 k_s / D_c \neq 0\) 该条件不正确。

若存在驻波解则波不能传播。驻波是 (7.93) 的定常解,即

\[ 0 = D_c \frac{\partial^2 c}{\partial x^2} - k_s c + \sum_n \delta(x - nL) f(c) \]

的解。在区间 \(nL < x < (n+1)L\) 上变为

\[ 0 = D_c \frac{\partial^2 c}{\partial x^2} - k_s c. \]

通过从 \(nL^-\)\(nL^+\) 积分得 \(x = nL\) 处的跳跃条件

\[ D_c c_x\big|_{nL^+}^{nL^+} + L f(c_n) = 0, \]

其中 \(c_n = c(nL)\)。解 (7.100) 得

\[ c(x) = \frac{(c_{n+1} - c_n \cosh\beta) \sinh\left(\frac{\beta}{L}(x - nL)\right)}{\sinh\beta} + c_n \cosh\left(\frac{\beta}{L}(x - nL)\right), \]

\(nL < x < (n+1)L\),其中 \(c_n = c(nL)\)\(\beta^2 = k_s L^2 / D_c\),故

\[ c_x(nL^+) = \frac{(c_{n+1} - c_n \cosh\beta)}{\beta / L}\, \frac{1}{\sinh\beta}. \]

类似地

\[ c_x(nL^-) = -\frac{(c_{n-1} - c_n \cosh\beta)}{\beta / L}\, \frac{1}{\sinh\beta}. \]

故 (7.101) 即差分方程

\[ \frac{k_s}{\beta \sinh\beta}\,(c_{n+1} - 2 c_n \cosh\beta + c_{n-1}) + f(c_n) = 0, \]

这是 \(c_n\) 的差分方程。

求非线性差分方程的解一般非平凡。但若 \(f(c)\) 是分段线性的(如释放函数 (7.95)),可找到解析解。寻求 (7.105) 形如

\[ c_n = \begin{cases} a \mu_0^{-n}, & n \le 0, \\ C - b \mu_f^n, & n > 0, \end{cases} \]

的解,其中稳态解 \(C\) 满足

\[ \frac{k_s}{\beta \sinh\beta}\,(C - 2C \cosh\beta + C) + A(c_e - C) = 0, \]

\[ C = \frac{c_e}{A_f / (2\cosh\beta - 2) + 1}, \quad A_f = \frac{A \beta \sinh\beta}{k_s}. \]

\(\mu_0\)\(\mu_f\) 都小于 1,满足二次方程

\[ \mu_j - 2\lambda_j + \frac{1}{\mu_j} = 0, \quad j = 0, f, \]

其中

\[ \lambda_0 = \cosh\beta, \quad \lambda_f = \cosh\beta + \frac{A_f}{2}, \]

\[ \mu_j = \lambda_j - \sqrt{\lambda_j^2 - 1}. \]

易知 \(\mu_0 = e^{-\beta}\)\(\mu_f\) 的表达式更复杂。

通过考察 \(n = 0\)\(n = 1\) 处的差分方程确定标量 \(a\)\(b\)。因假设 \(c_0 < c^*\)\(c_1 > c^*\),故 \(f(c_0) = 0\)\(f(c_1) = A(c_e - C + b\mu_f)\)。经一些代数运算得

\[ a = C \left(\frac{\mu_f - 1}{\mu_f - 1/\mu_0}\right), \quad b = C \left(\frac{1 - 1/\mu_0}{\mu_f - 1/\mu_0}\right). \]

这些驻波存在条件为

\[ a \le c^*, \quad C - b\mu_f \ge c^*, \]

故画出 \(a\)\(C - b\mu_f\) 有揭示意义。图 7.18 给出 \(a/c_e\)(下曲线)与 \(C/c_e - (b/c_e) \mu_f\)(上曲线)对 \(\beta = \sqrt{k_s L^2 / D}\) 的曲线,对 \(A/k_s = 5\)。该图的解释是——固定 \(\beta\) 值下,若 \(c^*/c_e\) 位于曲线 \(a/c_e\) 之下则有行波;若 \(c^*/c_e\)\(\beta\) 位于两条曲线之间则存在驻波,排除传播可能。故一般而言 \(\beta\) 越大,释放位点必须越激发(\(c^*/c_e\) 必须越小)才能发生传播。换言之,从成簇受体进行的离散释放使传播可能性低于相同 Ca²⁺ 量在空间连续均匀释放的情形。注意对固定 \(c^*/c_e\)\(\beta\) 充分增大总会导致传播失败。增大 \(\beta\) 对应增大释放位点间距 \(L\) 或摄取速率 \(k_s\),或减小 Ca²⁺ 扩散系数 \(D_c\)

Sneyd and Sherratt 1997 用类似方法探讨了此问题。

7.5.1 Fire-Diffuse-Fire Model

确定何时发生传播失败已很困难,在离散释放模型中找到传播速度则更困难。Keener 2000b 用了一种超出本书范围的方法。另一种方法是使用不同的模型。引起广泛关注的模型之一是 fire-diffuse-fire 模型(Pearson and Ponce-Dawson 1998; Keizer 等 1998; Ponce-Dawson 等 1999; Coombes 2001; Coombes and Bressloff 2003; Coombes and Timofeeva 2003; Coombes 等 2004)。

在该模型中一旦 [Ca²⁺] 达到释放位点的阈值 \(c^*\),该位点立即触发,瞬间释放固定量 \(\sigma\) 的 Ca²⁺。故 Ca²⁺ 波通过释放位点的顺序触发传播,每个释放位点响应于邻近位点扩散来的 Ca²⁺——故名 fire-diffuse-fire。

假设 Ca²⁺ 服从反应-扩散方程

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} + \sigma \sum_n \delta(x - nL) \delta(t - t_n), \]

其中 \(L\) 同前为释放位点间距。虽然该方程看似线性,但表象具有欺骗性。\(t_n\)\(c\) 首次在第 \(n\) 个释放位点达到阈值 \(c^*\) 的时间,发生时第 \(n\) 个位点释放量 \(\sigma\)。故 \(t_n\) 以复杂方式依赖 \(c\)

单个位点(如位点 \(i\))触发产生的 Ca²⁺ 分布为

\[ c_i(x, t) = \frac{\sigma H(t - t_i)}{\sqrt{4\pi D_c (t - t_i)}} \exp\left(-\frac{(x - iL)^2}{4 D_c (t - t_i)}\right), \]

其中 \(H\) 为 Heaviside 函数。这是扩散方程在 \(x = i, t = t_i\)\(\delta\) 函数输入的基本解,可在任何关于偏微分方程解析解的标准书(参见 Keener 1998 或 Kevorkian 2000)中找到。叠加各位点的解

\[ c(x, t) = \sum_i c_i(x, t) = \sigma \sum_i \frac{H(t - t_i)}{\sqrt{4\pi D_c (t - t_i)}} \exp\left(-\frac{(x - iL)^2}{4 D_c (t - t_i)}\right). \]

注意因瞬时释放,\(c(x, t)\) 在任何释放位点上不是时间的连续函数。

现假设位点 \(i = N, N-1, \ldots\) 已在已知时刻 \(t_N > t_{N-1} > \cdots\) 触发。下一触发时刻 \(t_{N+1}\)\(x_{N+1}\)\(c\) 首次达到阈值 \(c^*\) 的时刻确定,即

\[ c((N+1)L, t_{N+1}^-) = c^*, \quad \frac{\partial}{\partial t} c((N+1)L, t_{N+1}^-) > 0. \]

\(t_{N+1}\) 必须满足

\[ c^* = \sigma \sum_{i \le N} \frac{1}{\sqrt{4\pi D_c (t_{N+1} - t_i)}} \exp\left(-\frac{L^2 (N+1-i)^2}{4 D_c (t_{N+1} - t_i)}\right). \]

稳定传播波对应 \(t_i - t_{i-1} = \text{常数} = \tau\) 对所有 \(i\) 成立,即每位点在其左侧邻居触发后固定时间 \(\tau\) 触发。注意所得波不以恒定分布传播,但有明确波速 \(L/\tau\)。若这样的 \(\tau\) 存在则 \(t_{N+1} - t_i = \tau(N+1-i)\)\(\tau\) 是方程的解

\[ \frac{c^* L}{\sigma} = \sum_{n=1}^\infty \frac{1}{\sqrt{4\pi n \eta}} \exp\left(-\frac{n}{4\eta}\right) \equiv g(\eta), \]

其中 \(\eta = D_c \tau / L^2\) 为无量纲延迟。

为求 \(\eta\) 需反演该方程。\(g(\eta)\) 曲线如图 7.19A。可证明 \(0 \le g(\eta) \le 1\)\(g\) 单调,\(g \to 0\)\(\eta \to 0\)\(g \to 1\)\(\eta \to \infty\)。故 (7.121) 的解存在仅当 \(c^* L / \sigma < 1\)。即簇间距离或阈值过大或释放量过小时无传播。但 \(c^* L / \sigma < 1\) 时 (7.121) 的唯一解必有,故存在传播波。

通过适当反转图 7.19A 的坐标轴易画出延迟作为 \(c^* L / \sigma\) 的函数。故图 7.19B 中画出无量纲延迟 \(\eta\) 作为 \(c^* L / \sigma\) 的函数。同样易画出无量纲速度 \(1/\eta\) 作为 \(c^* L / \sigma\) 的函数(未画出)。结果是 \(c^* L / \sigma \to 0\) 时速度无穷大,\(c^* L / \sigma \ge 1\) 时速度为零。

若 (7.121) 的解存在则波速与 \(D_c\) 成正比——因速度 \(L/\tau = D / (\eta L)\)。该结果令人不安有两个原因。第一,空间均匀反应-扩散系统中波的传播速度通常与 \(D\) 的平方根成正比。由简单尺度论证可知——因扩散系数 \(D_c\) 单位为(长度)²/时间,距离变量可由 \(\sqrt{D_c k}\)\(k\) 为典型时间常数)缩放以去除所有距离单位,故波速与 \(\sqrt{D_c k}\) 成正比。第二,这与上一节显示 \(D_c\) 足够小时存在传播失败的结果矛盾。对 fire-diffuse-fire 模型传播成功或失败与 \(D_c\) 无关。

此不匹配可由 fire-diffuse-fire 模型仅允许 Ca²⁺ 释放但不允许摄取来解释——故 Ca²⁺ 瞬态不现实地单调增加。该 fire-diffuse-fire 模型的缺陷易通过加线性空间均匀摄取项补救(Coombes 2001),故模型变为

\[ \frac{\partial c}{\partial t} = D_c \frac{\partial^2 c}{\partial x^2} - k_s c + \sigma \sum_n \delta(x - nL) \delta(t - t_n). \]

此修正模型的分析几乎与 \(k_s = 0\) 情形相同。基本解由 \(k_s\) 略作修正——位点 \(i\) 触发产生的 Ca²⁺ 分布为

\[ c_i(x, t) = \frac{\sigma H(t - t_i)}{\sqrt{4\pi D_c (t - t_i)}} \exp\left(-\frac{(x - iL)^2}{4 D_c (t - t_i)} - k_s(t - t_i)\right). \]

按前述论证得知存在传播解的条件是 (7.124) 方程

\[ \frac{c^* L}{\sigma} = \sum_{n=1}^\infty \frac{1}{\sqrt{4\pi n \eta}} \exp\left(-\frac{n}{4\eta} - \beta^2 n\right) \equiv g_\beta(\eta) \]

有解,其中 \(\eta = D_c \tau / L^2\) 为无量纲延迟,\(\beta^2 = k_s L^2 / D_c\)

\(\beta^2\) 大小的估计可变化很大。例如用 \(k_s = 143/\text{s}\)(表 7.4 中 \(1/\alpha_2\))、\(L^2 = 30\,\mu\text{m}^2\)\(D_c = 20\,\mu\text{m}^2/\text{s}\),得 \(\beta^2 = 172\);而用 \(k_s = 3.78/\text{s}\)(Friel 模型表 7.5 中 \(k_4\))、\(L^2 = 4\,\mu\text{m}^2\)(适合心肌细胞 Ca²⁺ 释放)、\(D_c = 25\,\mu\text{m}^2/\text{s}\),得 \(\beta^2 = 0.6\)。不论如何 \(\beta\) 的效应显著。\(g_\beta(\eta)\) 对若干 \(\beta\) 值的曲线如图 7.20。特别地若 \(\beta \neq 0\) 则函数 \(g_\beta(\eta)\) 非单调递增,而有最大值 \(g_{\max}(\beta)\),是 \(\beta\) 的递减函数。进一步 \(g_\beta(\eta) \to 0\)\(\eta \to 0\)\(\eta \to \infty\)。若 \(c^* L / \sigma > g_{\max}(\beta)\) 则 (7.124) 无解,传播失败。另一方面若 \(c^* L / \sigma < g_{\max}(\beta)\) 则 (7.124) 有两个解;物理上有效的解是其中较小者,对应 \(c(x, t)\) 首次达到 \(c^*\) 的时刻。

对较大 \(\beta\) 值,\(g_{\max}(\beta)\) 与指数函数 \(\exp(-\beta)\) 吻合良好。这给出传播失败的近似判据——\(c^* L / \sigma > \exp(-\beta)\) 则传播失败。回顾 \(\beta = \sqrt{k_s L^2 / D_c}\) 可见该判据对 \(D_c\) 的依赖。

通过反转图 7.20 的坐标轴,可得无量纲延迟作为 \(c^* L / \sigma\) 的函数。注意 \(\beta \neq 0\)\(\beta = 0\) 时显著定性差异。\(\beta \neq 0\) 时传播在有限(而非无穷)延迟处停止。这暗示传播在正(而非零)速度处失败。

7.6 Calcium Puffs and Stochastic Modeling

上述所有模型中 Ca²⁺ 释放被建模为确定性的。但目前已知这并非总是适当。事实上每次 Ca²⁺ 振荡或波由若干随机基本释放事件(称 puff)构成,每个对应来自单个或一小群 IPR 的 Ca²⁺ 释放。低 IP₃ 浓度下发生来自单簇的点状释放;高 [IP₃] 下这些局部释放事件被协调成全局胞内波(Yao 等 1995; Parker 等 1996a,b; Parker and Yao 1996)。puff 详细研究主要在 Xenopus 卵母细胞(Marchant 等 1999; Sun 等 1998; Callamaras 等 1998; Marchant and Parker 2001)与 HeLa 细胞(Thomas 等 2000; Bootman 等 1997a,b)中进行。

图 7.21 给出一些典型实验结果。上图显示对单激光脉冲的 Ca²⁺ 释放响应(细胞先负载可光释放形式的 Ca²⁺)。Ca²⁺ 脉冲是局部的,不扩展形成全局波——因背景 IP₃ 浓度过低无法支持波。但通过光释放 IP₃ 提高背景 IP₃ 浓度后(图 B),胞质变得更具激发性,发生自发 Ca²⁺ puff(例如图顶部中心的白点)。若自发释放碰巧足够大则形成传播波,从释放位点向其他位点扩展,如图右下角所示。图 C 显示较高背景 IP₃ 浓度下施加 Ca²⁺ 脉冲的响应。此时介质足够激发使 Ca²⁺ 脉冲能引发全局传播波。

这些实验结果提出几个重要建模问题。第一,何时适用确定性模型、何时必须纳入随机行为?第二,如何最佳建模通过少量 IPR 的随机 Ca²⁺ 释放?第三,如何建模这种局部释放向胞内波等全局事件的协调?尽管这些问题从卵母细胞与 HeLa 细胞实验工作中最为明显,对其他细胞类型的 Ca²⁺ 振荡研究也重要。如 Falcke 2004 所指出,随机效应看起来如此基本且广泛,以至于对确定性方法的一般适用性提出质疑。

7.6.1 Stochastic IPR Models

建模 Ca²⁺ puff 的基本假设是每个 IPR 释放事件可建模为随机 Markov 过程,而 Ca²⁺ 扩散与其他泵与 Ca²⁺ 通量的作用可确定性建模。给定 IPR 的 Markov 状态模型(例如图 7.8 或图 7.10),可在每时间步选一个随机数,用该随机数决定该时间步的状态变化。这最高效地使用 Gillespie 方法(第 2.9.3 节)完成。

一般每个跃迁速率是 \(c\)\(p\) 的函数,故跃迁概率随浓度变化而不断变化。因此 \(c\)\(p\) 必须在每时间步更新——这种更新是确定性的,通过解 \(c\)\(p\) 的反应-扩散方程完成。当模拟的 IPR 模型处于开放态时,Ca²⁺ 的反应-扩散方程有附加的 IPR 通量。当 IPR 处于其他状态时该通量缺失。故得到一个由随机变化输入驱动的 \(c\) 反应-扩散方程。

第一个 IPR 随机模型由 Swillens 等 1998 提出。其模型是 18 态模型——IPR 可结合 0 或 1 个 IP₃、0、1 或 2 个激活 Ca²⁺、0、1 或 2 个失活 Ca²⁺。仅当 IPR 结合一个 IP₃、两个 Ca²⁺ 于激活位点、无 Ca²⁺ 于失活位点时受体开放。模型稳态开放概率受 Bezprozvanny 等 1991 实验数据约束,并假设 Ca²⁺ 沿径向从通道口扩散。Ca²⁺ 在通道口可累积到高浓度,这些局部浓度用于 IPR 模型的随机仿真。

该模型的仿真显示两件事。第一,通道开放以突发形式发生——因 Ca²⁺ 从通道口扩散得足够慢,允许重新结合激活位点。第二,通过与实验观察的 puff 幅度分布比较,Swillens 等 1999 显示典型簇约含 25 个受体,簇内 IPR 间隔可能不超过 12 nm。

因对 IPR 门控方案进行直接随机仿真计算量大,Shuai and Jung 用 Langevin 方程近似八态模型(第 7.2.5 节)的随机版本(Shuai and Jung 2002a,b, 2003),该方法已由 Fox and Lu 1994 用于 Hodgkin–Huxley 方程。虽然每簇 25 个受体使 Langevin 方程方法不如直接随机仿真精确,但定性行为与直接仿真吻合良好。

Local Concentrations. 任何 IPR 随机模型中不可能忽略通道口发生的局部高浓度。忽略该因素要么得到与现实无关的 IPR 模型,要么得到不呈现现实 Ca²⁺ 振荡的模型。故对要纳入全细胞 Ca²⁺ 振荡模型的 IPR 随机模型,必须设法将通道口微域的 Ca²⁺ 浓度与胞质中总体 Ca²⁺ 浓度联系起来。Huertas and Smith 2007 与 Bentele and Falcke 2007 提出了两种方法。

该问题与如何在心肌细胞中最佳建模 Ca²⁺ 释放的问题相似(第 12.2.4 节)。那里 RyR 将 Ca²⁺ 释放到非常受限的域——diadic cleft——并经历远高于胞质其余部分的 Ca²⁺ 浓度。故许多近期心肌细胞 Ca²⁺ 动力学模型不使用单一胞质 Ca²⁺ 域,而纳入各种 Ca²⁺ 微域,复杂度各异。

7.6.2 Stochastic Models of Calcium Waves

所有随机建模结果中最引人入胜的是 Falcke 2003a,b 的工作——研究 Xenopus 卵母细胞中 puff 向波的转变随 [IP₃] 增加。低 [IP₃] 下仅观察到 puff;从每簇释放的 Ca²⁺ 不足以刺激邻近簇释放 Ca²⁺,故响应纯粹局部。但 [IP₃] 增加时每个 IPR 的敏感性以及释放的 Ca²⁺ 量都增加。这允许从成核位点发展出全局波。但在 [IP₃] 大大增大前,这些全局事件罕见,许多情形下仅形成传播短距离后即死的夭折波。波间时间间隔 \(T_{\text{av}}\) 及其标准差 \(\Delta T_{\text{av}}\) 随 [IP₃] 增加而下降。最后高 [IP₃] 下全局波规则发生,周期明确。

Falcke 2003a,b 显示所有这些行为可由第 7.2.5 节八态模型的随机版本复现。低 [IP₃] 下连续波间的长时间间隔几乎完全是随机动力学的产物。每个时间间隔内有机会一个 IPR 触发,刺激整簇触发从而启动全局波——但因簇间间隔与每个 IPR 的低敏感性这些事件罕见。故此时 \(T_{\text{av}}\)\(\Delta T_{\text{av}}\) 都大。反之 [IP₃] 大时每个 IPR 对 Ca²⁺ 更敏感,每个受体的通量更大。故单个 IPR 触发几乎总足以刺激全局波。此时波间周期不由随机效应设定,而由 IPR 的内禀动力学设定——受体重新激活并准备传播另一波所需时间。这些结果在图 7.22 中最易看到——图 7.22 给出 Marchant and Parker 2001 的实验数据与 Falcke 2003a 的对应仿真。\(\Delta T_{\text{av}}/T_{\text{av}}\) 比值在很宽 \(T_{\text{av}}\) 范围内保持近似常数。

特别有趣的是这些结果显示在仿真中周期全局波发生于模型确定性版本为非振荡的 IP₃ 浓度下。即振荡波未必是振荡动力学的产物。它们可源于随机过程——该过程不时引起某簇触发足够强从而启动全局波。若平均触发时间间隔的标准差足够小,所得准周期响应可看似底层极限环的产物,即使不存在这样的极限环。这些结果对由少量随机 IPR 引起的 Ca²⁺ 波与振荡的确定性方法的适用性提出严肃质疑。然而其含义尚未被充分理解(Keener 2006)。

7.7 Intercellular Calcium Waves

Ca²⁺ 波不仅在单细胞内传播,也从细胞到细胞传播,形成可穿越许多细胞的胞间波。最早的胞间波例子之一由 Sanderson 等 1990, 1994 发现——上皮细胞培养中机械刺激(例如用微吸管戳单细胞)可启动胞内 Ca²⁺ 增加的波,从细胞到细胞传播形成胞间波。气道上皮细胞的典型实验结果如图 7.4。上皮细胞培养形成薄层细胞,由间隙连接相连。当培养中间某细胞被机械刺激时,被刺激细胞中 Ca²⁺ 迅速增加。一两秒延迟后,被刺激细胞的邻居也显示 Ca²⁺ 增加,该增加在培养中依次传播。胞内波穿过每个细胞,在细胞边界延迟,然后启动相邻细胞的胞内波。胞间波通过胞内波的顺序传播而移动。特别有趣的是——在无胞外 Ca²⁺ 时被刺激细胞无响应,但胞间波仍传播到培养中其他细胞。故被刺激细胞中 Ca²⁺ 上升似乎对波传播并非必需。Ca²⁺ 上升也非引发胞间波的充分条件——例如培养中上皮细胞有时呈现自发胞内 Ca²⁺ 振荡,这些振荡并不从细胞到细胞传播。但机械刺激的胞间波确实传播过自发振荡的细胞。

Charles 等 1991 与 Cornell-Bell 等 1990 也研究了胶质细胞培养中的胞间 Ca²⁺ 波。过去几年累积的越来越多证据表明这种胶质细胞间或胶质细胞与神经元间的胞间通信在脑信息处理中起重要作用(Nedergaard 1994; Charles 1998; Vesce 等 1999; Fields and Stevens-Graham 2002; Lin and Bergles 2004)。

正如不同细胞类型中存在多种胞间 Ca²⁺ 波,其传播机制也存在相应的多样性。但两种基本机制占主导:通过胞外信使扩散传播、通过胞内信使经间隙连接扩散传播。有时两种机制共同驱动胞间波(参见 Young and Hession 1997)。最常见胞内信使为 IP₃ 或 Ca²⁺(或两者),但更多胞外信使已被涉及,包括 ATP、ADP、一氧化氮。

胞间 Ca²⁺ 波模型较少。最早的由 Sneyd 等 1994, 1995a, 1998 研究机械诱导波的机制,Young 1997 提出另一早期模型。该基本模型的更新版本已用于研究肝细胞(Höfer 1999; Höfer 等 2001, 2002; Dupont 等 2000)与胰腺腺泡细胞(Tsaneva-Atanasova 等 2005)的胞间耦合;Jung 等 1998 与 Ullah 等 2006 采用不同方法。胞内与胞外信使之间相互作用的研究更少(Bennett 等 2005 与 Iacobas 等 2006 是两个近期例子),关于这类模型如何行为还有许多待发现。

7.7.1 Mechanically Stimulated Intercellular Ca2+ Waves

Sanderson 及其同事(Boitano 等 1992; Sanderson 等 1994; Sneyd 等 1994, 1995a,b)提出了上皮细胞中机械刺激胞间 Ca²⁺ 波的模型(图 7.23)。他们提出机械刺激引起被刺激细胞中产生大量 IP₃,该 IP₃ 通过被动扩散在培养中移动,经间隙连接从细胞到细胞。因 IP₃ 从 ER 释放 Ca²⁺,IP₃ 从细胞到细胞的扩散导致相应胞间 Ca²⁺ 波。实验结果表明细胞间 Ca²⁺ 移动对波传播不起主要作用,故模型假设细胞间 Ca²⁺ 移动可忽略。放松该假设对模型行为影响不大——因 IP₃ 经间隙连接的移动决定胞间波性质。

模型中上皮细胞培养被建模为正方形细胞网格。假设 IP₃ 经被动扩散移动且以饱和动力学降解。故若 \(p\) 为 [IP₃]:

\[ \frac{\partial p}{\partial t} = D_p \nabla^2 p - V_p \frac{p}{k_p + p}. \]

\(p \ll k_p\)\(p\) 以时间常数 \(1/V_p\) 衰减。Ca²⁺ 也假设经被动扩散,但由 IP₃ 从 ER 释放并由 Ca²⁺ ATPase 泵回 ER。方程为

\[ \frac{\partial c}{\partial t} = D_c \nabla^2 c + J_{\text{IPR}} - J_{\text{serca}} + J_{\text{in}}, \]
\[ \tau_h \frac{\mathrm{d}h}{\mathrm{d}t} = \frac{k_2^2}{k_2^2 + c^2} - h, \]
\[ J_{\text{IPR}} = k_f \mu(p) h \frac{b + (1-b) c}{k_1 + c}, \quad J_{\text{serca}} = \frac{\gamma c^2}{k_{2\gamma} + c^2}, \quad J_{\text{in}} = \beta, \quad \mu(p) = \frac{p^3}{k_\mu^3 + p^3}. \]

这是 Ca²⁺ 动力学的早期模型,反映在其相对简单性。IPR 模型基于 Atri 等 1993;\(J_{\text{IPR}}\)\(p\)\(c\) 与慢变量 \(h\)(未因 Ca²⁺ 失活的 IPR 比例)的函数。IP₃ 引起的开放被假设为瞬时,Ca²⁺ 激活亦瞬时;故项 \(\mu(p) [(b + (1-b) c)/(k_1 + c)]\)\(p\)\(c\) 的递增函数。失活变量 \(h\) 的稳态是 Ca²⁺ 的递减函数,以时间常数 \(\tau_h\) 达到该稳态。故 Ca²⁺ 快速激活后跟缓慢失活的性质被纳入模型方程。如常 \(J_{\text{serca}}\) 表示 ER 膜上 Ca²⁺ ATPase 对胞质 Ca²⁺ 的移除,以系数 2 的 Hill 方程建模(基于 Lytton 等 1992 数据);\(J_{\text{in}}\) 为不特指的 Ca²⁺ 漏入胞质,来自胞外或来自 ER。模型参数值见表 7.7。

该模型与本章早前讨论的模型有一重要区别——此处无描述 ER 中 Ca²⁺ 浓度的变量,即无 \(c_e\) 变量。这等价于假设 ER 耗竭可忽略,故 \(c_e\) 为常数。

内部边界条件以 IP₃ 从细胞到细胞的通量给出。若细胞 \(n\) 有 [IP₃] \(= p_n\),假设 IP₃ 从细胞 \(n\) 到细胞 \(n+1\) 的通量为 \(F(p_n - p_{n+1})\)\(F\) 为常数,称渗透率)。最初将 IP₃ 注入单细胞,允许其从细胞到细胞扩散,由此生成胞间 Ca²⁺ 波。图 7.24 给出模型方程在二维中数值解的密度图。可看到胞间波在细胞网格上扩展并随后随 IP₃ 降解而退去。如对基于被动扩散的过程所预期,胞内波速(即胞间波穿过单细胞的速度)随距被刺激细胞的距离而下降,到达时间与胞间延迟随距被刺激细胞的距离指数增加。对所选 \(F\) 值范围 1 至 8 \(\mu\text{m}\,\text{s}^{-1}\),模型与上皮、内皮与胶质细胞的实验数据吻合良好(Demer 等 1993; Charles 等 1992)。但模型最重要的预测是获得这种吻合所需的 \(F\) 值。若 \(F\) 低于约 \(1\,\mu\text{m}\,\text{s}^{-1}\),胞间波传播太慢与实验数据不符。因 \(F\) 值未知,该预测提供检验 IP₃ 被动扩散底层假说的方法。

7.7.2 Partial Regeneration

前述模型提出的主要问题之一是 IP₃ 分子能否在不降解的情况下扩散穿越多个细胞。一个 IP₃ 分子不太可能存活足够久以引起传播距离达几百微米的胞间波。但若 IP₃ 的生成是可再生的——例如通过 PLC 被增加的 [Ca²⁺] 激活——则波原则上可无限传播,如同神经轴突上的动作电位(第 6 章)。因多数细胞类型波最终停止,这一事实初看排除了 IP₃ 生成的可再生机制。

该问题由 Höfer 等 2002 解决——他们显示部分可再生机制可使胞间 Ca²⁺ 波传播比简单扩散可能距离更远,但非无限;波最终仍停止。其模型用于研究大鼠纹状体星形胶质细胞培养中的远程胞间 Ca²⁺ 波——该准备中已知胞间波主要通过胞内信使的间隙连接扩散传播。模型方程与前一节类似,但有两处主要例外。第一,模型包含 \(c_e\) 的动力学;第二,模型假设 IP₃ 由某 PLC 亚型(PLCδ)生成的速率是 [Ca²⁺] 的递增函数。故

\[ \frac{\partial p}{\partial t} = D_p \nabla^2 p + \frac{\nu_7 c^2}{K_{\text{Ca}}^2 + c^2} - k_9 p, \]

其中 \(K_{\text{Ca}} = 0.3\,\mu\text{M}\)\(k_9 = 0.08\,\text{s}^{-1}\)\(\nu_7\) 在 0 与 0.08 间变化。

结果如图 7.25。当 Ca²⁺ 对 IP₃ 生成的正反馈小时(即 PLCδ 最大活性 \(\nu_7\) 小时),胞间 Ca²⁺ 波行为如同源自被刺激细胞 IP₃ 被动扩散,与机械刺激胞间 Ca²⁺ 波模型相同。反之当 \(\nu_7\) 大时波变为完全可再生,无限传播。然而数值计算显示存在 \(\nu_7\) 值的中间范围给出有限再生波——波传播距离超过若仅由 IP₃ 被动扩散控制时所能达到的距离,但最终停止,因此非完全可再生。

7.7.3 Coordinated Oscillations in Hepatocytes

周期性胞间 Ca²⁺ 波传播跨越整个肝小叶(Nathanson 等 1995; Robb-Gaspers and Thomas 1995),该行为在较小的簇(双倍体与三倍体)耦合肝细胞中也被反映。波传播机制似乎与机械刺激波提出的机制显著不同。第一,肝细胞中波仅在每个细胞被刺激时(即每个细胞有增加 IP₃ 时)传播。第二,波是周期性的。尽管有这些差异,已知波的协调依赖胞内信使的间隙连接扩散;无间隙连接耦合时不出现协调波活动,每个细胞独立于邻居振荡(Tordjmann 等 1997)。也已知肝小叶中存在激素受体密度梯度,从门静脉周围区向中央静脉周围区递增。该梯度在三倍体肝细胞中也被反映——用 noradrenaline 刺激时,三倍体中细胞按顺序响应,受体密度最高的细胞首先响应。因三倍体中每个细胞对胞内光释放 IP₃ 同步响应,故不同延迟源自不同 IP₃ 生成速率。

基于这些观察,Dupont 等 2000 提出肝细胞簇中协调波为相位波——每个细胞以略不同于邻居的相位振荡。在该模型中,胞间同步由 IP₃ 经间隙连接的扩散提供。因每个细胞以略不同速率生成 IP₃,胞间 IP₃ 扩散起到平滑胞间 IP₃ 梯度的作用,允许振荡更大程度的同步并出现相位波。这种同步不会持续很久——因细胞本质上解耦,模型中无驱动 Ca²⁺ 振荡的 IP₃ 振荡,故无胞间 Ca²⁺ 扩散时细胞在更长尺度上表现为解耦振荡器。然而同步持续时间足以解释实验观察。Höfer 1999 提出略不同的模型——细胞由胞间 Ca²⁺ 扩散耦合——也由 Dupont 及其同事检验。按后者的观点,基于胞间 Ca²⁺ 扩散的模型不足以解释所有实验数据。

胰腺腺泡细胞中胞间波的类似研究也已完成(Tsaneva-Atanasova 等 2005)——结论是胞间 Ca²⁺ 扩散单独能在长时间周期上同步胞间相位波,而胞间 IP₃ 扩散的功能是使胞间 IP₃ 梯度最小化,从而让 Ca²⁺ 扩散更容易同步振荡器。

然而这些结果仍存在相当大争议。肝细胞或胰腺腺泡细胞中 IP₃ 振荡对驱动 Ca²⁺ 振荡有多重要尚未明确。要阐明这些胞间波与耦合振荡的详细机制仍有大量工作要做。

7.8 Appendix: Mean Field Equations

第 7.3 节中为研究空间分布系统中的 Ca²⁺ 动力学,给模型方程加了扩散。但这样做使用了某些重要近似技术,本附录加以描述。

7.8.1 Microdomains

假设物质 \(u\) 在三维区域中扩散与反应,其中一两个维度与扩散长度尺度相比很小。(若所有维度都比典型扩散长度尺度小,合理假设区域充分混合,全细胞模型适用。)该区域可以是心肌细胞中连接 SR 与质膜之间的区域(diadic cleft)、两片平面 ER 之间的窄胞质区、或 ER 与线粒体之间的窄胞质区。这些微域在心肌细胞与骨骼肌中已知非常重要,在 Cajal 间质细胞、平滑肌、神经元等细胞中也被认为重要(尽管细节仍不确定)。

在此域中

\[ \frac{\partial u}{\partial t} = D \nabla^2 u + f, \]

边界条件为 \(\mathbf{n} \cdot D \nabla u = J\) 在域边界上,\(\mathbf{n}\) 为单位外法向量。为便于讨论假设 (7.133) 采用无量纲时间单位,故 \(D\) 单位为长度²。

现假设例如域为长圆柱管,横截面积 \(A\)\(D\) 小。因 \(A\)\(D\) 小,预期解在每个横截面上近乎均匀,仅与其平均值略有变化。故寻求描述平均值演化的方程——亦称均场方程。

为利用长度尺度差异,将 Laplacian 算子拆分为两部分:

\[ \nabla^2 u = u_{xx} + \nabla_y^2 u, \]

其中 \(x\) 为沿管长的坐标,\(y\) 表示管横截面坐标。管侧面边界条件也用横截面坐标系表示为 \(\mathbf{n} \cdot D \nabla_y u = J\)。此处忽略管端的边界条件。

均场方程的快速非形式推导是——定义 \(\bar{u}\)\(u\) 的横截面平均值

\[ \bar{u} = \frac{1}{A} \int_\Sigma u\, \mathrm{d}A, \]

其中 \(\Sigma\) 为横截面域,然后对 (7.133) 关于 \(\Sigma\) 积分得

\[ \frac{\partial \bar{u}}{\partial t} = D \bar{u}_{xx} + \frac{D}{A} \int_\Sigma \nabla_y^2 u\, \mathrm{d}A + \frac{1}{A} \int_\Sigma f\, \mathrm{d}A. \]

应用散度定理

\[ \int_\Sigma \nabla_y^2 u\, \mathrm{d}A = \int_{\partial\Sigma} \mathbf{n} \cdot \nabla_y u\, \mathrm{d}S \]

\[ \frac{\partial \bar{u}}{\partial t} = D \bar{u}_{xx} + \frac{S}{A} \bar{J} + \bar{f}, \]

其中

\[ \bar{J} = \frac{1}{S} \int_{\partial\Sigma} J\, \mathrm{d}S, \]

\(S\) 为横截面周长。(7.138) 即控制 \(\bar{u}\) 行为的方程。

虽然 (7.138) 给出正确答案,我们中某些人可能希望有更系统的推导。此外若 \(J\)\(f\) 依赖 \(u\),则不清楚如何确定 \(\bar{J}\)\(\bar{f}\)。为此引入小参数 \(\epsilon\),其中 \(\epsilon^2 = A/D\),引入缩放坐标 \(y = \sqrt{A}\, \xi\),(7.133) 变为

\[ \frac{\partial u}{\partial t} = D u_{xx} + \frac{1}{\epsilon^2} \nabla_\xi^2 u + f, \]

边界条件为

\[ \mathbf{n} \cdot \nabla_\xi u = \epsilon J / \sqrt{D} \]

。现寻求 (7.140) 的形如

\[ u = \bar{u} + \epsilon w \]

的解,其中 \(w\) 在每个横截面中要求平均值为零(即

\[ \int_{\Sigma_\xi} w\, \mathrm{d}A_\xi = 0 \]

),\(\bar{u}\)\(\xi\) 无关。重要之处在于 \(\bar{u}\)\(w\) 不独立于 \(\epsilon\),这不是解的幂级数表示,而是由投影算子

\[ P u = \int_{\Sigma_\xi} u\, \mathrm{d}A_\xi \]

对解的分解;故 \(\bar{u} = P u\)\(\epsilon w = u - P u\)

现将 (7.141) 代入 (7.140) 然后应用投影算子 \(P\)\(I - P\),求 \(\bar{u}\)\(w\) 的方程。首先应用 \(P\)(即对 \(\Sigma_\xi\) 积分)得

\[ \frac{\partial \bar{u}}{\partial t} = D \bar{u}_{xx} + \frac{1}{\epsilon} \frac{S_\xi}{\sqrt{D}} \bar{J} + \bar{f}, \]

其中

\[ \bar{J} = \frac{1}{S_\xi} \int_{\partial\Sigma_\xi} J\, \mathrm{d}S_\xi \]

\(S_\xi\) 为以 \(\xi\) 为单位度量的管周长。接下来应用算子 \(I - P\)

\[ \frac{1}{\epsilon} \nabla_\xi^2 w = \frac{1}{\epsilon} \frac{S_\xi}{\sqrt{D}} \bar{J} + \bar{f} - f + \epsilon \frac{\partial w}{\partial t} - \epsilon D w_{xx}. \]

函数 \(w\) 必须满足沿管侧的边界条件

\[ \mathbf{n} \cdot \nabla_\xi w = J / \sqrt{D} \]

\(J\)\(f\) 均为已知函数,不依赖 \(u\),则完成——因 (7.143) 即精确的均场方程,不再需要近似。

然而若 \(J\)\(f\) 依赖 \(u\)(典型情形),则需更多了解 \(w\)。特别地若 \(w\) 有界且为一阶量,则

\[ \bar{J} = \frac{1}{S_\xi} \int_{\partial\Sigma_\xi} J(u)\, \mathrm{d}S_\xi = \frac{1}{S_\xi} \int_{\partial\Sigma_\xi} \left[J(\bar{u}) + \epsilon J_u(\bar{u}) w + O(\epsilon^2)\right] \mathrm{d}S_\xi = J(\bar{u}) + O(\epsilon^2), \]

类似地 \(\bar{f} = f(\bar{u}) + O(\epsilon^2)\)。故均场方程为

\[ \frac{\partial \bar{u}}{\partial t} = D \bar{u}_{xx} + \frac{1}{\epsilon} \frac{S_\xi}{\sqrt{D}} J(\bar{u}) + f(\bar{u}) + O(\epsilon), \]

或用原始参数表示

\[ \frac{\partial \bar{u}}{\partial t} = D \bar{u}_{xx} + \frac{S}{A} J(\bar{u}) + f(\bar{u}) + O\!\left(\sqrt{\frac{A}{D}}\right). \]

仍需确认 \(w\) 有界且为一阶量。为此寻求 (7.144) 形如幂级数

\[ w = w_1 + \epsilon w_2 + O(\epsilon^2) \]

\(w\) 解(固定 \(\bar{u}\),即忽略 \(\bar{u}\)\(\epsilon\) 的隐式依赖)。显然 \(w_1\) 必须满足

\[ \nabla_\xi^2 w_1 = \frac{S_\xi}{\sqrt{D}} \bar{J}, \]

边界条件为

\[ \mathbf{n} \cdot \nabla_\xi w_1 = J(\bar{u})/\sqrt{D} \]

。现令 \(W(\xi)\) 为边界值问题

\[ \nabla_\xi^2 W = S_\xi \text{ 在 } \Sigma \text{ 上}, \]
\[ \mathbf{n} \cdot \nabla_\xi W = 1 \text{ 在 } \partial\Sigma \text{ 上}, \quad \int_\Sigma W(\xi)\, \mathrm{d}A_\xi = 0 \]

的基本解(由 Poisson 方程标准理论可知其存在)。故

\[ w_1 = \frac{J(\bar{u})}{\sqrt{D}} W(\xi). \]

这至少可令我们满意地确认 \(w\)\(\epsilon\) 的良性有界函数,证毕。

假设感兴趣的域不是长细圆柱,而是位于距离 \(L\) 隔开的两片平面二维膜之间,其中 \(L^2 \ll D\)。并假设 \(D u_z = -J_0\)\(D u_z = J_1\) 分别在下膜与上膜上,其中 \(z\) 表示垂直空间坐标。该情形下可用同样方法(练习 19)证明均场方程为

\[ \frac{\partial \bar{u}}{\partial t} = D \nabla^2 \bar{u} + \frac{1}{L}\bigl(J_1(\bar{u}) + J_0(\bar{u})\bigr) + f(\bar{u}) + O\!\left(\frac{L}{\sqrt{D}}\right), \]

其中

\[ \bar{u} = \frac{1}{L} \int_0^L u\, \mathrm{d}z \]

\(\nabla^2\) 表示二维 Laplacian 算子。

7.8.2 Homogenization; Effective Diffusion Coefficients

如第 7.3 节所述,要将胞质与 ER 边界微观结构的细节纳入考虑既不切实际也无大用。相反,需要用有效扩散系数找到 Ca²⁺ 浓度的均场描述。避免微观结构细节的需要在许多其他语境下明显。例如如第 12 章所述,需要不依赖细胞结构及其连通性的心肌组织动作电位传播方程。

均化是实现这一目标非常强大的技术。本节展示如何使用均化找到带有效扩散系数的平均或均场方程。同样的技术将在第 12 章中调用以找到心肌组织的有效电导。该技术也允许推导 Ca²⁺ 的双域方程与心肌动作电位的双域方程。

作为热身问题,假设物质在一维区域中反应并扩散,扩散系数在空间中快速变化。具体说假设 \(u\) 由反应-扩散方程

\[ \frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\!\left[D\!\left(\frac{x}{\epsilon}\right) \frac{\partial u}{\partial x}\right] + f(u) \]

控制。这里 \(x\) 无量纲,\(D(x)\) 是以 1 为周期的一阶量函数,\(\epsilon\) 小。预期 \(u\) 应有某种平均或均场行为,特征长度尺度为一阶量,有相对该均场的一阶小量级小变化。

为探讨此可能性,引入两个变量

\[ z = x, \quad \xi = \frac{x}{\epsilon}, \]

视为独立变量。由链式法则

\[ \frac{\partial}{\partial x} = \frac{\partial}{\partial z} + \frac{1}{\epsilon} \frac{\partial}{\partial \xi}, \]

原始偏微分方程 (7.153) 变为

\[ \frac{\partial u}{\partial t} = \frac{\partial}{\partial z}\!\left[D(\xi)\!\left(\frac{\partial u}{\partial z} + \frac{1}{\epsilon} \frac{\partial u}{\partial \xi}\right)\right] + \frac{1}{\epsilon} \frac{\partial}{\partial \xi}\!\left[D(\xi)\!\left(\frac{\partial u}{\partial z} + \frac{\partial u}{\partial \xi}\right)\right] + f(u). \]

虽然该方程明显比 (7.134) 更复杂,其结构本质上相同。故后续计算与前一节相同。即寻求 (7.156) 形如

\[ u = \bar{u} + \epsilon w \]

的解,其中 \(\bar{u}\)\(\xi\) 无关,\(w\)\(\xi\) 的周期函数,平均值为零,

\[ \int_0^1 w\, \mathrm{d}\xi = 0 \]

。但相对于应用投影算子然后求这些的幂级数解,直接对 (7.156) 求幂级数解稍更方便(本情形等价)。即设 \(\bar{u} = u_0 + \epsilon u_1 + O(\epsilon^2)\)\(w = w_1 + \epsilon w_2 + O(\epsilon^2)\),代入 (7.156),收集 \(\epsilon\) 同次幂项,导出层次方程

\[ \frac{\partial}{\partial \xi}\!\left[D(\xi)\!\left(\frac{\partial w_1}{\partial \xi} + \frac{\partial u_0}{\partial z}\right)\right] = 0, \]
\[ \frac{\partial}{\partial \xi}\!\left[D(\xi) \frac{\partial w_2}{\partial \xi}\right] = \frac{\partial u_0}{\partial t} - \frac{\partial}{\partial z}\!\left[D(\xi)\!\left(\frac{\partial u_0}{\partial z} + \frac{\partial w_1}{\partial \xi}\right)\right] - \frac{\partial}{\partial \xi}\!\left[D(\xi)\!\left(\frac{\partial w_1}{\partial z} + \frac{\partial u_1}{\partial z}\right)\right] - f(u_0). \]

第一个方程由直接积分易解。定义 \(W(\xi)\) 为具有零均值的周期函数,满足微分方程

\[ \frac{\mathrm{d}W}{\mathrm{d}\xi} = \frac{R}{\bar{R}} - 1, \]

其中 \(R = 1/D\)

\[ \bar{R} = \int_0^1 R\, \mathrm{d}\xi \]

。则

\[ w_1 = -W(\xi) \frac{\partial u_0}{\partial z}. \]

函数 \(W(\xi)\) 决定解的小尺度结构。

接下来考察 (7.159) 并观察 \(w_2\) 能为周期函数仅当该方程右端对 \(\xi\) 有零平均。故必有

\[ \frac{\partial u_0}{\partial t} = \int_0^1 \frac{\partial}{\partial z}\!\left[D(\xi)\!\left(\frac{\partial u_0}{\partial z} + \frac{\partial w_1}{\partial \xi}\right)\right] \mathrm{d}\xi + f(u_0) = D_{\text{eff}} \frac{\partial^2 u_0}{\partial z^2} + f(u_0), \]

其中

\[ D_{\text{eff}} = \frac{1}{\bar{R}}. \]

(7.162) 即所求的均场方程,\(D_{\text{eff}}\) 为有效扩散系数。

完整三维问题的均化遵循同样的步骤。假设物质在三维空间的子区域中反应并扩散,该子区域细分为小周期单元 \(\Sigma\)(称为微观结构),每个包含在小长方盒中。\(\Sigma\) 又细分为胞质区 \(\Sigma_c\) 与 ER 区 \(\Sigma_e\)(图 7.26A)。\(\Gamma_m\) 为 ER 与胞质之间的膜边界,\(\Gamma_b\) 为 ER 与盒子壁的交。因假设微域周期重复,一盒的 ER 通过 \(\Gamma_b\) 与相邻盒的 ER 连接。

假设长方盒长度为 \(l\)(虽然盒子不必为立方体),其中 \(l\) 远小于问题的自然长度尺度——此情形下为 \(\lambda = \sqrt{D}\)(扩散长度)。故存在两个自然长度尺度 \(l\)\(\lambda\),故 \(\epsilon = l/\lambda \ll 1\) 是自然的小参数。因 \(l\) 与扩散长度尺度相比较小,合理假设解在每个微观结构单元中近乎均匀,微观结构只引起对背景解的小扰动。

现假设物质 \(u\)\(\Sigma_c\) 中按 (7.133) 扩散与反应,越过 \(\Gamma_m\) 的通量为 \(\mathbf{n} \cdot D \nabla u = J\)。另一变量(如 \(v\))将在 \(\Sigma_e\) 中反应与扩散,但对两个变量进行完整分析不必要。再次方便地假设这些方程采用无量纲时间单位,故 \(D\) 单位为长度²。进一步假设 \(J\)\(l\) 成比例。注意越过 \(\Gamma_m\) 的总通量与表面积成正比,故每单位体积的总通量与 \(J/l\) 成正比。为使 \(\epsilon \to 0\) 时该值有界,假设 \(J = l j\)

标准过程是引入两个变量

\[ z = \frac{x}{\lambda}, \quad \xi = \frac{z}{\epsilon}, \]

\(\xi\) 表示微观尺度上的空间变量,\(z\) 表示原始或长尺度上的空间变量。将 \(z\)\(\xi\) 视为独立变量,写

\[ \nabla_x = \frac{1}{\lambda}\!\left(\frac{1}{\epsilon} \nabla_\xi + \nabla_z\right), \]

下标表示梯度导数所取的空间变量。在这些变量下 (7.133) 变为

\[ \frac{\partial u}{\partial t} = \nabla_z^2 u + \frac{2}{\epsilon} \nabla_\xi \nabla_z u + \frac{1}{\epsilon^2} \nabla_\xi^2 u + f, \]

边界条件

\[ \mathbf{n} \cdot \!\left(\nabla_z u + \frac{1}{\epsilon} \nabla_\xi u\right) = \epsilon j \text{ 在 } \Gamma_m \text{ 上}. \]

注意 \(\mathbf{n}\) 为无量纲量,仅含边界法向的方向信息,故不随空间尺度变化而变化。

再次得到结构与 (7.134) 相同的偏微分方程。故与前一样,寻求 (7.166) 形如

\[ u = \bar{u}(z, t, \epsilon) + \epsilon w(\xi, z, t, \epsilon) \]

的解,其中 \(\bar{u}\) 为背景或均场解,仅在大空间尺度上变化,与 \(\xi\) 无关;而 \(\epsilon w\)\(\xi\) 上的周期校正项,在 \(\xi\) 变量上具有零平均值。下一步寻求 \(\bar{u}\)\(w\) 的幂级数表示

\[ \bar{u} = u_0 + \epsilon u_1 + O(\epsilon^2), \quad w = w_1 + \epsilon w_2 + O(\epsilon^2). \]

代入 (7.166) 与 (7.167) 并令 \(\epsilon\) 同次幂系数相等,给出可顺序求解的层次方程。领导阶方程为

\[ \nabla_\xi^2 w_1 = 0 \text{ 在 } \Sigma_c \text{ 内}, \]

边界条件

\[ \mathbf{n} \cdot (\nabla_z u_0 + \nabla_\xi w_1) = 0 \text{ 在 } \Gamma_m \text{ 上}, \]

下一个层次方程为

\[ \frac{\partial u_0}{\partial t} = \nabla_z^2 u_0 + 2 \nabla_\xi \nabla_z w_1 + \nabla_\xi^2 w_2 + f \text{ 在 } \Sigma_c \text{ 内}, \]

边界条件

\[ \mathbf{n} \cdot (\nabla_z w_1 + \nabla_\xi w_2) = j \text{ 在 } \Gamma_m \text{ 上}. \]

第一个方程可由 Bensoussan 等 1978 方法求解——令

\[ w_1 = \mathbf{W}(\xi) \cdot \nabla_z u_0, \]

其中 \(\mathbf{W}(\xi)\) 为满足向量微分方程

\[ \nabla_\xi^2 \mathbf{W} = 0 \text{ 在 } \Sigma_c \text{ 内}, \quad \mathbf{n} \cdot (\nabla_\xi \mathbf{W} + \mathbf{I}) = 0 \text{ 在 } \Gamma_m \text{ 上} \]

的基本解。为唯一确定 \(\mathbf{W}(\xi)\),还要求

\[ \int_{\Sigma_c} \mathbf{W}(\xi)\, \mathrm{d}V_\xi = 0. \]

\(u_0\) 的均场方程可通过在 \(\Sigma_c\) 上积分 (7.172) 得到

\[ V_c \frac{\partial u_0}{\partial t} = \int_{\Sigma_c} \nabla_z \cdot (\nabla_z u_0 + \nabla_\xi w_1)\, \mathrm{d}V_\xi + \int_{\Sigma_c} \nabla_\xi \cdot (\nabla_z w_1 + \nabla_\xi w_2)\, \mathrm{d}V_\xi + \int_{\Sigma_c} f\, \mathrm{d}V_\xi \]
\[ = \int_{\Sigma_c} \nabla_z \cdot (\nabla_z u_0 + \nabla_\xi \mathbf{W} \nabla_z u_0)\, \mathrm{d}V_\xi + \int_{\partial\Sigma_c} \mathbf{n} \cdot (\nabla_\xi w_2 + \nabla_z w_1)\, \mathrm{d}S_\xi + \int_{\Sigma_c} f\, \mathrm{d}V_\xi \]
\[ = \nabla_z \cdot \left[\int_{\Sigma_c} (\mathbf{I} + \nabla_\xi \mathbf{W}) \nabla_z u_0\, \mathrm{d}V_\xi\right] + \int_{\Gamma_m} j(u_0)\, \mathrm{d}S_\xi + V_\Sigma f(u_0). \]

注意 \(\partial\Sigma_c\) 的积分约化为 \(\Gamma_m\) 上的积分,因要求 \(w_1\)\(w_2\)\(\xi\) 上为周期函数。

(7.178) 是相容性条件,为 \(w_2\) 存在所需。换言之,(7.172) 是 \(w_2\) 的微分方程,若 \(w_2\) 为周期的,仅当相容性或可解性条件成立时才可解(这是 Fredholm 选择定理的应用)。所需相容性条件通过对 \(\Sigma_c\) 积分 (7.172) 得 (7.178)。

故得

\[ \frac{\partial u_0}{\partial t} = \nabla_z \cdot D_{\text{eff}} \nabla_z u_0 + \frac{S_c}{V_c} j(u_0) + f(u_0), \]

其中 \(V_c\)\(S_c\) 分别表示 \(\Sigma_c\) 的(无量纲)体积与表面积,\(D_{\text{eff}}\)(无量纲)有效扩散张量为

\[ D_{\text{eff}} = \frac{1}{V_c} \int_{\Sigma_c} (\mathbf{I} + \nabla_\xi \mathbf{W})\, \mathrm{d}V_\xi. \]

故对任何给定的周期几何,有效扩散张量 (7.180) 可通过解 (7.175) 与 (7.176) 计算。Goel 等 2006 计算了多种可能 ER 微观结构下的有效扩散系数。对图 7.26A 的规则结构,胞质中 Ca²⁺ 的有效扩散系数随 ER 体积分数从 0 增至 0.9 下降 60%。图 7.26B 更复杂结构给出非各向同性扩散——垂直方向的扩散系数比另两个方向大约 50%。

7.8.3 Bidomain Equations

此处描述的均化技术可应用于任何具有周期微观结构的空间域。故例如可应用于相互缠绕的胞内与胞外空间,或相互缠绕的胞质与 ER 空间。每种情形下可找到定义于空间各处的均场描述——故两个空间有效共在,故视为双域。唯一(轻微)复杂之处在于跨共同边界的通量必须在幅度上相等、方向相反。故进入胞内空间的通量即离开胞外空间的通量,进入胞质的通量即离开 ER 的通量。此外尽管域共享共同界面,它们的体积互补,故两个空间的表面积与体积比不同。

本章个人批注

本章体量非常大——25 个内容节,是 Keener 全书中节数最多的一章。读完后最强烈的感受是 Ca²⁺ 动力学的建模几乎没有任何"标准模型"可言,每个细胞类型、每个 IP₃ vs RyR、每个空间均匀 vs 空间不均匀情形都有专门的工具包组件选择,这正是为什么 Keener 一开始就强调"信号工具包"的模块化哲学。Falcke 2004 综述几乎 200 页似乎确实是必要的——光这一章就触及了从 IPR 八态 Markov 到 fire-diffuse-fire 跳跃波到微域均化的全光谱方法。

八态 IPR 模型的"对 Ca²⁺ 序贯激活/失活"作为振荡机制的核心假说很美,但 Keener 自己立刻补刀说"某些近期模型提示 IPR 动力学重要性不如先前所想,ER 耗竭可能也起重要作用"——这种"先讲主流、再讲怀疑"的话术是 Keener 全书一贯的诚实风格。同样地,DeYoung-Keizer 八态模型虽漂亮但被实验发现其 IPR 开放速率仅在一个数量级内变化(质量作用不行)所否定,于是 Sneyd-Dufour 引入饱和结合动力学——这是经典的"模型被实验逼着进化"的案例。

我比较欣赏的几个细节:(1) 开放细胞 vs 封闭细胞的 \(c_t\) 重构,把膜与 ER 通量在形式上分开,让我想到 Ch6 中关于 cAMP 的 open/closed cell 区分——本质是同一个数学操作。(2) 8.2 节均化的两层架构(先 microdomains 再 homogenization),非常清晰地把"小尺度结构 → 大尺度方程"的物理图像和数学对象对应起来。(3) Sneyd-Dufour 顺序模型中对"Ca²⁺ 介导激活与自发失活的内禀竞争"的描述,介于 IPR 门控的复杂性给我留下了具体印象。

比较困扰我的几点:(1) 第 7.4.1 节中有效扩散系数公式 \(D_{\text{eff}} = (D_c + D_b \theta)/(1+\theta)\) 的推导有一处 "\(-2 D_b \theta / ((K+c)(1+\theta)) |\nabla c|^2\)" 的非线性平流项——Keener 说它源自可移动 buffer 输运 Ca²⁺,但没给出直觉解释为何是非线性平流而非线性扩散。(2) 第 7.5 节离散释放位点的驻波分析中 \((7.114)\) 系数 \(a\)\(b\) 的公式我读了两遍才看清楚是几何序列解的极限匹配——这部分若能在每步加上"为什么这么分组"的物理直觉会更有帮助。(3) 第 7.6 节关于随机版本的八态模型在确定性版本为非振荡的 IP₃ 浓度下给出周期波——这段非常反直觉,我需要再读一遍 Falcke 2003a 才能确信理解对了。

Falcke 2004 综述确实需要找时间精读。

与上下章的衔接(一段话)

第 7 章紧接第 6 章(兴奋介质与波)之后:第 6 章建立了行进波、分岔分析、稳定性等反应-扩散系统的抽象数学工具,第 7 章则把这些工具应用到 Ca²⁺ 这一最具实验可见性的具体系统上——特别是 7.3.2 节明确说"在 s, p 平面上的 C-U 结构与 FitzHugh–Nagumo 方程和 Hodgkin–Huxley 方程相同,是激发系统的通有结构",这是第 6 章抽象理论的具体验证。第 7 章也大量回引前面的章节:第 1.4 节 Hill 方程(SERCA 泵)、第 2 章 Markov 模型与细致平衡(SERCA MacLennan 模型、八态 IPR 模型)、第 3.5.4 节亲和/效能(Sneyd-Dufour 饱和结合 IPR 模型)、第 5 章 Hodgkin–Huxley 方程(IPR 失活变量 \(h\) 与 HH 的类比、CICR 与 Na⁺ 通道的相似)、第 6 章波传播理论(第 7.3.2 节直接用 \(\xi = x + st\))。第 7 章向前则通向第 8 章(细胞迁移,Ca²⁺ 信号在迁移中的作用)、第 9 章(胰腺 β-细胞阵发振荡,引用了 Rinzel 1985 的封闭/开放细胞方法)、第 12 章(兴奋-收缩耦合,Ca²⁺ 在心肌细胞中的特殊作用,本章多次预告)。第 7 章本身几乎自足:它完整展示了从 IPR/RyR 分子机制到全细胞振荡、到空间行波、到随机 puff、再到胞间波与组织均化的全尺度建模链条,是 Keener 全书方法论的集大成演示。