钙离子通道和钙信号震荡模型研究

Modeling of Ca²⁺ Channels and Ca²⁺ Signal Oscillations

摘要: 钙离子是细胞内广泛存在的第二信使，控制着几乎所有的生命过程。本文以钙信号的生物背景为基础，首先介绍了构建钙离子信号网络模型的数学方法，包括生化反应建模、参数拟合等；接着我们介绍了关于1，4，5-三磷酸肌醇受体(IP₃R)通道和线粒体钙单向转运体(MCU)通道的各种模型；然后我们讨论了关于钙离子信号震荡的各种点源模型，包括内质网释放的钙信号震荡模型和内质网-线粒体相互作用调控的钙信号震荡模型；最后我们建议未来对钙信号相关的细胞信号网络建模展开进一步的研究。

Abstract: The calcium ions (Ca²⁺) are the second messengers which are ubiquitous in cells, regulating almost all the cellular processes. In this review, based on the biological background of Ca²⁺ signaling, we first introduce the mathematical methods for the construction of Ca²⁺ signaling models, including modeling of biochemical reactions and parameter fitting. Next, we discuss in detail the modeling of Ca²⁺ channels, including the 1, 4, 5-trisphosphate receptor channel (IP₃R) and mitochondrial Ca²⁺ signaling uniporter (MCU). Then we discuss various point source models of Ca²⁺ signal oscillations controlled by endoplasmic reticulum (ER) or the cross talk between ER and Mitochondria. Possible trends of modeling of Ca²⁺ signal in cellular signaling networks are also suggested.

1. 细胞钙离子信号

钙离子(Ca²⁺)是细胞内普遍存在的第二信使 [1] [2] ，它调节着细胞的多种生理过程，如心脏的活动、基因的转录、大脑信息的处理以及记忆的存储等 [3] [4] ，是细胞内重要的信号。

在静息状态下，细胞质内的钙离子浓度([Ca²⁺])维持在很低的水平，大约为100 nM以内 [2] [5] ，大部分的Ca²⁺主要存储于内质网(Endoplasmic Reticulum, ER)、肌浆网(Sarcoplasmic Reticulum, SR)和线粒体(mitochondria)中。内质网是真核细胞中重要的细胞器。1945年，K. R. Porter等人 [6] 首次在培养的细胞中观察到细胞质的内质部分有网状结构，因此将其命名为内质网。肌浆网主要存在于心肌细胞和骨骼肌细胞中，它是一种特殊的内质网。内质网是主要的钙库之一。它可以存储大部分的Ca²⁺，其浓度可以达到100~900 μM [7] ，有时甚至可以达到1 mM [8] 。因为内质网中含有一类Ca²⁺结合蛋白，这类结合蛋白对Ca²⁺有着高容量、低亲和力的特点。每个钙结合蛋白分子可与30个左右的Ca²⁺结合，大大提高了内质网存储Ca²⁺的能力。

线粒体也是存储Ca²⁺的容器之一。它与内质网、细胞外基质等协同合作，共同维持细胞内[Ca²⁺]的平衡。由于线粒体对Ca²⁺的亲和力不高 [9] ，因此静息状态下线粒体中的[Ca²⁺]不高，与细胞质内的[Ca²⁺]相差不大 [10] 。当细胞溶质内的[Ca²⁺]升高，为了维持细胞溶质内的钙平衡，一部分Ca²⁺又会迅速进入线粒体，之后又会释放出去，使线粒体成为钙离子的缓冲区。线粒体存储和释放的Ca²⁺可以改变细胞内的钙信号 [11] [12] 。

静息状态下，细胞内的钙离子浓度维持在较低水平，但当细胞受到刺激信号时，细胞外的Ca²⁺内流和内质网/肌浆网中Ca²⁺外流，使细胞质中的[Ca²⁺]迅速升高，形成不同的钙信号，从而执行其所需的细胞功能 [13] 。[Ca²⁺]的升高主要是由于“开反应”决定的 [4] [14] 。在“开反应”中，细胞外的Ca²⁺主要通过位于细胞膜上的门控通道进入细胞质中，这些门控通道主要包括：电压门控的钙通道(VOC)、受体门控的钙通道(ROC)和钙库控制的钙通道(SOC)；内质网中的Ca²⁺主要通过内质网膜上的IP₃R通道，雷诺定受体等Ca²⁺通道流入细胞质中 [4] 。

当钙信号执行完其生理功能，细胞质中的[Ca²⁺]通过“关反应”，使[Ca²⁺]降低至静息态水平 [5] 。在“关反应”中，细胞质内的Ca²⁺可通过细胞膜上的Ca²⁺泵(PMCA)和Na⁺-Ca²⁺交换体(NCX)流到细胞外，也可以通过内质网上的Ca²⁺泵重新进入内质网中 [15] [16] 。除此之外，线粒体上的Ca²⁺单向转运体(uniporter)可以快速结合Ca²⁺以促进钙恢复过程，且位于线粒体膜上的NCX可以使流入线粒体的Ca²⁺重新进入细胞质中 [4] 。

细胞膜上和内质网膜上都含有不同的Ca²⁺通道。在静息状态下，细胞内的Ca²⁺主要存储于内质网(ER)中。IP₃R通道是内质网上Ca²⁺的主要通道之一。IP₃R的大致结构可以分为三部分：位于细胞质中的N端，位于内质网膜附近的C端，以及两端的连接区域 [17] [18] 。IP₃R主要是由4个相同的亚基构成 [19] ，大概的分子量约为1200 kDa [14] [18] 。对于IP₃R通道的开关，主要是由1，4，5-三磷酸肌醇(IP₃)分子和Ca²⁺决定的。IP₃对IP₃R通道的调控是单向的，即IP₃与IP₃受体结合，可以改变受体的结构，增加其对Ca²⁺的敏感性。而Ca²⁺对IP₃R通道的调控是双向的，即在低浓度的Ca²⁺情况下，随着Ca²⁺浓度的升高，可以增加IP₃R对Ca²⁺敏感性及通道的开放概率，通过钙致钙释放机制促进Ca²⁺的进一步释放；而在高浓度的情况下，Ca²⁺会抑制通道的打开并促进通道的关闭以保证正常生理活动的进行 [20] - [23] 。

