【brainpy学习笔记】突触模型2——化学突触的生理学模型、电突触模型

Fellyhosn

已于 2023-06-24 14:46:17 修改

阅读量1k

点赞数 5

分类专栏：计算神经科学文章标签：学习笔记 python 神经网络

于 2023-05-05 16:32:12 首次发布

本文链接：https://blog.csdn.net/Fellyhosn/article/details/130492993

版权

计算神经科学专栏收录该内容

11 篇文章 45 订阅

订阅专栏

参考书目：《神经计算建模实战——基于brainpy》吴思

书接上文：

【brainpy学习笔记】突触模型1——化学突触的现象学模型https://blog.csdn.net/Fellyhosn/article/details/130412287?spm=1001.2014.3001.5501

3 化学突触的生理学模型

前文已经提到化学突触的两种建模思路，即根据突触前的神经递质到达突触后神经元产生的电流形状建模和根据突触后膜各个离子通道的动力学性质进行建模两种思路。前者称为现象学模型，后者称为生理学模型。下面进行生理学模型的介绍。

3.1 建模离子通道的开放与关闭

神经元的电导模型一章已经介绍了基本的离子通道建模，即用概率s表示离子通道开放的概率，从宏观上来模拟离子通道的各个闸门打开/关闭的比例。需要区别的是，突触上的离子通道不再是由电压调控而是配体调控的，离子通道从开放到关闭的转换速率依赖于突触间隙中的神经递质的浓度（设为[T]），因此将开放到关闭的转换速率记作 $\alpha[T]$ ；从关闭到开放是由配体的解离导致，与[T]无关，记作 $\beta$ 。 s的微分方程为： $\frac{ds}{dt}=\alpha[T](1-s)-\beta s$

下面探究[T]的取值，[T]在突触前神经元发放动作电位后升高，之后逐渐降为0，若忽略升高和降低的时间，则其变化可以近似于δ函数，即 $[T]=\displaystyle\sum_{t^{(f)}}(t-t^{(f)})$ 。下面将该模型应用到具体的突触模型之中

3.2 AMPA模型

根据AMPA受体的性质，一般采用电导模式对其建模。分别在25/50/75/100/160ms时产生动作电位进行模拟，代码如下：

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as plt

class AMPA(bp.dyn.TwoEndConn):
    def __init__(self,pre,post,conn,g_max=0.02,E=0.,alpha=0.98,beta=0.18,delay_step=0,
                 T_duration=0.5,T=0.5,method='exp_auto',name=None):
        super(AMPA,self).__init__(pre=pre,post=post,conn=conn,name=name)
        self.g_max = g_max
        self.E = E
        self.alpha = alpha
        self.beta = beta
        self.delay_step = delay_step
        self.T_duration = T_duration
        self.T = T
        self.pre2post = self.conn.require('pre2post')
        self.s = bm.Variable(bm.zeros(self.post.num))
        self.g = bm.Variable(bm.zeros(self.post.num))
        self.spike_arrival_time = bm.Variable(bm.ones(self.post.num)*-1e7)
        self.delay = bm.LengthDelay(self.pre.spike,delay_step)
        self.integral = bp.odeint(method=method,f=self.ds)
    def ds(self,s,t,TT):
        return self.alpha*TT*(1-s)-self.beta*s
    def update(self,tdi):
        delayed_pre_spike = self.delay(self.delay_step)
        self.delay.update(self.pre.spike)
        self.spike_arrival_time.value = bm.where(delayed_pre_spike,tdi.t,self.spike_arrival_time)
        TT = ((tdi.t - self.spike_arrival_time) < self.T_duration) * self.T
        self.s.value = self.integral(self.s,tdi.t,TT,tdi.dt)
        self.g.value = self.g_max * self.s
        self.post.input += self.g * (self.E - self.post.V)
