\[ \newcommand{\bs}{\boldsymbol} \newcommand{\bsX}{\boldsymbol{X}} \newcommand{\bf}{\mathbf} \newcommand{\msc}{\mathscr} \newcommand{\mca}{\mathcal} \newcommand{\T}{\text{T}} \newcommand{\rme}{\mathrm{e}} \newcommand{\rmi}{\mathrm{i}} \newcommand{\rmj}{\mathrm{j}} \newcommand{\rmd}{\mathrm{d}} \newcommand{\rmm}{\mathrm{m}} \newcommand{\rmb}{\mathrm{b}} \newcommand{\and}{\land} \newcommand{\or}{\lor} \newcommand{\exist}{\exists} \newcommand{\sube}{\subseteq} \newcommand{\lr}[3]{\left#1 #2 \right#3} \newcommand{\intfy}{\int_{-\infty}^{+\infty}} \newcommand{\sumfy}[1]{\sum_{#1=-\infty}^{+\infty}} \newcommand{\vt}{\vartheta} \newcommand{\ve}{\varepsilon} \newcommand{\vp}{\varphi} \newcommand{\Var}{\text{Var}} \newcommand{\Cov}{\text{Cov}} \newcommand{\edef}{\xlongequal{def}} \newcommand{\prob}{\text{P}} \newcommand{\Exp}{\text{E}} \newcommand{\t}[1]{\text#1} \newcommand{\N}{\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \newcommand{\versionofnewcommand}{\text{260125}} \]

Dendrites and Synapses

Synapses

Basic Concepts of Synapses

Activation of a pre-synaptic neuron results in a release of neurotransmitters into the synaptic cleft. After the diffusion in the cleft, transmitters activate receptors located in the post-synaptic membrane.

receptors in post-synaptic membrane can be classified into:

​ ionotropic receptors: have a direct influence on the state of an associated ion channel.

​ metabotropic receptors: control the state of an ion channel through signal pathways such as G protein and related certain second messengers.

Both cases lead to the opening of certain ion channel, which makes a current, including:

​ EPSC: excitatory post-synaptic transmembrane current, and

​ IPSC: inhibitory post-synaptic transmembrane current

Mathematical Modeling of Ion Channels in Synapse Membrane

神经递质在突触后膜作用的影响, 可以被视为是后膜电导 \(g_{syn}(t)\) 的变化. 进而, 在突触后膜:

\[ I_{\text{syn}}(t)=g_{\text{syn}}(t)(u(t)-E_{\text{syn}}) \]

其中, \(E_{syn}\) 就是 Reversal Potential / Equilibrium Potential. \(u(t)\) 是膜电位. 对于抑制性突触, \(E_{\text{syn}}\) 常设为 \(-75 \text{mV}\); 兴奋性突触的 \(E_{\text{syn}}\) 则约为 \(0\).

假设忽略神经递质在突触间隙扩散所需的时间, 那么, 就可以将动作电位传递到突触的时间记为 \(t^f\). 假设神经递质是量子释放的, 即每次动作电位引起的神经递质释放量都相等, 且后膜电导变化的响应随时间呈指数衰减, 则:

\[ g_{\text{syn}}(t)=\sum_f \bar{g}_{\text{syn}}{e}^{-(t-t^f)/\tau} \]

但这样还不够! 由于突触后膜无法对未来时刻到达的神经递质响应, 必须乘上一个单位阶跃函数保证因果性.

\[ g_{\text{syn}}(t)=\sum_f \bar{g}_{\text{syn}}{e}^{-(t-t^f)/\tau}\cdot\Theta(t-t^f) \]

再加上一个 smooth rise 项, 刻画电导升高的过程:

\[ g_{\text{syn}}(t)=\sum_f \bar{g}_{\text{syn}}\left[1-e^{-(t-t^f)/\tau_{\text{rise}}}\right]{e}^{-(t-t^f)/\tau}\Theta(t-t^f) \]

根据实验观察, post-synaptic current 由两部分构成, 一部分衰减很快, 另一部分衰减很慢. 因此, 将指数衰减项写作两个电流的加权和:

\[ g_{\text{syn}}(t)=\sum_f \bar{g}_{\text{syn}}\left[1-e^{-(t-t^f)/\tau_{\text{rise}}}\right]\left[a{e}^{-(t-t^f)/\tau_{\text{fast}}}+(1-a){e}^{-(t-t^f)/\tau_{\text{slow}}}\right]\Theta(t-t^f) \]

我对此的认知是, 建模的精细化是可以无穷尽的, 不过不是只有精细化的建模才可称为智慧, 如何用简单的形式近似也是智慧.

基于离子通道开放比例的模型 (Destexhe et al.)

接下来介绍一个相对 Biological 的看待 \(g_{\text{syn}}(t)\) 的角度: 假设突触后膜电导的上限是 \(g_{\text{max}}\), 离子通道开放的比例为 \(R(t)\), 那么就有:

\[ g_{\text{syn}}(t)=g_{\text{max}}R(t) \]

假设单位时间内有比例为 \(\alpha\) 的空载受体结合了神经递质; 而神经递质与受体结合后, 单位时间内有比例为 \(\beta\) 的递质解离, 导致通道重新关闭.

\[ \text{R+T}\underset{\beta}{\overset{\alpha}{\rightleftharpoons}}\text{RT} \]

则, 考察 \(\text{d}R/\text{d}t\), 应当有:

\[ \frac{\text{d}R}{\text{d}t}=\alpha(t)(1-R)-\beta R \]

其中, \(\alpha\) 应当和神经递质在受体周围的浓度有关, 因此可以当作是时间的函数 \(\alpha(t)\). \(\beta\) 则一般是一个常数. 为了使形式更加简明, 我们把 \(\alpha(t)\) 写作 \(\alpha N(t)\). 其中 \(\alpha\) 变成常数, 和 \(\beta\) 的形式对应; \(N(t)\) 则由下式给出:

\[ N(t)=\int_{0}^{\infty}\gamma(s)S_j(t-s)\text{d}s \]

其中, \(\gamma(t)\) 与神经递质在受体周围的浓度随时间的变化有关, \(S_j(t-s)\) 则是先前提到过的突触 \(j\) 的 spike train. (Destexhe et al. 在 1994 年的文章中说, \(\gamma(t)\) 就是神经递质的浓度)

那么, \(\text{d}R/\text{d}t\) 的动力学方程可以写作:

\[ \frac{\text{d}R}{\text{d}t}=\alpha N(t)(1-R)-\beta R \]

Destexhe 的文章里提到:

There is evidence from both the neuromuscular junction (Anderson and Stevens 1973) and excitatory central synapses (Colquhoun et al. 1992) that the concentration of transmitter in the cleft rises and falls very rapidly.

因此, 几乎只需要考虑 \(N(t)=N_{\text{max}}\) 和 \(N(t)=0\) 的情况.

\(N(t)=N_{\max}\) 时:

\[ \begin{aligned} \frac{\text{d}R}{\text{d}t} &=\alpha N_{\max}(1-R)-\beta R\\ &=-(\alpha N_{\max}+\beta)R+\alpha N_{\max}\\ \end{aligned} \]

这是个简单的一阶线性常微分方程, 其解为:

\[ R(t)=\left[R(t_0)-\frac{\alpha N_{\max}}{\alpha N_{\max}+\beta}\right]e^{-(\alpha N_{\max}+\beta)(t-t_0)}+\frac{\alpha N_{\max}}{\alpha N_{\max}+\beta} \]

设

\[ \begin{aligned} R(\infty) &=\frac{\alpha N_{\max}}{\alpha N_{\max}+\beta}\\ \frac{1}{\tau_R} &=\alpha N_\max+\beta \end{aligned} \]

则:

\[ R(t)=\left[R(t_0)-R(\infty)\right]e^{-{(t-t_0)}/{\tau_R}}+R(\infty) \]

Inhibitory Synapses

The effect of fast inhibitory neurons in the central nervous system of higher vertebrates is almost exclusively conveyed by a neurotransmitter called γ-aminobutyric acid (GABA).

inhibitory synapses 后膜上的 GABA 受体可以分为两类: \(\text{GABA}_{\text{A}}\) 和 \(\text{GABA}_{\text{B}}\). 其中: \(\text{GABA}_{\text{A}}\) 本身就作为 \(\text{Cl}^-\) 离子通道, reversal potential \(E_{\text{syn}}\approx -75\text{mV}\). 而 \(\text{GABA}_\text{B}\) 则作为 G 蛋白偶联受体激活 \(\text{K}^+\) 通道 (GIRK, G protein-coupled inwardly rectifying potassium channels), \(E_{\text{syn}}\approx -100\text{mV}\).

根据受体种类不同, inhibitory synapses 可以相应地分为 \(\text{GABA}_{\text{A}}\) 和 \(\text{GABA}_{\text{B}}\) 突触. 两突触的相关参数如下:

	\(\text{GABA}_{\text{A}}\)	\(\text{GABA}_{\text{B}}\)
\(\tau_{\text{rise}}\)	\(\approx 1\ \text{ms}\)	\(25-50\ \text{ms}\)
\(a\)	\(1\)	\(0.8\)
\(\tau_{\text{fast}}\)	\(\approx 6\ \text{ms}\)	\(100-300\ \text{ms}\)
\(\tau_{\text{slow}}\)	-	\(500-1000\ \text{ms}\)

Excitatory Synapses

Excitatory synapses 更多地用 glutamate 作为神经递质. 突触后膜上的受体同样有许多类别, 这些分类一般是基于对何种药物敏感. 例如, NMDA-sensitive 的受体就是 NMDA receptors; AMPA-sensitive 的受体就是 AMPA 受体.

AMPA (2-Amino-3-(3-hydroxy-5-methyl-isoxazol-4-yl)propanoic acid)

AMPA receptor channel 允许 \(\text{Na}^+\) 和 \(\text{K}^+\) 通过, 而 \(\text{Ca}^{2+}\) 不可以. 它对递质的响应非常快, \(\tau_{\text{rise}}\approx 0\), 因此不必乘上 smooth rise 项. 另外, decay time constant \(\tau_{\text{fast}}\) 约为 \(2-5\ \text{ms}\), 也很快, 并且没有 \(\tau_{\text{slow}}\) 项.

NMDA (N-methyl-D-aspartate)

NMDA receptor 激活与否不光与递质的结合有关, 还与膜电位有关. 其机制为: 在静息电位时, 细胞膜外的\(\text{Mg}^{2+}\) 会抑制离子通道活性; 膜电位约为 \(-50\ \text{mV}\) 时, \(\text{Mg}^{2+}\) 的抑制作用被移除, 离子通道持续开放约 \(10-100\ \text{ms}\). \(\text{Mg}^{2+}\) 的效应可以这样刻画:

\[ g_{\infty}(u,[\text{Mg}^{2+}]_\text{o})=(1+\beta\ \text{e}^{\alpha u}[\text{Mg}^{2+}])^{-1} \]

把它乘进 \(g_{\text{NMDA}}(t)\) 的表达式中, 就得到:

\[ g_{\text{NMDA}}(t)=\bar{g}_{\text{NMDA}}(t)\left[1-e^{-(t-t^f)/\tau_{\text{rise}}}\right]{e}^{-(t-t^f)/\tau_{\text{decay}}}g_{\infty}(u,[\text{Mg}^{2+}]_\text{o})\Theta(t-t^f) \]

其中, \(\tau_{\text{rise}}\) 约为 \(3-15\ \text{ms}\), \(\tau_{\text{decay}}\) 约为 \(40-100\ \text{ms}\). 无论是 rise 还是 decay, 都显著地慢于 AMPA 受体. 这意味着 NMDA 受体可以作为一个低通滤波器, 或是信号整合元件.

另外, NMDA 受体对 \(\text{Ca}^{2+}\) 的 permeability 要比对 \(\text{Na}^+\) 和 \(\text{K}^+\) 高得多. 具体有什么作用, 详见相关的生物学研究吧, 不在这里赘述了.

Rapid Synaptic Dynamics

众所周知, 突触之间的联系是有记忆的. Long-term synaptic plasticity 会在 Chapter 19 给出, 这里先介绍一个随时间迅速衰减的 short-term synaptic plasticity.

短时的突触可塑性 (synaptic plasticity) 可以用 \(\bar{g}_{\text{syn}}\) 的变化来表达. \(\bar{g}_{\text{syn}}\) depends on the fraction \(P_{\text{rel}}\) of presynaptic sites releasing neurotransmitter. Facilitation and depression can both be modeled as presynaptic processes that modify \(P_{\text{rel}}\). With each presynaptic spike, the number of available presynaptic release sites changes. 假设 \(P_{\text{rel}}\) 以时间参数 \(\tau_P\) 衰减到 \(P_0\), 那么,

\[ \begin{aligned} &\frac{\text{d}P_{\text{rel}}}{\text{d}t}=-\frac{P_{\text{rel}}-P_0}{\tau_P}+f_F(1-P_\text{rel})\sum_f \delta(t-t^f)\\ &\frac{\text{d}P_{\text{rel}}}{\text{d}t}=-\frac{P_{\text{rel}}-P_0}{\tau_P}-f_DP_\text{rel}\sum_f \delta(t-t^f) \end{aligned} \]

where \(f_F\) denotes the degree of facilitation and \(f_D\) denotes depression.

Dendritic Tree

对于一个神经元, 它从树突接受来自上一神经元的信号, 而后在轴丘 (axon hillock) 产生动作电位. 从树突到轴丘, 电压以怎样的形式传导? 这一部分将会对该问题给予一些解释.

Derivation of the Cable Equation

在树突上选取一小段 \(\mathrm{d}x\), 设该段的电容为 \(c\ \rmd x\), 注入细胞的电流 \(i_{\text{ext}}\ \rmd x\), 该部分离子通道电流为 \(i_{\text{ion}}\ \rmd x\); 根据基尔霍夫定律:

\[ i(t,x\text+{d}x)-i(t,x)=c\ \rmd x\cdot\frac{\partial}{\partial t}u(t,x)+ \sum_{\text{ion}}i_\text{ion}(t,x)\ \rmd x-i_{\text{ext}}(t,x)\ \rmd x \]

假设该段的电阻为 \(r_L\ \rmd x\), 则对应的电压:

\[ u(t,x+\rmd)-u(t,x)=r_L\ i(t,x)\ \rmd x \]

上面两个式子, 两边同时除以 \(\rmd x\):

\[ \begin{aligned} &\frac{\partial}{\partial x}u(t,x)=r_L\ i(t,x)\\ &\frac{\partial}{\partial x}i(t,x)=c\frac{\partial}{\partial t}u(t,x)+ \sum_{\text{ion}}i_\text{ion}(t,x)-i_{\text{ext}}(t,x) \end{aligned} \]

从而得出 Cable Equation:

\[ \frac{\partial^2}{\partial x^2}u(t,x)=c\ r_L\frac{\partial}{\partial t}u(t,x)+ r_L\sum_{\text{ion}}i_\text{ion}(t,x)-r_L\ i_{\text{ext}}(t,x) \]

Cable Equation for Passive Dendrite

对于 passive dendrite, 只考虑膜的 leaky conductance \(g_l=1/r_T\) 和 leaky reversal potential \(E_l\), 即:

\[ \sum_{\text{ion}}i_{\text{ion}}(t,x)=g_l(u-E_l) \]

Let \(\lambda^2=r_T/r_L\), 作为电缆的空间常数, 在生理学上叫 electrotonic length scale, 电紧张长度; \(\tau=r_T\ c\), 作为膜的时间常数, 那么由 cable equation:

\[ \lambda^2\cdot\frac{\partial^2}{\partial x^2}u(t,x)=\tau\cdot\frac{\partial}{\partial t}u(t,x)+ [u(t,x)-E_l]-r_T\ i_{\text{ext}}(t,x) \]

这样的形式其实还不够简洁优美, 我们进行一个坐标轴尺度的伸缩 (无量纲化):

\[ \begin{aligned} x&\to&\hat{x}&=x/\lambda\\ t&\to&\hat{t}&=t/\tau \end{aligned} \]

以及电流电压的尺度变化 (此时电流单位变为电压单位):

\[ \begin{aligned} i&\to&\hat{i}&=\sqrt{r_T\ r_L}\ i\\ i_{\text{ext}}&\to&\hat{i}_{\text{ext}}&=r_T\ i_{\text{ext}}\\ u&\to& \hat{u}&=u-E_l \end{aligned} \]

省略所有的 hat 上标, 就有:

\[ \frac{\partial^2}{\partial x^2}u(t,x)=\frac{\partial}{\partial t}u(t,x)+u(t,x)-i_{\text{ext}}(t,x) \]

或者移项写作:

\[ \frac{\partial}{\partial t}u(t,x)=\frac{\partial^2}{\partial x^2}u(t,x)-u(t,x)+i_{\text{ext}}(t,x) \]

上式中, 第一项为 diffusion term, 第二项为 decay term, 第三项为 source term.

\[ \left[\frac{\partial}{\partial t}-\frac{\partial^2}{\partial x^2}+1\right]u(x,t)=i_{\text{ext}}(t) \]

Stationary Solutions of the Cable Equation

首先, 只需考虑这样的简单的静态情况:

\[ \frac{\partial}{\partial t}u(t,x)=0 \]

此时,

\[ \frac{\partial^2}{\partial x^2}u(t,x)-u(t,x)=-i_{\text{ext}}(t,x) \]

其通解为:

\[ u(t,x)=c_1\ \sinh(x)+c_2\ \cosh(x) \]

如果给予边界条件 \(u(\pm\infty)=0\), 则:

\[ u(t,x)=\int\frac{1}{2}e^{-|x-x^\prime|}i_{\text{ext}}(x^\prime)\ \rmd x^\prime \]

如果 input 电流是个在 \(x=0\) 处产生的冲激, \(i_{\text{ext}}(t,x)=\delta(x)\), 则:

\[ u(t,x)=\frac{1}{2}e^{-|x|} \]

更复杂的随时间动态的分析将在接下来的小节中给出.

Green's Function of the Passive Cable

Green's Function

为了得出 Cable Equation 对任意 input 的解, 还是先考虑最简单的情况: 如果输入电流 \(i_{\text{ext}}(t,x)\) 是 \(t=0\) 时刻在 \(x=0\) 短脉冲, 即

\[ i_{\text{ext}}(t,x)=\delta(t)\ \delta(x) \]

那么就有:

\[ \frac{\partial}{\partial t}u(t,x)-\frac{\partial^2}{\partial x^2}u(t,x)+u(t,x)=\delta(t)\ \delta(x) \]

接下来, 两边对空间变量 \(x\) 做 Fourier 变换:

首先, 时域上有 \(\omega T=2\pi\), 空间上有 \(k\lambda=2\pi\), 即 \(k\) 是空间上的频率. 因此, 我们这样书写 Fourier 变换的结果.

\[ u(t,\rmj k)=\int_{-\infty}^{+\infty}u(t,x)\rme^{-\rmj kx}\rmd x \]

于是, 相应地就有:

\[ \begin{aligned} &\int_{-\infty}^{+\infty}\left[ \frac{\partial}{\partial t}u(t,x)-\frac{\partial^2}{\partial x^2}u(t,x)+u(t,x) \right]\rme^{-\rmj kx}\rmd x\\ =\ &\frac{\partial}{\partial t}u(t,\rmj k)-(\rmj k)^2u(t,\rmj k)+u(t,\rmj k)\\ =\ &\frac{\partial}{\partial t}u(t,\rmj k)+k^2u(t,\rmj k)+u(t,\rmj k) \end{aligned} \]

这里用到了如下引理:

\[ x(t)\overset{\mathcal{F}}\leftrightarrow X(\rmj\omega)\\ \Rightarrow x^\prime(t) \overset{\mathcal{F}}\leftrightarrow \rmj\omega \cdot X(\rmj\omega) \]

又:

\[ \intfy\delta(t)\ \delta(x)\ \rme^{-\rmj kx}\rmd x=\delta(t)\cdot1 \]

则, 原方程化为如下简单易解的微分方程:

\[ \frac{\partial}{\partial t}u(t,\rmj k)+(1+k^2)u(t,\rmj k)=\delta(t) \]

其解为:

\[ \begin{aligned} u(t,\rmj k)&=C\rme^{-(1+k^2)t}+\rme^{-(1+k^2)t}\int_{-\infty}^{t}\delta(\tau)\ \rme^{(1+k^2)\tau}\rmd \tau\\ &=(C+1)\rme^{-(1+k^2)t}\Theta(t) \end{aligned} \]

为了后续计算的方便, 取常数 \(C=-1+\frac{1}{\sqrt{2\pi}}\), 那么就有:

\[ u(t,\rmj k)=\frac{1}{\sqrt{2\pi}}\exp[-(1+k^2)\ t]\Theta(t) \]

做个 Fourier 逆变换给它变回去! 解得:

\[ u(t,x)=\frac{\Theta(t)}{\sqrt{4\pi t}}\exp\left[ -t-\frac{x^2}{4t} \right]\overset{\text{def}}{=}G_{\infty}(t,x) \]

\(G_{\infty}(t,x)\) 被称为 Green's Function.

Checking the Green's Property

容易计算得:

\[ \left[\frac{\partial}{\partial t}-\frac{\partial^2}{\partial x^2}+1\right]G_\infty(x,t)=\frac{1}{\sqrt{4\pi t}}\exp\left(-t-\frac{x^2}{4t}\right)\delta(t) \]

而 \(t\to 0\) 时, 恰好就有:

\[ \frac{1}{\sqrt{4\pi t}}\exp\left(-t-\frac{x^2}{4t}\right)=\delta(x) \]

Superposition of Green's Functions

General solution 的形式可以由线性叠加得出:

\[ u(t,x)=\int_{-\infty}^{t}\rmd t^\prime \intfy \rmd x^\prime G_\infty(t-t^\prime,x-x^\prime)i_{\text{ext}}(t^\prime, x^\prime) \]

将这个结果代入原方程中, 的确是 make sense 的.

Nonlinear Extensions to the Cable Equation

上述的讨论都是基于 leaky 的离子通道, 没有考虑具有复杂动力学特性的, 类似 HH model 里的离子通道. 把这些东西再加进去的话, 恐怕是只能给出数值解了. 考虑到其消耗算力巨大, 且没有太多理论上的扩展, 在此不做说明.

Axons

对于 unmyelinated (无髓鞘) axons, 用前文介绍的 Cable Equation 就可以描述了.

本节主要讨论 myelinated axons. 假设 nodes of Ranvier 之间的距离都为 \(L\), 第 \(n\) 个 node 处的膜电位 \(u_n\) 可以由基尔霍夫定律导出:

\[ C\frac{\rmd u_n}{\rmd t}=\frac{1}{Lr_L}[(u_{n+1}-u_n)-(u_n-u_{n-1})]-\sum_{\text{ion}}I_{\text{ion,n}}(t) \]

不过为什么在有髓鞘轴突能够进行跳跃式的传播, 我认为书中没有清晰的解释. 有一个直觉是, diffusion term 的空间常数足够大, 以至于只靠电位在轴向的快速 diffusion 就可以保证到达下一郎飞结时仍能激发动作电位.

Compartmental Models

书中这一部分主要介绍了树突的区室化模型, 不过没有详细深入地讲. 介绍的 software tools 也都感觉太老旧了 (这是可以说的吗). 据说有些机器学习的模型参考了这个, 不过真的能 work 吗... 尤其是放在当今的时代.