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

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

书接上文：

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

3 化学突触的生理学模型

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

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

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

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

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 pltclass 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_maxself.E = Eself.alpha = alphaself.beta = betaself.delay_step = delay_stepself.T_duration = T_durationself.T = Tself.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*sdef 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.Tself.s.value = self.integral(self.s,tdi.t,TT,tdi.dt)self.g.value = self.g_max * self.sself.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受体的影响：

，其中。用x代表受体离子通道通量大小，采用二级动力学方程对s建模：，[T]的计算方法同上。代码实现如下：

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as pltclass 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_maxself.E = Eself.c_Mg = c_Mgself.alpha1 = alpha1self.beta1 = beta1self.alpha2 = alpha2self.beta2 = beta2self.T = Tself.T_duration = T_durationself.delay_step = delay_stepself.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*sdef dx(self,x,t,T):return self.alpha2*T*(1-x)-self.beta2*xdef update(self,tdi):t,dt=tdi.t,tdi.dtdelayed_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,bself.s.value,self.x.value = self.integral(self.s,self.x,tdi.t,T,tdi.dt)self.g.value = self.g_max * self.sself.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的建模公式为：

搭建GABAb模型如下：

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as pltclass 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_maxself.E = Eself.alpha = alphaself.beta = betaself.T_0 = T_0self.T_dur = T_durself.k1 = k1self.k2 = k2self.k_d = k_dself.delay_step = delay_stepself.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*rdef dG(self,G,t,r):return self.k1*r - self.k2*Gdef update(self,tdi):t,dt = tdi.t,tdi.dtdelayed_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_0self.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.sself.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 电突触模型概述

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

电突触的数学模型没有微分方程，由等效电路可得出其表达式为：

4.2 电突触模型的代码实现

import brainpy as bp
import brainpy.math as bm
import numpy as np
import matplotlib.pyplot as pltclass 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 =gself.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.dtinputs = 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，观察膜电位变化以及间隙连接中电流的变化：

当第 0 个神经元膜电位升高时，第 1 个神经元将通过间隙连接接收到一个正向的电流，于是膜电位也随之升高；第 0 个神经元发放脉冲也会引发第 1 个神经元发放脉冲。此外，因为只有一个间隙连接存在，因此两个神经元接收到的电流之和总是为 0。

当外部输入电流为5.0时，不能让神经元产生动作电位，因为电突触像一个漏电装置，对第0个神经元来说，向外流出的电流i小了外部输入电流的部分效果，因此需要更大的外部电流才能激发动作电位。

5 总结

静态突触的建模在生理结构上分为化学突触和电突触模型。化学突触从建模方法上又可以分为现象学模型和生理学模型。本章节介绍了现象学模型中的电压跳变模型、指数衰减模型、双指数衰减模型等，生理学模型中的AMPA\NMDA\GABAa\GABAb模型以及电突触模型。