视动反应中PF-PC突触钙级联机制：物质状态预测

Calcium Cascade of PF-PC Synapses in Optokinetic Response: Prediction of Material State

摘要: 小脑的运动学习记忆所涉及的PF-PC突触权重变化被称为钙级联机制，这种机制利用化学受体磷酸化和去磷酸化的动态调节作用来模拟突触权重的变化。鉴于钙级联涉及多种化学物质组成的复杂动力学系统，在计算刺激信号是随机输入的视动反应时耗时较大，本文创新地采用与常微分方程性质相似的长短期记忆网络来预测短周期连续状态序列，从而在保证准确性的同时，有效减少了计算成本。不仅对小脑运动学习记忆形成有着理论意义，而且为小脑控制器提供工程上的应用价值。

Abstract: The change of PF-PC synaptic weight involved in cerebellar motor learning and memory is called calcium cascade. This mechanism uses the dynamic regulation of chemical receptor phosphorylation and dephosphorylation to simulate the change of synaptic weight. Since the calcium cascade involves a complex dynamic system composed of a variety of chemicals, it takes a lot of time to calculate the optokinetic response when the stimulus signal is random input. In this paper, the short-term and long-term memory network with similar properties to ordinary differential equations is innovatively used to predict the short-term continuous state sequence, so as to ensure the accuracy and effectively reduce the computational cost. It not only has theoretical significance for the formation of cerebellar motor learning and memory, but also provides engineering application value for cerebellar controller.

1. 引言

在小脑皮质中，平行纤维与浦肯野细胞间的突触连接强度变化，如长时程抑制和增强，对小脑学习至关重要。小脑长时程抑制(Long-Term Depression, LTD)发生在平行纤维和攀缘纤维同时刺激浦肯野细胞时，导致突触后钙离子浓度升高，突触效能减弱[1] [2]。这一过程的关键在于钙离子和肌醇1,4,5-三磷酸(IP₃)共同激活IP₃受体钙释放门控通道[3]。相反，长时程增强(Long-Term Potentiation, LTP)则增加兴奋性突触权重，通常发生在单独激活平行纤维时，此时IP₃受体通道未激活。LTP与LTD在小脑突触中需保持平衡，以支持运动学习的灵活性和逆转。

在海马体等结构中，平行纤维介导的LTP通过特定的生化机制实现。蛋白质磷酸酶对AMPA受体的去磷酸化作用增强突触传递效率，形成记忆痕迹。有研究指出诱发LTP的是钙调素依赖性蛋白激酶II (CaMKII)与N-甲基-D-天冬氨酸受体(N-methyl-D-aspartic acid receptor, NMDAR)亚单位2B (GluN2B)的结构性结合[4]。CaMKII有两种类型分别是αCaMKII和βCaMKII。αCaMKII通过突触传输介导的钙内流被激活，其随后的磷酸化在突触可塑性中起到核心作用[5]。而LTD则涉及AMPA受体的磷酸化和内吞作用，降低突触传递强度。这两种机制共同构成神经网络中学习与记忆的基础。

当平行纤维(Parallel Fiber, PF)激活先于攀缘纤维(Climbing fiber, CF)激活50~200 ms时，在平行纤维–浦肯野细胞(Parallel fibers-Purkinje cells, PF-PC)突触处诱导小脑长时程抑制所必需的大Ca²⁺信号，这与小脑学习理论一致。然而，大的Ca²⁺信号和/或LTD也可以通过单独的大规模PF刺激或通过笼状Ca²⁺+或IP₃的光解来诱导[6]。IP₃水平和钙离子动力学精确调控突触可塑性。此外，CaMKII作为钙信号系统的关键介质，参与介导小脑中的突触可塑性，促进学习。尽管突触后可塑性的信号通路已得到广泛研究，但多个生化通路如何共同影响双向可塑性仍需进一步探索。

尽管βCaMKII在小脑内表达丰富，其在学习过程中的精确功能和机制尚待深入挖掘。利用基因敲除技术，科学家成功培育了βCaMKII缺陷的小鼠，为探索该蛋白在PC突触可塑性中的作用提供了宝贵工具。研究表明，βCaMKII对PF-PC突触的可塑性变化具有关键作用。正常情况下，PF-PC突触可响应特定刺激产生LTD，优化运动学习。然而，在βCaMKII缺失的小鼠中，原本应诱发LTD的刺激却导致了LTP，突显了βCaMKII在维持正常LTD中的重要性[7] [8]。但这一转变背后的生化机制仍不清楚。

Woerden等人[9]推测，αCaMKII与βCaMKII间的生化差异可能是可塑性方向转换的关键。βCaMKII独特的F-actin结合能力可能使其与AMPA受体磷酸化过程相隔离。Pinto [10]等人进一步通过数学模型深入剖析了这一过程，揭示了F-actin与βCaMKII的相互作用如何精确调控PF-PC突触的可塑性方向。这一动态调控系统不仅影响长时程抑制或增强的发生，更对小脑神经网络的整体活动状态起着微调作用。通过不断的研究与探索，我们有望更加清晰地揭示βCaMKII在小脑学习中的核心作用及其精确调控机制。