RYR通道主要是存在于肌浆网膜上，它也是Ca²⁺的主要通道之一。RYR通道的结构和IP₃R通道相似，都是四聚体，但是RYR通道的分子量较大。它的总分子量约有2 MD，每个亚基大约有550 kDa，是至今发现的肌质网上最大的离子通道 [24] 。对于Ca²⁺对RYR通路的影响，其实和IP₃R通路相似，即较低的[Ca²⁺]会促进通道的开放，而较高的[Ca²⁺]会导致通道的关闭。RYR通路与IP₃R通路一样，控制着Ca²⁺的释放，因而控制着大量的生理活动。

钙结合蛋白对钙信号传导有着重要的作用。在开反应中，Ca²⁺流入细胞内，钙结合蛋白主要以钙效应器和缓冲蛋白的形式和Ca²⁺相互作用 [25] 。在细胞内大多数的缓冲蛋白是固定不动的，只有大约四分之一的缓冲蛋白以低于Ca²⁺扩散速度进行缓慢地移动。在“开反应”中，缓冲蛋白主要与Ca²⁺结合，而在“关反应”中，Ca²⁺与缓冲蛋白分离。这个过程，可以改变Ca²⁺信号的幅值和恢复时间，以满足Ca²⁺信号的多样性。进入细胞内的Ca²⁺，只有一小部分会以自由离子的形式在细胞内扩散，大多数Ca²⁺会以较快的速度结合到缓冲蛋白上，少部分会与效应器结合。

在“关反应”中，钙泵发挥了巨大作用，它将Ca²⁺输送到细胞外或者泵到内质网中以维持细胞内[Ca²⁺]达到静息水平 [4] ，保持胞外和内质网中钙的高浓度 [26] 。细胞内的钙泵主要由两种，一种是位于细胞膜上将Ca²⁺由细胞内转运到细胞外的钙泵(PMCA)，另一种是位于内质网上将Ca²⁺由细胞内运送到内质网中的钙泵(SERCA)。钙泵是一种蛋白酶，它可以催化膜内的ATP水解，释放出能量，将细胞内的钙离子输送出去以维持细胞内钙离子的低浓度。这个过程中钙泵消耗的能力比较少，但是输送Ca²⁺的效率很高。

Ca²⁺信号的动力学多样性是由于其空间及时间特性变化而得到的 [2] [27] 。由于IP₃R通道在内质网膜上的集团化非均匀化分布，每个集团内的IP₃R通道大约为几个到几十个 [28] [29] 。根据IP₃R通道的打开的数量不同，Ca²⁺信号呈现出不同的时空特性。Ca²⁺信号中最根本的信号称为Blip [30] ，其形成主要是在较低的IP₃浓度下，只有少数的IP₃R通道结合IP₃分子，导致只有单个IP₃R通道打开，引起细胞内局域[Ca²⁺]的上升 [31] ；Ca²⁺信号中基本的信号为puff，主要是在中等IP₃浓度下，位于同一个集团内的一些IP₃R通道协同打开，释放出Ca²⁺。Ca²⁺信号中全局信号称为钙波(wave)，主要是在高IP₃浓度下，大量的IP₃R通道结合IP₃分子，基本事件释放的Ca²⁺扩散至附近的集团中，通过钙致钙释放机制使得各个集团内的IP₃R通道大量打开，形成胞内钙波。

为了详细地了解钙信号在生物系统中的作用，生物实验为我们提供了很多的洞见。但是由于钙信号调控机制的复杂性、生物实验技术和方法的限制性，我们还不能完全了解钙信号在生物系统中各种调控机制的详细过程。因此，基于实验过程的理论建模给我们提供了一个新的视角来研究钙信号的调控机制，大大提高了研究效率和定量性，计算建模已成为一个对复杂系统描述及预测的强大工具 [32] 。

现在已经有很多模型讨论不同生物过程中的钙动力学过程 [33] - [37] 。本文主要针对钙离子通道和钙信号震荡模型进行综述。本文第二章主要介绍了相关的数学基础，第三章综述了钙离子通道，主要是IP₃R通道和MCU通道的建模研究，第四章综述了内质网调控的Ca²⁺模型，第五章综述了内质网–线粒体耦合调控钙信号系统的动力学模型，最后为总结和展望。

2. 相关数学基础

数学建模就是将生物系统之间的相互作用用数学语言表现出来。下面就介绍几个建立生物模型的基本数学工具。

2.1. 微分方程建模

生物化学反应主要是通过常微分方程(Ordinary Differential Equations, ODEs)来表示，因为ODEs可以比较准确地描述系统动力学的变化，即各个参与反应物质的浓度随时间的变化 [38] 。ODEs构建的基本方法主要是质量作用定律。质量作用定律是19世纪G.M.古德贝格和P.瓦格提出的，主要内容是：化学反应速率与反应物的有效质量成正比 [39] 。其中有效质量是指反应物的有效浓度。

假设有两种反应物A，B相互作用，生成产物AB。其反应方程为：

， (2.1)

其中[A]，[B]分别代表反应物A，B的浓度，[AB]代表生成物AB的浓度；k₁和k₋₁分别为结合速率常数和解离速率常数。结合速率r₁和解离速率r₂分别表示为：

， (2.2)

。 (2.3)

对于反应物A，B的浓度变化主要是在正反应过程中减少，在逆反应过程中增多，而生成物AB的浓度变化则与A，B相反；用公式可以表示为：

， (2.4)

。 (2.5)

当反应达到平衡状态，即结合速率与解离速率相等，反应物与生成物地浓度都保持不变，可以得到：

， (2.6)

其中，K_d称为平衡常数(equilibrium constant)，又称为解离常数(dissociation constant)。K_d越大，说明反应物A，B的结合能力越弱，生成物AB越少，AB的解离速度更快。

2.2. 马尔科夫过程