def run_syn(syn_model, title, run_duration=200., sp_times=(10, 20, 30), **kwargs):
    neu1 = bp.neurons.SpikeTimeGroup(1, times=sp_times, indices=[0] * len(sp_times))
    neu2 = bp.neurons.HH(1, V_initializer=bp.init.Constant(-70.68))
    syn1 = syn_model(neu1, neu2, conn=bp.connect.All2All(), **kwargs)
    net = bp.dyn.Network(pre=neu1, syn=syn1, post=neu2)

    # 进行模拟
    runner = bp.dyn.DSRunner(net, monitors=['pre.spike', 'post.V', 'syn.g', 'post.input','syn.s'])
    runner.run(run_duration)

    # 可视化
    fig, gs = bp.visualize.get_figure(11, 1, 0.5, 6.)
    ax = fig.add_subplot(gs[0:1, 0])
    plt.plot(runner.mon.ts, runner.mon['pre.spike'], label='pre.spike',color='r')
    plt.legend(loc='upper right')
    plt.title(title)
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['right'].set_visible(False)

    ax = fig.add_subplot(gs[1:3, 0])
    plt.plot(runner.mon.ts, runner.mon['syn.g'], label='g', color='g')
    plt.legend(loc='upper right')
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax = fig.add_subplot(gs[3:6, 0])
    plt.plot(runner.mon.ts, runner.mon['post.input'], label='PSC', color='b')
    plt.legend(loc='upper right')
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax = fig.add_subplot(gs[6:9, 0])
    plt.plot(runner.mon.ts, runner.mon['post.V'], label='post.V',color='k')
    plt.legend(loc='upper right')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax = fig.add_subplot(gs[9:11, 0])
    plt.plot(runner.mon.ts, runner.mon['syn.s'], label='s', color='m')
    plt.legend(loc='upper right')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.xlabel('Time [ms]')
    plt.show()
run_syn(AMPA,title='AMPA Receptor',sp_times=[25,50,75,100,160],g_max=0.2)

50/100ms给的刺激并没有导致post.V大幅度的升高，是因为post.input(PSC)在经过了25/75ms的动作电位发放后降低为负值，对其积分后的值的减小，即为V的增量减小，导致V无法达到翻转膜电位。（此时处于一个“不应期”？）

3.3 GABAa模型

GABAa受体是抑制型受体，突触的翻转膜电位低于静息膜电位，设置翻转膜电位为-80mV，其余均采用与AMPA模型相同的步骤，编写代码如下：

class GABAa(AMPA):
    def __init__(self,pre,post,conn,g_max=0.2,E=-80.,alpha=0.53,beta=0.18,delay_step=0,
                 T_duration=1.,T=1.,method='exp_auto'):
        super(GABAa,self).__init__(pre,post,conn, g_max=g_max, E=E, alpha=alpha,
                                   beta=beta, delay_step=delay_step, T_duration=T_duration, T=T, method=method)
run_syn(GABAa,title='GABAa$-\mathrm{A}$ Synapse Model',sp_times=[25,50,75,100,160])

由图可见，GABAa传递的是抑制性信号，每次突触电导g增加意味着突触后膜电位post.V的降低。

3.4 NMDA模型

NMDA受体的特殊之处在于离子通道的开放与关闭不仅受神经递质的调节，还受膜外镁离子的影响。若突触后膜的膜电位较低，NMDA受体的离子通道会被镁离子堵住，无法让神经递质通过。只有突触后神经元的膜电位去极化后，膜内外电场发生转变，镁离子才会打开离子通道。

通过电导模式对NMDA模型进行建模，在原始表达式中增加一项，以刻画镁离子对NMDA受体的影响：