2. 实现方法

2.1. PF-PC钙级联通路

钙级联通路结构如图1分成四个部分，一是平行纤维刺激诱导产生IP₃，IP₃打开内质网的IP₃R介导的门控通道；二是通过攀缘纤维的激活，引发浦肯野细胞膜的去极化反应，Ca²⁺流入浦肯野细胞的突触棘部；三是Ca₄CaM参与CaMKII活跃化的化学反应，使得α-氨基-3-羟基-5-甲基-4-异恶唑丙酸受体(α-amino-3-hydroxy-5-methyl-4-isoxazole-propionicacid receptor, AMPAR)磷酸化(如梯形区域所示)；四是Ca₄CaM参与蛋白磷酸酶2B (PP2B)的转化，促进AMPAR的去磷酸化(如椭圆区域所示)。决定磷酸化还是去磷酸化的关键在于作为钙库的内质网有无通过IP₃R门控通道释放大量的钙，把打开IP₃R门控通道比作化学反应，那么IP₃是催化剂，减少“化学反应”所需的“活化能”，即浦肯野细胞树突中的Ca²⁺。

Figure 1. Molecular pathways of synaptic plasticity in Purkinje cells

图1. 浦肯野细胞突触可塑性的分子通路

对于长时程抑制，PF刺激的去向分成mGluR通路和AMPAR通路。在代谢型谷氨酸受体(mGluR)通路中，PF末端释放谷氨酸(Glu)刺激mGluR，诱导质膜内肌醇磷脂(磷酸肌醇)水解产生IP₃，IP₃扩散到细胞质基质中；而在AMAPR通路中，平行纤维的刺激触发AMPAR激活，导致去极化并开启电压门控Ca通道(VGCC)通道，Ca²⁺流入胞质溶胶。PF刺激后紧接着CF刺激产生去极化，打开VGCC通道Ca²⁺内流进入浦肯野细胞树突棘。在IP₃的帮助下Ca²⁺的量已经达到打开IP₃受体钙门控通道的阈值，大量的Ca²⁺被释放出来，与CaM结合成Ca₄CaM后，促进CaMKII转为激活状态导致AMPAR磷酸化。CaMKII有四种状态，依次是抑制、绑定、磷酸化和自激活。尽管这四种状态均有可能与肌动蛋白相互作用，但只有绑定、磷酸化和自激活三种状态可以促进AMPAR的磷酸化反应，如虚线框所示[10]。图中褐色箭头代表在CaMKII的状态切换中，Ca₄CaM作为反应物，绿色箭头代表Ca₄CaM作为生成物。

对于长时程增强，PF刺激的去向相同，但由于缺乏后续的CF刺激，进入浦肯野细胞树突棘的Ca²⁺没有达到IP₃R通道的打开阈值，这些少量的钙参与Ca₄CaM介导的PP2B的激活，促进AMAPR去磷酸化。

2.2. 钙级联通路通信机制

参考Pinto [10]的动力学模型，CaMKII存在以下质量守恒关系，这里使用Pinto的简化表达以及符号，CaMKII用W表示，[ ]代表括号内物质的浓度，相关参数参照表1：

$\left[{\text{W}}_{\text{i}}]=[{\text{W}}_{\text{tot}}]-[{\text{W}}_{\text{b}}]-[{\text{W}}_{\text{p}}]-[{\text{W}}_{\text{a}}]-[{\text{W}}_{\text{i}}\text{Ac}]-[{\text{W}}_{\text{b}}\text{Ac}]-[{\text{W}}_{\text{p}}\text{Ac}]-[{\text{W}}_{\text{a}}\text{Ac}\right]$ (1)

其中W_tot代表CaMKII的总浓度。CaMKII_b和CaMKII_bAc分别磷酸化成CaMKII_p和CaMKII_pAc的浓度变化率计算如下：

${\text{V}}_{\text{a}}=\frac{{\text{K}}_{\text{a}}{\left({\text{c}}_{\text{b}}\left[{\text{W}}_{\text{b}}\right]\right)}^2+\left({\text{c}}_{\text{b}}\left[{\text{W}}_{\text{b}}\right]\right)\left({\text{c}}_{\text{p}}\left[{\text{W}}_{\text{p}}\right]\right)+\left({\text{c}}_{\text{b}}\left[{\text{W}}_{\text{b}}\right]\right)\left({\text{c}}_{\text{a}}\left[{\text{W}}_{\text{a}}\right]\right)}{{\left[{\text{W}}_{\text{tot}}\right]}^2}$ (2)

其中c_b、c_p和c_a是与每个活动状态的激酶活性成比例的加权因子。K_a是经验立方函数，用来模拟相邻的自磷酸化，K_a的计算如下：