马尔科夫过程是以俄国数学家安德烈·马尔可夫命名的随机过程，马尔科夫过程被认为是“无记忆”的，因为系统可以仅仅根据其目前的状态进程对未来作出预测，并不考虑之前的状态。

在生物系统建模中，微分方程可以用来表示生化反应过程。但是，对于随机性比较强的反应过程，微分方程就不能具体描述。因此，马尔科夫随机过程就成为解释生物过程中随机过程的一种重要的手段。对于钙信号模型，由于钙门控通道的随机性特点，其通道亚基的激活、失活和恢复过程就可以用马尔科夫过程来描述。

2.3. 参数优化

由于生物建模过程中有许多未定的参数，为了尽量增大模型的输出结果与实验结果的达到更好的拟合效果，尽量还原各种生物学反应。对于模型来说，参数优化是一个重要的课题。本节主要对模拟退火算法、遗传算法和免疫算法作简单的介绍。

2.3.1. 模拟退火算法

模拟退火算法成为了一种求解大规模组合优化的随机性方法。它是由美国物理学家N. Metropolis等人在1953年提出的 [40] ，1983年，S. Kirkpatrick等人将模拟退火思想引入到组合优化的领域中 [41] 。

模拟退火算法是基于蒙特卡洛迭代法，其思想来源于固体物质的退火过程。因为在自然规律中，物质总是趋向于最低能态，因为最低能态是最稳定的状态。

模拟退火算法的基本思想是先把固体加热至足够的温度，使固体中的所有粒子处于无序的状态，然后慢慢地将温度降低，使粒子活动逐渐有序，在常温时粒子达到最低能态，即寻找到最优解。为了防止陷入局部最优解，常常使用Metropolis准则。

Metropolis准则是根据概率的大小判断是否接受新的状态。在高温下，可以接受与当前状态能量差较大的新状态；而在低温下，只接受与当前状态能量差较小的新状态。

模拟退火算法具有可并行性、扩展性和通用性，它可以高效的解决几乎所有的组合优化过程。

2.3.2. 遗传算法

随着达尔文的进化论和孟德尔的遗传学说被广泛应用之后，在计算机领域，人们将这两种学说应用于算法中，得到了遗传算法。遗传算法是由美国Holland教授首先提出的 [42] ，是一种通过模拟自然进化过程随机搜索全局最优解的方法 [43] 。

遗传算法将问题的求解过程模拟为群体的生存过程，通过一代代的优胜劣汰的筛选进化，找出问题的最优解。它的基本思想是，首先随机产生初代种群，然后根据进化理论，逐代进行优胜劣汰的选择，产生越来越接近的最优解。在每一代中，根据适应度的大小来选择个体，并运用遗传算子对个体进行交叉和变异，改变个体的基因，从而产生新的解的解集。

由于遗传算法具有使用方便、可以并行计算、自适应强等优点，因而在函数优化、组合优化、人工智能等领域被广泛应用。但是遗传算法也存在一定的不足之处，比如遗传算子可能引起系统退化，遗传算子相对固定，缺少灵活性等。

2.3.3. 免疫算法

因为遗传算法存在着一定的不足，因此人们在保留遗传算法优点的基础上，为了避免其缺陷，引入人体的免疫系统机制，形成了免疫算法。免疫算法是人们模仿生物免疫系统基本机制，模拟免疫系统对外来病菌的识别能力而设计出来的多峰值搜索方法 [43] ，最早的免疫系统的模型是在1973年由Jerne提出的 [44] 。

免疫算法的基本思想是：对于具体要解决的问题(抗原)，首先是必须对问题进行具体分析，提取问题最基本的特征信息(疫苗)，对特性信息进行分析处理，得到待解决问题的一种方案(抗体)，然后对抗原和抗体进行亲和力分析，如果没有得到满意的结果，免疫算子进行具体的操作，形成新的抗体。免疫算子主要通过接种疫苗和免疫选择两个步骤完成，接种疫苗是为了确保提高亲和力，而免疫选择是为了防止群体退化 [45] 。

3. 通道的建模研究

由于细胞内含有多种钙离子通道，因此，对于钙离子通道的建模也出现了很多，本文主要以内质网上的IP₃R通道和线粒体上的MCU通道为主。下面就这两种通道的建模研究进行介绍，这里主要侧重介绍IP₃R通道的建模。

3.1. DeYong-Keizer IP₃R通道模型模型

1992年，DeYong和Keizer在PNAS上发表了关于IP₃R离子通道的模型 [46] 。根据内质网上IP₃R通道的激活和失活的定量测量，他们建立了一个动力学模型来反映IP₃R通道的特性。

模型主要的动力学过程包含两个假设。第一个是假设每个IP₃R通道主要是由3个完全相同的亚基构成。第二个是假设每个亚基由3个结合位点，一个IP₃结合位点，两个Ca²⁺结合位点，分别为Ca²⁺激活位点和抑制位点。每个结合位点有两种状态，激活(用1表示)或失活(用0表示)，因此每个亚基就有8个不同的状态(图1)，且亚基的8种状态可以相互转换，转换速率分别用不同的，表示。当IP₃结合位点和激活Ca²⁺结合位点被占据，抑制Ca²⁺结合位点未被占据，即亚基处于110态时，我们称亚基处于激活态；当一个通道中3个亚基全部处于激活态(110态)时，IP₃R通道打开。根据质量作用定律，描述亚基动力学的方程可以很容易的给出。

在此，我们主要列举了000态和110态的微分方程，其他6态就不详细描述：

，(3.1)

， (3.2)

其中，，这样，通过7个微分方程，可以很好地描述IP₃R通道的开放与关闭的动力学过程。

通道的开放概率可以表示为：

。(3.3)

DeYong-Keizer模型很好的模拟出了实验得到的IP₃R通道的开放概率与[IP₃]和[Ca²⁺]的结果。

Figure 1. The transformation of 8 states in subunits of IP₃R channels. a_1-5, b_1-5 donate the on-rates and off-rates, respectively. Ca donates [Ca²⁺] and I donates [IP₃]

图1. IP₃R通道中亚基8个状态之间的相互转换示意图。a_1-5，b_1-5分别表示各个转换速率，Ca表示Ca²⁺浓度，I表示IP₃浓度