$I_{syn}=g(t)b(V,[Mg^{2+}])(V_{post}-E_{syn})$ ，其中 $b=(1+e^{-0.062V}\cdot \frac{[Mg^{2+}]}{3.57})^{-1}$ 。用x代表受体离子通道通量大小，采用二级动力学方程对s建模： $\left\{\begin{matrix} \frac{ds}{dt}=\alpha_1x(1-s)-\beta_1s \\ \frac{dx}{dt}=\alpha_2[T](1-s)-\beta_2s \end{matrix}\right.$ ，[T]的计算方法同上。代码实现如下：

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as plt

class NMDA(bp.dyn.TwoEndConn):
    def __init__(self,pre,post,conn,g_max=0.02,E=0.,c_Mg=1.2,alpha1=2.,beta1=0.01,alpha2=0.2,
                 beta2=0.5,delay_step=2,T=1.,T_duration=1.,method='exp_auto'):
        super(NMDA,self).__init__(pre=pre,post=post,conn=conn)
        self.g_max = g_max
        self.E = E
        self.c_Mg = c_Mg
        self.alpha1 = alpha1
        self.beta1 = beta1
        self.alpha2 = alpha2
        self.beta2 = beta2
        self.T = T
        self.T_duration = T_duration
        self.delay_step = delay_step
        self.pre2post = self.conn.require('pre2post')
        self.x = bm.Variable(bm.zeros(self.post.num))
        self.s = bm.Variable(bm.zeros(self.post.num))
        self.g = bm.Variable(bm.zeros(self.post.num))
        self.b = bm.Variable(bm.zeros(self.post.num))
        self.spike_arrival_time = bm.Variable(bm.ones(self.post.num)* -1e7 )
        self.delay = bm.LengthDelay(self.pre.spike,delay_step)
        self.integral = bp.odeint(method=method,f=bp.JointEq(self.ds,self.dx))
    def ds(self,s,t,x):
        return self.alpha1*x*(1-s)-self.beta1*s
    def dx(self,x,t,T):
        return self.alpha2*T*(1-x)-self.beta2*x
    def update(self,tdi):
        t,dt=tdi.t,tdi.dt
        delayed_pre_spike = self.delay(self.delay_step)
        self.delay.update(self.pre.spike)
        self.spike_arrival_time.value = bm.where(delayed_pre_spike,t,self.spike_arrival_time)#神经递质到达突触后膜的时间
        T = ((t-self.spike_arrival_time)<self.T_duration)*self.T
        #更新s,x,g,b
        self.s.value,self.x.value = self.integral(self.s,self.x,tdi.t,T,tdi.dt)
        self.g.value = self.g_max * self.s
        self.b.value = 1/(1+bm.exp(-0.062*self.post.V)*self.c_Mg/3.57)
        #电导模式下计算突触后电流
        self.post.input += self.g * self.b * (self.E - self.post.V)
def run_syn(syn_model, title, run_duration=200., sp_times=(10,20,30), **kwargs):
    neu1 = bp.neurons.SpikeTimeGroup(1, times=sp_times, indices=[0]*len(sp_times))
    neu2 = bp.neurons.HH(1, V_initializer=bp.init.Constant(-70.68))
    inputs,dur=bp.inputs.section_input(values=[0.,8.,0.],
                                      durations=[130.,1.,69.],return_length=True)
    syn = syn_model(neu1, neu2, conn=bp.connect.All2All(), **kwargs)
    net = bp.dyn.Network(pre=neu1, syn=syn, post=neu2)
    # 进行模拟
    runner = bp.dyn.DSRunner(net, monitors=['pre.spike', 'post.V', 'syn.g', 'post.input','syn.b']
                             ,inputs=[('post.input',inputs,'iter')])
    runner.run(run_duration)

    # 可视化
    fig, gs = bp.visualize.get_figure(9, 1, 0.5, 6.)
    ax = fig.add_subplot(gs[0:1, 0])
    plt.plot(runner.mon.ts, runner.mon['pre.spike'], label='pre.spike',color='r')
    plt.legend(loc='upper right')
    plt.title(title)
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['right'].set_visible(False)

    ax = fig.add_subplot(gs[1:3, 0])
    plt.plot(runner.mon.ts, runner.mon['syn.g'], label='g', color='g')
    plt.legend(loc='upper right')
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax = fig.add_subplot(gs[3:5, 0])
    plt.plot(runner.mon.ts, runner.mon['post.input'], label='PSC', color='b')
    plt.legend(loc='upper right')
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax = fig.add_subplot(gs[5:7, 0])
    plt.plot(runner.mon.ts, runner.mon['post.V'], label='post.V',color='k')
    plt.legend(loc='upper right')
    ax = fig.add_subplot(gs[7:9, 0])
    plt.plot(runner.mon.ts, runner.mon['syn.b'], label='b')
    plt.legend(loc='upper right')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.xlabel('Time [ms]')
    plt.show()
run_syn(NMDA,title='NMDA synapse',sp_times=[25,50,75,100,160])

为了观察b对模型的影响，在130ms时对突触后神经元施加一个脉冲电流，使其产生动作电位。结果显示，一开始时虽然突触前神经元接收到刺激使得g增加，但是b的值很小，因为受体几乎全被镁离子堵住，因此无法产生突触后电流。在130ms对突触后神经元施加脉冲电流时，突触后神经元的膜电位去极化，镁离子不再堵塞，b增大，从而post.V显著增大。

3.5 GABAb模型

GABAb受体与GABA结合后激活G蛋白，G蛋白结合到钾离子通道上使其开放，钾离子外流使得细胞膜超极化。用[G]表示胞内被激活的G蛋白浓度，其受两个部分调控：一是在r的刺激下以k1速率被激活，而是自身以k2的速率失活。用s表示钾离子通道开放概率，Kd表示G蛋白与钾离子通道的解离常数。G蛋白与钾离子通道结合的动力学方程由其配位方程给出，GABAb的建模公式为：

$\left\{\begin{matrix} \frac{dr}{dt}=\alpha[T](1-r)-\beta r \\ \frac{d[G]}{dt}=k_1r-k_2[G] \\ s = \frac{[G]^4}{[G]^4+K_d} \end{matrix}\right.$

搭建GABAb模型如下：

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as plt

class GABAb(bp.dyn.TwoEndConn):
    def __init__(self,pre,post,conn,g_max=1.,E=-95.,alpha=0.09,beta=0.0012,T_0=0.5,
                 T_dur=0.5,k1=0.18,k2=0.034,k_d=0.1,delay_step=2,method='exp_auto',**kwargs):
        super(GABAb,self).__init__(pre=pre,post=post,conn=conn,**kwargs)
        self.g_max = g_max
        self.E = E
        self.alpha = alpha
        self.beta = beta
        self.T_0 = T_0
        self.T_dur = T_dur
        self.k1 = k1
        self.k2 = k2
        self.k_d = k_d
        self.delay_step = delay_step
        self.pre2post = self.conn.require('pre2post')
        self.r = bm.Variable(bm.zeros(self.post.num))
        self.G = bm.Variable(bm.zeros(self.post.num))
        self.s = bm.Variable(bm.zeros(self.post.num))
        self.g = bm.Variable(bm.zeros(self.post.num))
        self.spike_time_arrival = bm.Variable(bm.ones(self.post.num)*-1e7)
        self.delay = bm.LengthDelay(self.pre.spike,delay_step)
        self.integral = bp.odeint(method=method,f=bp.JointEq(self.dr,self.dG))
    def dr(self,r,t,T):
        return self.alpha*T*(1-r) - self.beta*r
    def dG(self,G,t,r):
        return self.k1*r - self.k2*G
    def update(self,tdi):
        t,dt = tdi.t,tdi.dt
        delayed_pre_spike = self.delay(self.delay_step)
        self.delay.update(self.pre.spike)
        self.spike_time_arrival = bm.where(delayed_pre_spike,t,self.spike_time_arrival)
        T = ((t-self.spike_time_arrival)<self.T_dur)*self.T_0
        self.r.value,self.G.value = self.integral(self.r,self.G,t,T,dt)
        self.s.value = bm.power(self.G,4)/(bm.power(self.G,4)+self.k_d)
        self.g.value = self.g_max*self.s
        self.post.input += self.g*(self.E-self.post.V)
def run_syn(syn_model, title, run_duration=1000., sp_times=(130), **kwargs):
    neu1 = bp.neurons.SpikeTimeGroup(1, times=sp_times, indices=[0]*len(sp_times))
    neu2 = bp.neurons.HH(1, V_initializer=bp.init.Constant(-70.68))
    syn = syn_model(neu1, neu2, conn=bp.connect.All2All(), **kwargs)
    net = bp.dyn.Network(pre=neu1, syn=syn, post=neu2)
    # 进行模拟
    runner = bp.dyn.DSRunner(net, monitors=['pre.spike', 'syn.r', 'syn.g','syn.G'])
    runner.run(run_duration)

    # 可视化
    fig, gs = bp.visualize.get_figure(10, 1, 0.5, 6.)
    ax = fig.add_subplot(gs[0:2, 0])
    plt.plot(runner.mon.ts, runner.mon['pre.spike'], label='pre.spike',color='r')
    plt.legend(loc='upper right')
    plt.title(title)
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['bottom'].set_visible(False)
    ax.spines['right'].set_visible(False)

    ax = fig.add_subplot(gs[2:6, 0])
    plt.plot(runner.mon.ts, runner.mon['syn.r'], label='r', color='g')
    plt.plot(runner.mon.ts, runner.mon['syn.G'] / 4, label='G/4', color='b')
    plt.legend(loc='upper right')
    plt.xticks([])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax = fig.add_subplot(gs[6:10, 0])
    plt.plot(runner.mon.ts, runner.mon['syn.g'], label='g',color='m')
    plt.legend(loc='upper right')
    plt.xlabel('Time [ms]')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    plt.show()
run_syn(GABAb, title='GABAb synapse', sp_times=[100])

由于GABAb模型中各个变量的时间常数较大，g的变化较缓慢，因此进行模拟时只让突触前神经元产生一个脉冲。

由图可见，GABAb模型的各个参数变化十分缓慢，这种长时程的作用使得突触后电位长时间处于超极化状态，突触后神经元在这段时间内呈现出被抑制的状态。

3.6 生理学模型与现象学模型的比较

通过动力学构建模型的优势在于，变量s存在上限，其值不可能超过1，更加符合生物学意义。而现象学模型的电导接受高频刺激后会一直增加。生理学模型限定了g的上界，保证了模型的生物合理性。

4 电突触模型

4.1 电突触模型概述

在电突触中，两个神经元之间的连接结构称为连接间隙，两个神经元的细胞膜靠的非常近，允许离子自由通过且没有选择性。电流的流动方向是双向的，其产生是被动的。

电突触的数学模型没有微分方程，由等效电路可得出其表达式为： $I_1=g_c(V_1-V_2)$

4.2 电突触模型的代码实现

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as plt

class GapJunction(bp.dyn.TwoEndConn):
    def __init__(self,pre,post,conn,g=0.2,**kwargs):
        super(GapJunction,self).__init__(pre=pre,post=post,conn=conn,**kwargs)
        self.g =g
        self.pre_ids,self.post_ids=self.conn.require('pre_ids','post_ids')
        self.current = bm.Variable(bm.zeros(self.post.num))
    def update(self,tdi):
        t,dt=tdi.t,tdi.dt
        inputs = self.g*(self.pre.V[self.pre_ids] - self.post.V[self.post_ids])
        self.current.value = bm.syn2post(inputs,self.post_ids,self.post.num)
        self.post.input += self.current
def run_syn_GJ(syn_model,title,run_duration=100.,Iext=7.5,**kwargs):
    neu = bp.dyn.HH(2)
    syn = syn_model(neu,neu,conn = bp.connect.All2All(include_self=False),**kwargs)
    net = bp.dyn.Network(syn=syn,neu=neu)
    runner = bp.dyn.DSRunner(net,inputs=[('neu.input',bm.array([Iext,0.]))],
                             monitors=['neu.V','syn.current'],jit=True)
    runner.run(run_duration)

    fig,gs = plt.subplots(2,1,figsize=(6,4.5))
    plt.sca(gs[0])
    plt.plot(runner.mon.ts,runner.mon['neu.V'][:,0],label='neu0-V')
    plt.plot(runner.mon.ts,runner.mon['neu.V'][:,1],label='neu1-V')
    plt.legend(loc='upper right')
    plt.title(title)

    plt.sca(gs[1])
    plt.plot(runner.mon.ts, runner.mon['syn.current'][:, 0], label='neu1-current',color=u'#48d688')
    plt.plot(runner.mon.ts, runner.mon['syn.current'][:, 1], label='neu0-current',color=u'#d64888')
    plt.legend(loc='upper right')
    plt.tight_layout()
    plt.show()
run_syn_GJ(GapJunction,Iext=7.5,title='Gap Junction Model(Iext=7.5)')
run_syn_GJ(GapJunction,Iext=5.,title='Gap Junction Model(Iext=5.0)')

模型使用self.current 变量来记录通过间隙连接的电流（其不包括外部输入的电流）。将外部电流分别设置为7.5和5.0，观察膜电位变化以及间隙连接中电流的变化：