${\text{K}}_{\text{a}}={\text{<msup> K <mo>′</mo> </msup>}}_{\text{a}}\left({\text{aT}}_{\text{a}}+{\text{bT}}_{\text{a}}^2+{\text{cT}}_{\text{a}}^3\right)$ (3)

其中CaMKII自磷酸化的现象速率为${\text{<msup> K <mo>′</mo> </msup>}}_{\text{a}}$ , a, b和c为调整参数，T_a为活跃CaMKII的总比例：

${\text{T}}_{\text{a}}=\frac{\left[{\text{W}}_{\text{b}}]+[{\text{W}}_{\text{p}}]+[{\text{W}}_{\text{a}}\right]}{\left[{\text{W}}_{\text{tot}}\right]}$ (4)

W_b、W_p和W_a的浓度变化率计算如下：

$\begin{array}{c}\frac{\text{d}\left[{\text{W}}_{\text{b}}\right]}{\text{dt}}=-{\text{V}}_{\text{a}}\left[{\text{W}}_{\text{tot}}\right]+{\text{k}}_{\text{ib}}\left[{\text{W}}_{\text{i}}\right]\left[{\text{Ca}}_4\text{CaM}\right]-{\text{k}}_{\text{bi}}\left[{\text{W}}_{\text{b}}\right]\cdots \\ +{\text{k}}_{\text{bacb}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]-{\text{k}}_{\text{bbac}}\left[{\text{W}}_{\text{b}}\right]\left[\text{Ac}\right]-{\text{k}}_{\text{fphos}}\left[{\text{W}}_{\text{b}}\right]\left[\text{AMPAR}\right]\cdots \\ +\left({\text{k}}_{\text{bphos}}+{\text{k}}_{\text{actphos}}\right)\left[{\text{W}}_{\text{b}}\text{AMPAR}\right]\end{array}$ (5)

$\begin{array}{c}\frac{\text{d}\left[{\text{W}}_{\text{p}}\right]}{\text{dt}}=-{\text{V}}_{\text{a}}\left[{\text{W}}_{\text{tot}}\right]+{\text{k}}_{\text{ap}}\left[{\text{W}}_{\text{p}}\right]\left[{\text{Ca}}_4\text{CaM}\right]-{\text{k}}_{\text{pa}}\left[{\text{W}}_{\text{p}}\right]\cdots \\ +{\text{k}}_{\text{pacp}}\left[{\text{W}}_{\text{p}}\text{Ac}\right]-{\text{k}}_{\text{ppac}}\left[{\text{W}}_{\text{p}}\right]\left[\text{Ac}\right]-{\text{k}}_{\text{fphos}}\left[{\text{W}}_{\text{p}}\right]\left[\text{AMPAR}\right]\cdots \\ +\left({\text{k}}_{\text{kbphos}}+{\text{k}}_{\text{kactphos}}\right)\left[{\text{W}}_{\text{p}}\text{AMPAR}\right]\end{array}$ (6)

$\begin{array}{c}\frac{\text{d}\left[{\text{W}}_{\text{a}}\right]}{\text{dt}}=-{\text{k}}_{\text{ap}}\left[{\text{W}}_{\text{p}}\right]\left[{\text{Ca}}_4\text{CaM}\right]-{\text{k}}_{\text{pa}}\left[{\text{W}}_{\text{p}}\right]+{\text{k}}_{\text{aaca}}\left[{\text{W}}_{\text{a}}\text{Ac}\right]\cdots \\ -{\text{k}}_{\text{aaac}}\left[{\text{W}}_{\text{p}}\right]\left[\text{Ac}\right]-{\text{k}}_{\text{fphos}}\left[{\text{W}}_{\text{a}}\right]\left[\text{AMPAR}\right]\cdots \\ +\left({\text{k}}_{\text{kbphos}}+{\text{k}}_{\text{kactphos}}\right)\left[{\text{W}}_{\text{a}}\text{AMPAR}\right]\cdots -{\text{k}}_{\text{dephos}}\left[{\text{W}}_{\text{a}}\right]\end{array}$ (7)

其中[W_bAMPAR]、[W_pAMPAR]和[W_aAMPAR]分别是与AMPAR结合的W_b、Wp和Wa的浓度。此外CaMKII亚基可以通过相邻的活性Ac绑定亚基磷酸化生成W_bAc、W_pAc和W_aAc。

与Ac结合的CaMKII亚基的自磷酸化率V_ac是如下：

${\text{V}}_{\text{ac}}=\frac{{\text{K}}_{\text{ac}}{\left({\text{c}}_{\text{b}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]\right)}^2+\left({\text{c}}_{\text{b}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]\right)\left({\text{c}}_{\text{p}}\left[{\text{W}}_{\text{p}}\text{Ac}\right]\right)+\left({\text{c}}_{\text{b}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]\right)\left({\text{c}}_{\text{a}}\left[{\text{W}}_{\text{a}}\text{Ac}\right]\right)}{{\left[{\text{W}}_{\text{tot}}\right]}^2}$ (8)