3.2. Li-Rinzel IP₃R通道模型模型

1994年，Li和Rinzel在研究DeYong-Keizer模型的基础上，根据快慢时间尺度分离法将模型进行简化，形成了经典的Li-Rinze模型 [47] 。Li-Rinze模型将DYK模型中描述IP₃R通道的的7个动力学变量简化至只有1个变量的系统。

IP₃R通道的开放概率替代为：

， (3.4)

其中，，和项分别表示为：

， (3.5)

， (3.6)

， (3.7)

其中：

。 (3.8)

，和分别代表IP₃的结合，激活Ca²⁺的结合和抑制Ca²⁺的结合。由于IP₃和激活Ca²⁺的结合是一个快速的过程，而抑制Ca²⁺的结合是一个相对较缓慢的过程，因此用下列的微分方程表示：

， (3.9)

其中：

。 (3.10)

Li-Rinze模型在简化DeYong-Keizer模型的基础上，保留了DeYong-Keizer模型中大部分的动力学特征，明确了三个不同门控通道在时间尺度上的分层结构，详细地说明了IP₃R通道的动力学过程，从而大大降低了模型的计算量，但是模型的结果仍然与DeYong-Keizer 模型符合的很好。

Li-Rinze模型与DeYong-Keizer 模型都假设IP₃R通道含有3个相同的亚基，而Shuai等 [48] 在他们的模型中假定IP₃Rs通道是一个由4个相同且独立的IP₃R亚基组成的四聚体，其中的每个单体都由IP₃和Ca²⁺控制其激活或失活。他们假设当一个通道中至少有3个亚基处于激活态时IP₃Rs通道便处于开放状态。因此，IP₃R通道的开放概率可以表示为：

。(3.11)

3.3. 其他基于DeYong-Keizer模型简化的IP₃R模型

Shuai [49] 随后改进了DeYong-Keizer模型，建立了基于亚基或者IP₃R通道几种状态相互转换的简化的IP₃R模型。Shuai [49] 将DeYong-Keizer模型中IP₃R通道的8态模型，简化为6态模型(如图2所示)该模型得到的结果与DeYong-Keizer模型得到了很好的吻合。

Figure 2. The 6-state subunit structure of the channel model [49]

图2. 通道模型的6态亚基结构 [49]

其次，Shuai [49] 构建了关于IP₃R通道的7态模型，如图3所示。该模型是根据2007年，Shuai发表的9态模型简化而来 [50] 。模型进一步包括了一个构象变化，即在110态(IP₃位点和激活Ca²⁺位点被结合)的亚基是“非激活的”，必须通过构象变化到达A态，亚基才是激活态，可以促进通道的打开。图3是在9态模型上的简化，主要是去掉了101和001态。因为模型结果发现这两个态的变化对通道动力学的影响不大。

3.4. 整体单元IP₃R模型

上述DeYong-Keizer 模型和其简化和改进模型等，都假设IP₃R通道含有3个或4个相同且独立的IP₃R亚基。另一种建模方式是把IP₃R通道当作一个整体单元来建模，而不是针对其亚基建模。下面我们讨论该类整体单元模型。

1999年，Swillens等人 [51] 建立了一个18态的IP₃R通道模型(如图4所示)。该模型主要假定通道含有1个IP₃结合位点，两个激活Ca²⁺结合位点，两个抑制Ca²⁺结合位点。通道处于R₂₂₀状态时，即IP₃结合位点和两个激活Ca²⁺结合位点均被占据，则通道打开。

图5是Shuai [49] 建立的20态的IP₃R通道模型。在模型中，将IP₃R通道作为一个单元考虑，该单元有四个连续的IP₃结合位点，两个连续的激活Ca²⁺结合位点和一个抑制Ca²⁺结合位点。模型假设，只有当两个激活Ca²⁺结合位点被结合，抑制Ca²⁺结合位点未被结合的时候，即通道处于420状态时，IP₃R通道打开。该模型也可以很好的解释IP₃R通道的动力学特性。

3.5. Sneyd IP₃R通道模型

2001年，Sneyd等人 [52] 建立了一个关于IP₃R通道的10态模型(如图6(a)所示)。该模型虽然包含了一系列的状态，但是它的背景结构很简单(图6(b)所示)。R表示受体，它可以结合Ca²⁺到非激活态I₁，也可以结合IP₃分子到打开态O；打开态O可以到关闭态(S态)，也可以结合Ca²⁺并激活到激活态(A态)。A态可以结合Ca²⁺到抑制态I₂。

通过上述模型的研究，得到了Ca²⁺可以降低IP₃结合速率，同时可以提高受体与IP₃之间稳定状态的敏感性；Ca²⁺可以提高钙离子激活受体的速度，且速度高于抑制受体的速度等结论。

Figure 3. The 7-state subunit structure of the channel model [49]

图3. 通道模型的7态亚基结构 [49]

Figure 4. Swillens model of Ca²⁺ channel. The symbol Rijk refers to the state of the channel, to which i (0 or 2), IP3 molecules, j (0,1 or 2), Ca²⁺ ions at the activating sites, and k (0, 1 or 2), Ca²⁺ ions at the desensitizing sites are bound [51]

图4. Swillens钙通道模型。R_ijk表示通道的状态，i表示IP₃结合位点(0或2)，j表示激活Ca²⁺结合位点(0, 1, 2)，k表示抑制Ca²⁺结合位点(0, 1, 2) [51]

Figure 5. The state structure of the sequential binding IP₃R model [49]

图5. 连续结合的IP₃R通道模型的状态结构 [49]

(a) (b)

Figure 6. (a) The full IP₃R model. R, receptor; O, open; A, activated; S, shut; I, inactivated. c is [Ca²⁺]; p is [IP₃]. (b) Simplified diagram of the IP₃R model [52]

图6. (a) IP₃R通道模型。其中R代表受体；O是打开；A是激活；S是关闭；I是抑制。c是钙浓度，p是IP₃浓度。(b) IP3R模型的简化图 [52]

3.6. Siekmann IP₃R通道模型模型