${\text{K}}_{\text{ac}}={\text{<msup> K <mo>′</mo> </msup>}}_{\text{a}}\left({\text{aT}}_{\text{ac}}+{\text{bT}}_{\text{ac}}^2+{\text{cT}}_{\text{ac}}^3\right)$ (9)

${\text{T}}_{\text{ac}}=\frac{\left[{\text{W}}_{\text{b}}\text{Ac}]+[{\text{W}}_{\text{p}}\text{Ac}]+[{\text{W}}_{\text{a}}\text{Ac}\right]}{\left[{\text{W}}_{\text{tot}}\right]}$ (10)

W_iAc、W_bAc、W_pAc和W_aAc的计算如下：

$\frac{\text{d}\left[{\text{W}}_{\text{i}}\text{Ac}\right]}{\text{dt}}=-{\text{k}}_{\text{iaci}}\left[{\text{W}}_{\text{i}}\text{Ac}\right]+{\text{k}}_{\text{iiac}}\left[{\text{W}}_{\text{i}}\right]\left[\text{Ac}\right]+{\text{k}}_{\text{baciac}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]\mathrm{...}-{\text{k}}_{\text{iacbac}}\left[{\text{W}}_{\text{i}}\text{Ac}\right]\left[{\text{Ca}}_4\text{CaM}\right]$ (11)

$\begin{array}{c}\frac{\text{d}\left[{\text{W}}_{\text{b}}\text{Ac}\right]}{\text{dt}}=-{\text{V}}_{\text{ac}}\left[{\text{W}}_{\text{tot}}\right]+{\text{k}}_{\text{iacbac}}\left[{\text{W}}_{\text{i}}\text{Ac}\right]\left[{\text{Ca}}_4\text{CaM}\right]+{\text{k}}_{\text{bbac}}\left[{\text{W}}_{\text{b}}\right]\left[\text{Ac}\right]\cdots \\ -{\text{k}}_{\text{baciac}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]-{\text{k}}_{\text{bacb}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]\end{array}$ (12)

$\begin{array}{c}\frac{\text{d}\left[{\text{W}}_{\text{p}}\text{Ac}\right]}{\text{dt}}=-{\text{V}}_{\text{ac}}\left[{\text{W}}_{\text{tot}}\right]+{\text{k}}_{\text{aacpa}}\left[{\text{W}}_{\text{p}}\text{Ac}\right]\left[{\text{Ca}}_4\text{CaM}\right]+{\text{k}}_{\text{ppac}}\left[{\text{W}}_{\text{p}}\right]\left[\text{Ac}\right]\cdots \\ -{\text{k}}_{\text{pacaac}}\left[{\text{W}}_{\text{p}}\text{Ac}\right]-{\text{k}}_{\text{pacp}}\left[{\text{W}}_{\text{p}}\text{Ac}\right]\end{array}$ (13)

$\frac{\text{d}\left[{\text{W}}_{\text{a}}\text{Ac}\right]}{\text{dt}}={\text{k}}_{\text{aaca}}\left[{\text{W}}_{\text{a}}\text{Ac}\right]+{\text{k}}_{\text{aaac}}\left[{\text{W}}_{\text{a}}\right]\left[\text{Ac}\right]-{\text{k}}_{\text{pacaac}}\left[{\text{W}}_{\text{a}}\text{Ac}]-[{\text{W}}_{\text{a}}\text{Ac}\right]\left[{\text{Ca}}_4\text{CaM}\right]$ (14)

Ca₄CaM除了能够活化CaMKII之外，还具有与PP2B的非活性形式(PP2B_i)相结合的能力。一旦结合，原本不活跃的磷酸酶随即被激活并转化为其活性形式PP2B_ac。PP2B_ac对于AMPA受体的去磷酸化起介导作用。PP2B_i和PP2B_ac浓度的时间演化为

$\frac{\text{d}\left[\text{PP}2{\text{B}}_{\text{i}}\right]}{\text{dt}}=-{\text{k}}_{\text{ppia}}\left[\text{PP}2{\text{B}}_{\text{i}}\right]\left[{\text{Ca}}_4\text{CaM}\right]+{\text{k}}_{\text{ppai}}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]$ (15)

$\begin{array}{c}\frac{\text{d}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]}{\text{dt}}={\text{k}}_{\text{ppia}}\left[\text{PP}2{\text{B}}_{\text{i}}\right]\left[{\text{Ca}}_4\text{CaM}\right]-{\text{k}}_{\text{ppai}}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]-{\text{k}}_{{\text{f}}_{\text{dephos}}}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]\left[\text{AMPARP}\right]\cdots \\ +\left({\text{k}}_{{\text{b}}_{\text{dephos}}}+{\text{k}}_{{\text{act}}_{\text{dephos}}}\right)\left[\text{PP}2{\text{B}}_{\text{ac}}\text{AMPARP}\right]\end{array}$ (16)