2012年，Siekmann等建立了一种新的关于IP₃R通道的模型——“Park-Drive”模型 [53] 。Siekmann IP₃R模型是一个含有6个状态的马尔科夫模型，主要有两个模式，park模式和drive模式。每个模式都包含着开放和关闭状态(如图7所示)。整体模型包含两个部分，一部分是高活性部分(drive mode)，包括关闭状

Figure 7. The structure of the Siekmann IP₃R model. C represents closed state and O is open state; q is transition rates connecting two adjacent states [53]

图7. Siekmann IP₃R模型结构。C表示IP₃R通道关闭，O表示通道打开，q表示两种状态的转化速率 [53]

态C₁，C₂，C₃和打开状态O₆。在drive模式，通道几乎是打开的。另一部分是低活性部分(park mode)，主要有一个关闭状态C₄和一个打开状态O₅。在park模式，IP₃R通道几乎是关闭的。这个模型很好的模拟了单通道的随机动力学过程。

3.7. MCU通道模型

这里我们简单介绍Qi等 [48] 建立的关于线粒体上MCU通道的模型。MCU摄钙的驱动力为线粒体膜电位(线粒体基质与膜间隙之间巨大的电压差)。除此之外，实验上证明MCU摄钙还可由MICU1所控制。综合上述两点，可得到通过MCU的Ca²⁺流为

。 (3.12)

v_MCU表示通过MCU的最大流量，ΔΦ代表电压驱动力，R_MICU代表MICU1调控项，P_oMCU为MCU的打开概率。

模型假设在线粒体摄钙时膜电位为一恒定值(Ψ = 170 mV)。MCU驱动力使用Gunter等 [54] 给出的公式：

， (3.13)

其中F为法拉第常数，R为气体常数，T为绝对温度，b和Ψ₀为拟合参数 [54] 。

MICU1作为MCU的一个调控亚基，通过感受线粒体膜间隙之中的[Ca²⁺]而为MCU摄钙设定一个阈值 [55] 。实验结果表明只有当线粒体外的[Ca²⁺]超过3 μM时MCU才会显著摄钙。所以，MICU1对MCU活性的调控有一个陡峭的激活常数(k_MICU)，其表观希尔系数为4，得到

， (3.14)

这里的[Ca²⁺]_Mic为MCU感受到的钙微域中的[Ca²⁺]。

模型假设每个MCU单体有两个Ca²⁺结合位点，一个为激活位点，另一个为失活位点。所以每个MCU单体可以s₀₀，s₁₀，s₀₁和s₁₁这4种状态存在，这里第一个下标代表激活Ca²⁺位点，第二个下标代表失活Ca²⁺位点，1表示处于结合态，0表示处于空闲态(图8)。与IP₃Rs通道类似，只要有3个以上MCU单体处于激活态，即s₁₀态MCU通道便处于开放状态，于是

MCU的开放概率为：

， (3.15)

这里的x_ij表示一个单体处于激活态的概率。控制MCU单体状态转化的动力学方程组为：

， (3.16)

且x₀₀ + x₀₁ + x₁₀ + x₁₁ = 1。该MCU模型可以拟合通道电压控和受体控两种调控方式。

4. 内质网调控Ca²⁺模型的研究

在本综述中，我们将忽略钙信号的空间扩散行为，而关注于钙信号的震荡动力学，所以我们将只讨论钙信号震荡的点源模型，不考虑空间扩散动力学，所以其数学描述仅仅是常微分方程。本章主要讨论内质网调控Ca²⁺模型。

4.1. 全局钙震荡模型

在DeYong-Keizer模型和Li-Rinze模型中，细胞质中的钙浓度是随着时间震荡的，而IP₃浓度是不变的。细胞质内的钙浓度都是由IP₃R通道、漏通道和钙泵三项构成，其微分方程为：

Figure 8. The four-state dynamic model of MCU monomer [48]

图8. MCU的4态动力学模型 [48]

， (4.1)

其中，[Ca²⁺]表示细胞质内自由的Ca²⁺浓度，J₁表示由于IP₃R通道打开引起的Ca²⁺的外流，J₂表示由漏通道流出的Ca²⁺。J₁和J₂项是表示Ca²⁺的外流，引起细胞溶质中的[Ca²⁺]升高。J₃表示通过钙泵将Ca²⁺逆浓度梯度从细胞质中输送至内质网中，是引起细胞溶质内[Ca²⁺]减小的过程。

因此，J₁项由方程可以表示为：

， (4.2)

其中，v₁表示IP₃R通道的最大流量。c₁表示内质网与细胞质的体积比。表示内质网中的[Ca²⁺]。P_open则是前章所计算的通道开通几率。

J₂项可以用方程表示为：

， (4.3)

其中，v₂表示漏通道的最大流量。

内流的Ca²⁺主要是通过钙泵(SERCA)将Ca²⁺逆浓度梯度输送到内质网中，其表达式主要是通过希尔方程表示：

。 (4.4)

根据上述的数学模型，DeYong-Keizer模型和Li-Rinze模型都可以得到了钙的震荡曲线。图9是DeYong-Keizer模型得到的分岔图。

4.2. Rudiger-Shuai局域钙模型

2010年，Rudiger和Shuai在PRL上发表了一篇关于运用决定性和随机性模拟方法模拟单个集团内几个IP₃Rs通道释放局域钙信号的网络模型 [56] 。在模型中，9个IP₃R通道均匀地分布在立方体的一个表面上(如图10(a)所示)。由于IP₃R通道集团相对较大，并且Ca²⁺扩散速度相对较快，同一集团内的IP₃R通道往往处于不同浓度的Ca²⁺环境中，因此他们对于关闭和开放的IP₃R通道引入了不同的[Ca²⁺]。对于一个开放的IP₃R通道，其通道口的[Ca²⁺]可以达到100 μM甚至更高；而对于关闭的IP₃R通道取决于其他打开的IP₃R通道。

图10(b)可以看出，关闭通道口的平均[Ca²⁺]与打开通道数成线性关系。因此，可以用如下公式表示：

， (4.5)