其中[AMPARP]表示与磷酸化AMPA受体(AMPARP)结合的PP2Bac的浓度。AMPA受体(AMPAR)的计算如下：

$\begin{array}{c}\frac{\text{d}\left[\text{AMPAR}\right]}{\text{dt}}=-{\text{k}}_{{\text{f}}_{\text{phos}}}\left[{\text{W}}_{\text{b}}\right]\left[\text{AMPAR}\right]+{\text{k}}_{{\text{b}}_{\text{phos}}}\left[{\text{W}}_{\text{b}}\text{AMPAR}\right]-{\text{k}}_{{\text{f}}_{\text{phos}}}\left[{\text{W}}_{\text{p}}\right]\left[\text{AMPAR}\right]\cdots \\ +{\text{k}}_{{\text{b}}_{\text{phos}}}\left[{\text{W}}_{\text{p}}\text{AMPAR}\right]-{\text{k}}_{{\text{f}}_{\text{phos}}}\left[{\text{W}}_{\text{a}}\right]\left[\text{AMPAR}\right]+{\text{k}}_{{\text{b}}_{\text{phos}}}\left[{\text{W}}_{\text{a}}\text{AMPAR}\right]\cdots \\ +{\text{k}}_{{\text{cat}}_{\text{dephos}}}\left[\text{PP}2{\text{B}}_{\text{ac}}\text{AMPARP}\right]\end{array}$ (17)

其中W_bAMPARP、W_bAMPARP、W_bAMPARP和PP2B_acAMPARP的浓度变化率计算如下：

$\frac{\text{d}\left[{\text{W}}_{\text{b}}\text{AMPAR}\right]}{\text{dt}}={\text{k}}_{{\text{f}}_{\text{phos}}}\left[{\text{W}}_{\text{b}}\right]\left[\text{AMPAR}\right]-\left({\text{k}}_{{\text{b}}_{\text{phos}}}+{\text{k}}_{{\text{cat}}_{\text{phos}}}\right)\left[{\text{W}}_{\text{b}}\text{AMPAR}\right]$ (18)

$\frac{\text{d}\left[{\text{W}}_{\text{p}}\text{AMPAR}\right]}{\text{dt}}={\text{k}}_{{\text{f}}_{\text{phos}}}\left[{\text{W}}_{\text{p}}\right]\left[\text{AMPAR}\right]-\left({\text{k}}_{{\text{b}}_{\text{phos}}}+{\text{k}}_{{\text{cat}}_{\text{phos}}}\right)\left[{\text{W}}_{\text{p}}\text{AMPAR}\right]$ (19)

$\frac{\text{d}\left[{\text{W}}_{\text{a}}\text{AMPAR}\right]}{\text{dt}}={\text{k}}_{{\text{f}}_{\text{phos}}}\left[{\text{W}}_{\text{a}}\right]\left[\text{AMPAR}\right]-\left({\text{k}}_{{\text{b}}_{\text{phos}}}+{\text{k}}_{{\text{cat}}_{\text{phos}}}\right)\left[{\text{W}}_{\text{a}}\text{AMPAR}\right]$ (20)

$\frac{\text{d}\left[\text{PP}2{\text{B}}_{\text{ac}}\text{AMPARP}\right]}{\text{dt}}={\text{k}}_{{\text{f}}_{\text{dephos}}}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]\left[\text{AMPARP}\right]-\left({\text{k}}_{{\text{b}}_{\text{dephos}}}+{\text{k}}_{{\text{cat}}_{\text{dephos}}}\right)\left[\text{PP}2{\text{B}}_{\text{ac}}\text{AMPARP}\right]$ (21)

AMPARP的演化计算如下：

$\begin{array}{c}\frac{\text{d}\left[\text{AMPARP}\right]}{\text{dt}}={\text{k}}_{{\text{cat}}_{\text{phos}}}\left[{\text{W}}_{\text{b}}\text{AMPARP}\right]+{\text{k}}_{{\text{cat}}_{\text{phos}}}\left[{\text{W}}_{\text{p}}\text{AMPARP}\right]+{\text{k}}_{{\text{cat}}_{\text{phos}}}\left[{\text{W}}_{\text{a}}\text{AMPARP}\right]\cdots \\ -{\text{k}}_{{\text{f}}_{\text{dephos}}}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]\left[\text{AMPARP}\right]+{\text{k}}_{{\text{b}}_{\text{dephos}}}\left[\text{PP}2{\text{B}}_{\text{ac}}\text{AMPARP}\right]\end{array}$ (22)

长时程抑制和长时程增强诱导的协议涉及钙信号刺激，钙的浓度变化如下：

$\frac{\text{d}\left[\text{Ca}\right]}{\text{dt}}=-4{\text{k}}_{\text{on}}{\left[\text{Ca}\right]}^4\left[\text{CaM}\right]+4{\text{k}}_{\text{off}}\left[{\text{Ca}}_4\text{CaM}\right]+\text{φ}\left(\text{t}\right)-\text{κ}\left(\left[\text{Ca}\right]-\left[{\text{Ca}}_{\min }\right]\right)$ (23)

$\varphi \left(\text{t}\right)-\text{κ}\left(\left[\text{Ca}\right]-\left[{\text{Ca}}_{\min }\right]\right)$ 描述了钙动力学模型的动态，其中$\varphi \left(\text{t}\right)$ 表示每个最小时间(ms)的钙浓度输入，其值来自输入表，$\text{κ}$ 反映钙扩散、泵送、交换的转移率的术语，[Ca_min]是基底钙浓度。CaM的时间演化为：

$\frac{\text{d}\left[\text{CaM}\right]}{\text{dt}}=-{\text{k}}_{\text{on}}{\left[\text{Ca}\right]}^4\left[\text{CaM}\right]+{\text{k}}_{\text{off}}\left[{\text{Ca}}_4\text{CaM}\right]$ (24)

四个钙离子与CaM的结合生成Ca₄CaM，它的作用是激活PP2B、W_b、W_p、W_bAc和W_pAc。表示Ca₄CaM浓度演化的方程为

$\begin{array}{c}\frac{\text{d}\left[\text{Ca}4\text{CaM}\right]}{\text{dt}}={\text{k}}_{\text{on}}{\left[\text{Ca}\right]}^4\left[\text{CaM}\right]-{\text{k}}_{\text{off}}\left[{\text{Ca}}_4\text{CaM}\right]+{\text{k}}_{\text{ppai}}\left[\text{PP}2{\text{B}}_{\text{ac}}\right]\cdots \\ -{\text{k}}_{\text{ib}}\left[{\text{W}}_{\text{i}}\right]\left[{\text{Ca}}_4\text{CaM}\right]-{\text{k}}_{\text{ppia}}\left[\text{PP}2{\text{B}}_{\text{i}}\right]\left[{\text{Ca}}_4\text{CaM}\right]+{\text{k}}_{\text{bi}}\left[{\text{W}}_{\text{b}}\right]+{\text{k}}_{\text{pa}}\left[{\text{W}}_{\text{p}}\right]\cdots \\ -{\text{k}}_{\text{ap}}\left[{\text{W}}_{\text{a}}\right]\left[{\text{Ca}}_4\text{CaM}\right]-{\text{k}}_{\text{iacbac}}\left[{\text{W}}_{\text{i}}\text{Ac}\right]\left[{\text{Ca}}_4\text{CaM}\right]+{\text{k}}_{\text{baciac}}\left[{\text{W}}_{\text{b}}\text{Ac}\right]\cdots \\ +{\text{k}}_{\text{pacaac}}[{\text{W}}_{\text{p}}\text{Ac}\left]-{\text{k}}_{\text{aacpac}}\right[[{\text{W}}_{\text{a}}\text{Ac}]\left[{\text{Ca}}_4\text{CaM}\right]\end{array}$ (25)

训练后每个PF-PC突触权重的计算如下：

${\text{w}}_{\text{n}}\left(\text{Per}\right)={\text{w}}_0\left(\text{Per}\right)+{\text{w}}_0\left(\text{Per}\right)\text{*}\left(\left[{\text{AMPAR}}_{\text{n}}]-[{\text{AMPAR}}_0\right]\right)$ (26)

其中w_n(Per)表示第Per个周期下第n个突触的权重，[AMAPR_n]表示第n个周期结束时AMPAR的浓度。

Table 1. Relevant kinetic parameters

表1. 相关动力学参数

符号	意义	值
Wtot	非转基因小鼠的CaMKII总浓度	26uM
Ac	F-肌动蛋白总浓度	10 uM
Camin	Ca基础浓度	0.045 uM
CaM	CaM初始浓度	36 uM

续表

AMPAR	非磷酸化AMPAR初始浓度	0.5 uM
AMPARP	磷酸化AMPAR初始浓度	0.5 uM
PP2Bi	无活性PP2B初始浓度	26 uM
${\text{<msup> K <mo>′</mo> </msup>}}_{\text{a}}$	CaMKII自磷酸化速率	0.29/s
cb	b亚基型CaMKII激活系数	75%
cp	p亚基型CaMKII激活系数	100%
ca	a亚基型CaMKII激活系数	80%
a	拟合参数a	0.5
b	拟合参数b	1.956
c	拟合参数c	−1.8
κ	Ca移除速率	4000/s
kbi	Ca4CaM与Wb的分离速率	0.2/s
kib	Ca4CaM与Wi的结合速率	10/(uM*s)
kap	Ca4CaM与Wa的结合速率	10/(uM*s)
kpa	Ca4CaM与Wp的分离速率	0.004/s
kiacbac	Ca4CaM与WiAC的结合速率	10/(uM*s)
kbaciac	Ca4CaM与WbAC的分离速率	1/s
kaacpac	Ca4CaM与WaAC的结合速率	10/(uM*s)
kpacaac	Ca4CaM与WpAC的分离速率	0.02/s
kdephos	CaMKII衰退速率	0.0005/s
kiiac	F-肌动蛋白与Wi的结合速率	10/(uM*s)
kbbac	F-肌动蛋白与Wb的结合速率	10/(uM*s)
kkppac	F-肌动蛋白与Wp的结合速率	10/(uM*s)
kaaac	F-肌动蛋白与Wa的结合速率	10/(uM*s)
kiaci	F-肌动蛋白与WiAC的分离速率	30.1/s
kbacb	F-肌动蛋白与WbAC的分离速率	150.5/s
kpacp	F-肌动蛋白与WpAC的分离速率	1505/s
kaaca	F-肌动蛋白与WaAC的分离速率	301/s
kppia	Ca4CaM与PP2Bi的结合速率	0.15/(uM*s)
kppai	Ca4CaM与PP2Bac的分离速率	0.00042/s
kon	Ca与CaM的结合速率	2000/(uM^4*s)

续表

koff	Ca与CaM的分离速率	2.3 * 10^6/s
${\text{k}}_{{\text{f}}_{\text{phos}}}$	CaMKII磷酸化AMPAR的速率	0.5/(uM*s)
${\text{k}}_{{\text{b}}_{\text{phos}}}$		72.283/s
${\text{k}}_{{\text{cat}}_{\text{phos}}}$		6/s
${\text{k}}_{{\text{f}}_{\text{dephos}}}$	PP2B去磷酸化AMPAR的速率	0.5/(uM*s)
${\text{k}}_{{\text{b}}_{\text{dephos}}}$		72.283/s
${\text{k}}_{{\text{cat}}_{\text{dephos}}}$		6/s

2.3. RNN与ODE对比

钙级联动力学系统的求解涉及到常微分方程组(Ordinary Differential Equation, ODE)的计算，对于无法求出解析解的情况，通常人们使用欧拉法或龙格–库塔法求出数值解。ODE的结构如下：

$\text{<msup> y <mo>′</mo> </msup>}\left(\text{t}\right)=\text{f}\left(\text{y}\left(\text{t}\right),\text{t}\right)$ (27)

欧拉法的算式如下：

$\text{y}\left(\text{t}+\text{h}\right)=\text{y}\left(\text{t}\right)+\text{h<msup> y <mo>′</mo> </msup>}\left(\text{t}\right)$ (28)

其中h是步长，步长越短求解越精确，但同时有耗时的问题。有研究指出[11]，ODE的欧拉解法实际上是循环神经网络(Recurrent Neural Network, RNN)的一个特例。在RNN中对于一个向量序列的输入$\left({\text{x}}_1,{\text{x}}_2,{\text{x}}_3,\cdots ,{\text{x}}_{\text{t}}\right)$ ，输出时另一个序列$({\text{y}}_1,{\text{y}}_2,{\text{y}}_3,\cdots ,{\text{y}}_{\text{t}})$ ，满足以下递归关系：

${\text{y}}_{\text{t}}=\text{f}\left({\text{y}}_{\text{t}-1},{\text{x}}_{\text{t}},\text{t}\right)$ (29)

式(28)中的h是离散变量，式(4)~(29)中的t是连续变量。在式(4)~(28)中以h为时间单位，记$t=nh$ ，n是整数，那么式(4)~(28)变成：

$\text{y}\left(\left(\text{n}+1\right)\text{h}\right)=\text{y}\left(\text{nh}\right)+\text{hf}\left(\text{y}\left(\text{nh}\right),\text{nh}\right)$ (30)

此时(30)满足(29)的递归关系。实际上前向传播就是解ODE(RNN的预测过程)，反向传播是推断ODE的参数(RNN的训练过程)。

3. 实验与分析

3.1. 实验方案

以视动反应 (optokinetic eye movement response, OKR)实验作为例子，用该模型进行仿真，观察视动反应的运动学习记忆形成过程。刺激信号如图2所示，一次训练周期有2000 ms，刺激信号和误差信号的脉冲数量分布符合正弦信号的频率，分别是30 hz和3 hz，最后两者均在2000 ms时结束，这样的流程对应一次训练，共训练200周期。

3.2. 试验分析

为了方便观察AMPAR随训练次数和钙尖峰的变化，此处以延迟眨眼反射范式为例说明，周期仍为2000 ms，输入和误差信号采用固定频率。输入信号从200 ms开始，1200 ms结束，误差信号从450 ms开始，持续30 ms。用范围在0到30 uM内的多个钙尖峰分别计算200周期后AMPAR的含量，通过插值得出不同训练次数和不同钙尖峰AMPAR的变化。

Figure 2. Stimulus signal of optokinetic eye movement response

图2. 视动反射刺激信号

如图3所示，上部分图不同曲线代表了不同训练次数下AMPAR浓度随钙尖峰的变化，下部分图是不同钙尖峰下AMPAR浓度随训练次数的变化。上部分图中钙尖峰0到2 uM对AMAPR浓度几乎没有影响，因此认为该区间无可塑性；从2 uM开始，训练30次、100次和200次后的AMAPR浓度都远远高于初始值，随着训练次数上升，AMPAR浓度回到初始值时对应的钙尖峰从12.2 uM逐渐下降到9.1 uM，这几个从2 uM到不同钙尖峰结束值的区间分别是这几个不同训练次数曲线的LTP区间，然而训练5次后的曲线几乎无LTP。LTP区间后AMPAR浓度低于初始值，是LTD区间。因此可得出结论，在学习效果突显后，随着训练次数增加，可塑性方向转换的阈值逐渐下降。

下部分图中钙尖峰在2 uM和6.3 uM时，AMPAR浓度随着训练次数上升。而从8.3 uM开始，浓度有下降趋势，到了10 uM，浓度在训练110次后低于初始值。之后随着钙尖峰上升，可塑性方向转换所需的训练次数越来越少。

Figure 3. Changes of AMPAR with training times under different calcium spikes

图3. AMPAR在不同钙尖峰下随训练次数的变化

对于视动反应，本文使用优化了RNN结构的seq2seq进行预测。首先通过TorchDiff库的odeint求解器使用龙格–库塔法求解在第11次到第18次训练内16个平行纤维上的突触权重对应的钙级联动力学系统18种物质的量的变化，因此数据集时间长度是16,000 ms，有16个样本和18种特征。通过小脑神经元模型生成随机的时空活动序列和误差信号，计算出钙尖峰序列。之后取对应样本在第1次到第9次的状态变化作为seq2seq的训练集，在预处理阶段，由于所有特征均大于等于0，且部分特征的数值数量级差异大，这部分特征先做对数化处理再进行归一化。外部输入通过钙尖峰序列查表得出每毫秒对应的钙信号。共训练1000次，loss在0.00376处收敛。

图4从左到右分别是三条平行纤维对应的AMPARP、AMPAR和PP2B_i变化，其中AMPARP和AMPAR的含量用百分比表示，PP2B的量化单位是微摩尔(uM)，时间统一是8次训练总时长。其中同一种颜色的曲线表示来自同一条PF，散点是对应曲线的预测值。

Figure 4. Changes of chemical molecular content with training

图4. 化学分子含量随训练的变化

设定AMPARP和AMPAR起始值都是50%，随着训练次数AMPARP含量上升，AMPAR下降，PP2B_i呈上升趋势。红色样本的物质变化率较大，原因是这些样本来自不同的颗粒细胞簇，活跃细胞簇对应的动力学系统对外来的钙尖峰更为敏感。

AMPARP和AMPAR在钙级联中互相转化，因此当误差信号紧接着输入信号时，会导致内质网的IP₃受体门控通道的打开阈值暂时降低，在这个短暂的时间窗口由误差信号驱动释放在浦肯野细胞树突棘上的Ca²⁺，能够与CaM相结合，形成Ca₄CaM。接下来Ca4CaM参与CaMKII自磷酸化从而导致AMPAR的磷酸化，形成LTD。这时Ca₄CaM几乎没有参与由PP2B_i转化成PP2Bac介导的AMPAR去磷酸化，LTP受到抑制，PP2B_i含量上升。

4. 讨论

使用TorchDiff库或Scipy库的微分方程求解器时，由于钙级联动力学系统状态方程的复杂性、高耦合度，随着要计算的PF数增加，计算时长呈指数上升。求解器每算出1 ms的后状态量，需要更改作为外部输入的钙信号值才能计算下一毫秒的状态，另外方程组对精度的要求非常高，自适应步长求解器在每次计算下一毫秒的状态时的步长多达80,000，除了迭代容易出现异常值外，还使得计算10,000条PF在一毫秒的耗时异常漫长。在8G显存的NVIDIA显卡使用Torch框架进行运算，计算一毫秒后状态耗时6分钟，若以训练200个周期为参考，耗时约167天。对于同样运算条件下的seq2seq网络，每个细胞状态代表1 ms的时长，通过GPU和多Stream的使用提速训练200个连续周期，耗时缩短到18小时，在预测少量连续周期时保持较高的精准度。

5. 结语

本文介绍了PF-PF钙级联通路以及相关动力学的通信机制。长时程增强和长时程抑制的背后是AMPAR的磷酸化和去磷酸化，其受到浦肯野细胞树突棘中钙离子含量影响。IP₃R门控通道接收误差信号和输入信号分别诱导的钙刺激和IP₃是打开内质网钙库的关键。

在动力学方程组的计算上对比ODE和RNN的相似性。RNN的前向传播与解ODE的过程相同，RNN的反向传播是推断ODE的参数然后。用seq2seq对动力学系统做短时预测，分析部分状态量的变化，系统的高耦合度是状态量此消彼长的原因。

最后讨论两种计算方式的优缺点，龙格–库塔法相对精准但耗时长，seq2seq网络的耗时短但只能预测短连续周期的时间序列。