其中，表示某一关闭通道得到其他n个打开通道扩散的[Ca²⁺]的平均值。

根据上述公式，Rudiger-Shuai模型很好的模拟出了实验中所观测到的局域钙Puff信号幅度分布的结果，为之后的理论建模提供了新的思路。

4.3. Qi局域钙模型

2014年，Qi建立了一个双Ca²⁺浓度的模型研究IP₃R通道集团内的产生的局域钙blip (单个通道打开)和puff (多通道打开)信号动力学行为 [57] 。该模型由一个简单的钙动力学方程和一个马尔科夫过程决定。对于内置网上的Ca²⁺通道，在该模型中只考虑IP₃R通道对Ca²⁺动力学的影响。简化模型的同时，还考虑了Ca²⁺结合蛋白对Ca²⁺动力学的影响。对于IP₃R通道的开关动力学，主要是由修饰过的DeYong-Keizer模型来讨论，其中亚基的状态变换是由马尔科夫动力学决定。

Figure 9. The bifurcation diagram. The diagram indicates stable and unstable steady states, and gives the maximum and minimum cytosolic [Ca²⁺] along stable periodic orbits [46]

图9. 分岔图。图中显示了稳定和不稳定状态，给出了在稳定轨道周期内细胞溶质钙浓度的最大值和最小值 [46]

(a)(b)

Figure 10. Calcium release. (a) The box of dimension 8*8*5 m³ represents the cytosolic space. IP₃R channels are located in a grid on the surface. (b) Calcium concentrations obtained by averaging over closed channels for a given total number of open ones [56]

图10. Ca²⁺的释放过程。(a) 8*8*5 m³的立方体表示细胞质溶质，IP₃R通道分布于立方体的某一表面上。(b) 所有关闭通道口的平均[Ca²⁺]与打开通道数的关系 [56]

对于打开和关闭的IP₃R通道周围的Ca²⁺浓度分别采用不同的[Ca²⁺]。对于开放的IP₃R通道，它将产生一个很高的自抑制浓度，在模型中，设置[Ca²⁺] = 120 μM。而对于关闭的IP₃R通道，根据Rudiger-Shuai模型 [56] 的推导，在平衡状态下关闭通道平衡时周围的[Ca²⁺]可以表示为：

， (4.6)

其中，[Ca²⁺]_Rest为静息状态下细胞溶质中的[Ca²⁺]，∆C为打开通道的Ca²⁺扩散到邻近通道，引起[Ca²⁺]上升的平均值，N_Open为打开的通道数。

当一个通道由打开状态转向关闭状态，达到平衡状态的[Ca²⁺]_Equil之前，[Ca²⁺]会受到Ca²⁺的扩散和Ca²⁺缓冲蛋白的影响。因此，对于关闭状态的[Ca²⁺]，可以用微分方程表示为：

，(4.7)

其中，γ为Ca²⁺的衰减率。

根据上述模型，可以得到集团数目和[IP₃]的增多会减少blip的比例和持续时间。相邻局域钙信号之间的时间间隔及第一次释放事件的响应延迟时间与集团大N的倒数成线性关系等结论。图11是上述模型得到的[IP₃]对blip比例和持续时间的影响。

4.4. Cao局域钙模型

2013年，Cao发表了一篇在单通道数据的基础上得到的puff信号的随机模型 [58] 。该模型的IP₃R通道模型是引用了Siekmann的IP₃R模型 [53] 。Cao的模型除了考虑IP₃R通道对Ca²⁺动力学影响，还考虑了漏通道、快速缓冲蛋白和慢速缓冲蛋白dye(fluo-4)对钙信号的影响。在模型中，主要的动力学微分方程可以表示为：

， (4.8)

， (4.9)

其中，表示可以是细胞质内钙浓度升高的钙流。表示使Ca²⁺浓度减少的钙流，主要是通过Ca²⁺的扩散和钙泵的影响。表示内质网上漏通道的影响。和分别表示dye的浓度和与Ca²⁺结合的dye的浓度。表示打开的通道数。

Cao根据以上的模型，发现在小IP₃R集团内，puff的平均幅值与集团大小成线性关系，而对于较大IP₃R集团来说，呈非线性关系。Puff的衰减时间对集团大小敏感性比上升时间对其敏感性更强。

4.5. 全局钙和IP₃震荡模型

前面介绍的模型中，IP₃的浓度是不随时间改变的一个常数值。而对于有些类型的细胞，其IP₃的浓度是由钙离子浓度所影响的，这样钙离子浓度的震荡将导致IP₃震荡，这里我们也讨论其相关的建模。

Sneyd等人 [59] 在PANS上发表一篇关于IP₃浓度影响钙震荡的模型。该模型主要包含了一个钙震荡方程和一个IP₃震荡方程，表示如下：

， (4.10)

， (4.11)

，(4.12)

。 (4.13)

(a) (b)

Figure 11. (a) The blip fraction and (b) the blip duration against [IP₃] at different cluster size [57]

图11. (a) blip事件的比例和(b) blip持续时间在不同集团大小情况下随IP₃分子浓度的变化 [57]

其中，[Ca²⁺]表示全局钙浓度，[Ca²⁺]_ER表示内质网中的钙浓度，[IP₃]表示IP₃的浓度，n表示IP₃受体的比例。通过该模型，他们很好的模拟出了在胰腺腺泡细胞中得到的IP₃浓度影响钙震荡的实验结果(如图12所示)。

Shuai等人 [60] 在文章中用了另一种方程表示IP₃分子的生成。可以用方程表示为：

， (4.14)

其中，IP₃的降解J_D表示为：

， (4.15)

IP₃的生成可以表示为：

。 (4.16)

5. 内质网-线粒体调控Ca²⁺模型

前一章主要描述了内质网释放的Ca²⁺模型，分析了不同模型对内质网钙信号的研究。线粒体作为钙离子存储的另一个主要钙库，其存储和释放钙离子，也会引起钙信号的变化。研究线粒体参与调控的Ca²⁺模型相对较少，下面我们就内质网和线粒体耦合调控钙信号系统的动力学模型进行介绍。

5.1. Fall-Keizer模型

2001年，Fall和Keizer整合他们已经发表的多个钙信号模块，研究了线粒体对细胞质钙信号的影响，以及细胞质钙信号对线粒体氧化呼吸功能的影响 [61] 。他们的模型主要包含以下四个模块：1) 内质网钙信号成分，2) 线粒体钙信号成分，3) 线粒体代谢成分，和4) 细胞质钙信号及ATP代谢成分。在这四个模块之中，模块1与模块4有Ca²⁺流的交换，模块2与模块4也有Ca²⁺流的交换。在该模型中，线粒体上的MCU直接感受的是细胞质中的平均[Ca²⁺]。

Fall和Keizer对MCU采用了Monod-Wyman-Changeux模型的方法，即假设MCU为变构调节酶类

Figure 12. Model and experimental responses to IP₃ pulses of increasing magnitude [59]

图12. IP₃脉冲幅值增加引起模型和实验的结果 [59]

型的通道，其数学表达式为：

， (5.1)

其中，

。 (5.2)

这里，CAM和CAC分别为线粒体和细胞质中[Ca²⁺]，MWC为Monod-Wyman-Changeux模型的表达形式，公式(5.1)中包含了线粒体膜电位的作用，其中F为法拉第常数，R为气体常数，T为绝对温度，为线粒体膜电压。

该模型引入线粒体机制，使IP₃模型具有了双稳态特性，并且代谢底物的增加降低了震荡的频率；与线粒体相关的细胞溶质内钙信号，可以诱导代谢的变化。图13是模型得到的线粒体影响细胞溶质的钙浓度变化结果。

5.2. Marhl模型

Marhl等 [62] 于2000年构建了包含内质网、线粒体以及Ca²⁺结合蛋白的简单模型，主要通过数学手段模拟了一些复杂的钙震荡行为。

Figure 13. Effect of adding mitochondria to IP₃-mediated Ca²⁺ release model. The dotted line shows the cytoplasmic Ca²⁺ concentration ([CAC]) for the De Young- Keizer model tuned to bistable behavior. At the same IP₃ system parameter values, the solid line shows the effect of addition of mitochondria [61]

图13. IP₃调控Ca²⁺释放模型加入线粒体的影响。虚线是De Young-Keizer模型得到的细胞溶质钙浓度达到稳态值。在同样的参数下，实线是加入线粒体的影响结果 [61]

该模型主要集中于细胞溶质和内质网、线粒体及Ca²⁺结合蛋白之间的钙交换。模型中主要包括5个变量，细胞溶质内的自由[Ca²⁺]([Ca²⁺])，内质网中的自由[Ca²⁺]([Ca²⁺]_ER)，线粒体中的自由[Ca²⁺]([Ca²⁺]_m)，细胞溶质中自由的钙结合蛋白浓度([Pr])及结合到钙结合蛋白中的[Ca²⁺]([CaPr])。则细胞溶质内的[Ca²⁺]可以表示为：

。 (5.3)

他们利用简单的希尔函数来表示Ca²⁺从细胞质进入线粒体：

， (5.4)

其中K₂ = 0.8 μM，为线粒体摄钙激活常数。

根据上述模型，Marhl等发现80%从ER中释放的Ca²⁺会迅速被线粒体吸收，之后会缓慢的释放。该模型也可以很好的解释钙的震荡及其混沌过程。

5.3. Szopa模型

2013年Szopa等 [63] 在Marhl模型的基础上，考虑了线粒体与内质网紧密关联膜(mitochondria-as- sociated ER membrane, MAM)在MCU摄钙过程中所起的作用(见图14)。他们将进入线粒体的Ca²⁺流分成两部分：一部分为位于MAM区域的MCU所摄取的Ca²⁺，这里的MCU感受的是一个很高的微域[Ca²⁺]，他们假设MAM区域内的[Ca²⁺]与内质网内[Ca²⁺]相当；另一部分为位于MAM区域以外的MCU所摄取

Figure 14. The diagram of Szopa model. The arrows point the flow direction of Ca²⁺. Cyt, cytosol; Mit, mitochondria; Er, endoplasmic reticulum; Pr, calcium binding protein [62]

图14. Szopa模型示意图。箭头为Ca²⁺流动方向。Cyt，细胞质；Mit，线粒体；Er，内质网；Pr，Ca²⁺结合蛋白 [62]

的Ca²⁺，这里的MCU感受的是平均的细胞质[Ca²⁺]。此外，他们假设MCU有两种工作状态：一种为标准的MCU状态，用希尔系数为2的希尔函数所描述；一种为RaM (RaM为线粒体上可能存在的一种摄钙通道，其是否存在现在尚有争议)状态，用希尔系数为8的希尔函数所描述，即

(5.5)

和

(5.6)

其中，J_MAM为通过MAM区域内MCU的Ca²⁺流，J_in为通过其它非MAM区域内MCU的Ca²⁺流。

通过Szopa模型，他们发现数值稳定震荡存在相当大的一组参数值。某些参数可以使震荡消失，并且模型的轨迹倾向于稳定状态，且线粒体钙水平非常高 [62] 。

5.4. Qi内质网-线粒体模型

Qi等人 [48] 在2015年建立了关于内质网-线粒体的模型。该Ca²⁺模型是一个闭合的细胞模型(如图15所示)，并没有考虑细胞质与细胞外环境中Ca²⁺的交换。模型包含3个区室(compartment)：细胞质(Cyt)，内质网(ER)和线粒体(Mt)。

3个区室中游离Ca²⁺的动力学将由区室间的Ca²⁺流以及各区室中Ca²⁺与Ca²⁺结合蛋白(binding protein, BP)之间的Ca²⁺流所决定，因此该模型的基本数学表达式为3个耦合的微分方程：

(a) (b)

Figure 15. (a) The schematic diagram of the components and fluxes in the calcium signaling model. (b) IP3Rs and MCUs in close proximity form a microdomain. A high microdomain [Ca²⁺] ([Ca²⁺]_Mic) can be generated upon opening of IP3Rs which in turn leads to Ca²⁺ uptake into mitochondria through MCUs [48]

图15. (a)钙信号模型的成分和钙流示意图，(b)近距离的IP3Rs和MCUs构成一个微域。由于IP3Rs通道的打开，会生成一个高浓度微域钙(Ca²⁺] ([Ca²⁺]_Mic)，进一步导致MCUs摄取钙进入线粒体内。

， (5.7)

， (5.8)

。 (5.9)

在这里，[Ca²⁺]代表游离Ca²⁺的浓度，V表示3个区室的体积，jⁱⁿ代表从外向内的Ca²⁺流，而j^out代表从内向外的Ca²⁺流，j_Cai这3项表示Ca²⁺与BP解离的流。

下面通过模块化的方式描述内质网、线粒体和钙微域的动力学，其中内质网的Ca²⁺流主要是由IP₃R通道，漏通道和钙泵三项构成，与4.1节全局钙震荡模型相同，此处不再赘述。而线粒体的Ca²⁺流主要经Na⁺-Ca²⁺交换体排出到细胞质中，线粒体中Ca²⁺的外流就可以表示为

， (5.10)

其中，v_NCX是NCX的最大活性，k_Na和k_NCX都是NCX的激活常数。

本模型的重点是建立一个代表相对位置极近的内质网与线粒体之间的钙微域模型(图15(b))。钙微域内相互作用的钙信号系统成分为IP₃Rs通道和MCU。由IP₃Rs通道集团释放Ca²⁺产生的电流就可表示为：

(5.11)

其中，P_Open是IP₃Rs通道的开放概率；n_IP3R是IP₃Rs通道集团内IP₃Rs通道的数目。开放的IP₃Rs通道释放的Ca²⁺，由于扩散作用会形成一个从IP₃Rs通道集团的中心到附近细胞质的陡峭的浓度梯度，靠近释放点的[Ca²⁺]要远高于[Ca²⁺]_Cyt，称为[Ca²⁺]_Mic，其表达式为：

， (5.12)

根据上述的模型，他们重复出了实验上得到的钙震荡，得到了线粒体可以明显减少细胞质钙震荡的振幅，IP₃R和MCU的距离对钙震荡的调节作用等一系列结果。图16是上述模型得到的IP₃R和MCU的距离对钙震荡的影响。

6. 总结和展望

本文介绍了生物建模过程中用到的相关建模理论，综述了目前已构建的IP₃R和MCU通道模型、内质网钙离子模型及内质网–线粒体相互作用钙离子模型。这些理论研究加深了我们对钙离子信号动力学行为的了解，为进一步研究钙信号在细胞生理过程中扮演的重要角色提供了理论依据。

由于钙离子是细胞内普遍存在的第二信使，控制着细胞大量的生理过程。近年来，对钙离子信号的研究也越来越深入，且在应用上也有了长足的发展。

钙离子可以调控基因的表达机制 [64] - [66] 。实验表明细胞内的钙离子震荡可以增加基因表达的有效性和特异性 [67] ，且胞内的Ca²⁺和核内的Ca²⁺都通过不同的机制影响着基因的转录。将钙离子的震荡动力学过程与基因表达机制相结合，从理论上研究钙离子对基因表达机制调控的影响，尤其是钙震荡对基因表达的随机性作用都需要进一步的研究。

钙离子信号对细胞凋亡的调控有一定的作用。钙信号参与细胞的凋亡过程已经有很多文章指出 [68] - [71] ，但是对于钙信号作用于凋亡过程的具体机理仍不清楚。已知钙离子有三条主要的信号通路 [72] ：

(a) (b)(c) (d)

Figure 16. Ca²⁺ dynamics modulated by the IP₃R-MCU distance [48]

图16. IP₃R-MCU距离对Ca²⁺动力学的调节 [48]

内质网通路 [73] - [75] ，线粒体通路 [76] - [78] 及死亡受体通路 [79] - [81] 。这三条通路都与钙离子相关，且最终都可以激活Caspases-3，使细胞凋亡。对于钙信号如何与凋亡的信号通道相互作用以及钙信号在不同凋亡途径中的作用，仍需要进一步的探索。将上述的Ca²⁺动力学模型结合到细胞凋亡网络中进行研究，是一种新的思路。这样我们不仅可以从理论上清楚研究钙信号和细胞凋亡的具体耦合过程，也可以从应用上将钙离子信号作为靶向目标，为临床癌症治疗提供一定的理论依据。

钙离子结合介导的钙调控异常不仅会导致肿瘤的产生 [82] [83] ，也会引起神经性疾病 [84] 。研究表明，钙离子信号可以调节学习记忆的过程，细胞内Ca²⁺信号的改变可以影响神经元的活性和功能 [85] 。阿尔兹海默病(AD)和帕金森病(PD)等慢性疾病都以神经细胞退行性病变和衰老为特征。现在已经有很多的研究针对钙信号对AD和PD的影响 [86] - [89] ，在临床上也已经使用了一些钙拮抗剂对AD等慢性疾病进行治疗。因此，研究神经系统钙信号通路的变化，以及其对学习通路机制的影响，可以为神经性疾病寻求新的发现，同时也为相应的药物治疗提供理论基础和新的思路。

在心肌细胞中，钙离子可以精密调控细胞的各种生理过程 [90] 。心脏的正常活动取决于钙循环的动态平衡，但是对于钙循环的确切机制仍存在着很多争议。尤其是目前人们对于微观钙循环的研究仍处于探索阶段。

除了上述的钙离子信号的相关研究，钙离子信号对糖尿病、高血压等疾病也都有一定的影响。钙蛋白酶系统的异常也会引起糖尿病 [91] ，骨骼肌萎缩 [92] 等疾病。因此，钙离子及其相关蛋白在人体中发挥着巨大的作用。随着对这些情况下的钙信号进行建模，并讨论相关的动力学行为，将会帮助我们更清楚地认识钙信号调控各种生命活动具体机制，并为各种疾病提供新的治疗方案。

致谢

本课题获得如下基金项目支持：国家自然科学基金面上项目(11504214和31370830)，和福建省高校领军人才资助计划。

^*通讯作者。