Neuronal Dynamics:Finite cable model 实验结果展示
一、演示效果和代码
通过设置输入的类型,在代码运行完之后会在本地产生一个显示输出内容的gif文件,和一个包含u-x和u-t图的png图片
当设置type=1的时候,相当于设定在t=0ms时在x=1输入一个冲激脉冲,会输出如下显示u-x随时间t变化的动态图:
当设置type=2的时候,相当于设定在t=0~1ms的时段内,在x=1持续输入的一个脉冲,会输出如下显示u-x随时间t变化的动态图:
当设置type=3的时候,相当于设定在t=0ms时在x=1输入一个冲激脉冲,同时在t=0~0.1ms的时段内,在x=4持续输入的一个脉冲,会输出如下显示u-x随时间t变化的动态图:
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.pyplot import *
import matplotlib.animation
rc('font',size=16)
rc('font',family='SimHei')
plt.rcParams['axes.unicode_minus'] = False
class Cable():
def __init__(self,L,all_time):
self.L=L
self.all_time=all_time
self.dt=0.01
self.dx=0.1
self.I_list = [[0 for j in range(int(L/self.dx))] for i in range(int(all_time/self.dt))]
self.x_list = [i*self.dx for i in range(int(L/self.dx))]
self.t_list = [i*self.dt for i in range(int(all_time/self.dt))]
self.u_list = [[0 for j in range(int(L/self.dx))] for i in range(int(all_time/self.dt))]
self.type=0
def change(self,type):
self.type=type
if type==0:
self.I_list = [[0 for j in range(int(self.L/self.dx))] for i in range(int(self.all_time/self.dt))]
elif type==1:
self.I_list = [[(300 if (i*self.dt==0 and j==int(1/self.dx)) else 0) for j in range(int(self.L/self.dx))] for i in range(int(self.all_time/self.dt))]
elif type==2:
self.I_list = [[(10 if (i*self.dt<1 and j==int(1/self.dx)) else 0) for j in range(int(self.L/self.dx))] for i in range(int(self.all_time/self.dt))]
elif type==3:
self.I_list = [[(30 if (i!=0 and i*self.dt<0.1 and j==int(4/self.dx)) else (300 if (i*self.dt==0 and j==int(1/self.dx)) else 0)) for j in range(int(self.L/self.dx))] for i in range(int(self.all_time/self.dt))]
def get_G(self,t,x):
if t>0:
return np.exp(-t-x*x/4/t)/np.sqrt(4*np.pi*t)
return 0
def get_G_limit(self,t0,x0,t,x):
sum = 0
for n in range(-10,10):
sum+=self.get_G(t-t0,x-2*n*self.L-x0)+self.get_G(t-t0,x-2*n*self.L+x0)
return sum
def caculate_u(self):
t0max = int(self.all_time/self.dt)
x0max = int(self.L/self.dx)
for x in range(x0max):
for t in range(t0max):
#通用计算过程
for t0 in range(t):
for x0 in range(x0max):
if self.I_list[t0][x0]!=0:
self.u_list[t][x]+=self.get_G_limit(t0*self.dt,x0*self.dx,t*self.dt,x*self.dx)*self.I_list[t0][x0]
#type1用下列代码计算更快
# for t0 in range(0,1):
# self.u_list[t][x] += self.get_G_limit(t0*self.dt, 1, t*self.dt, x*self.dx) * 300
self.u_list[t][x] *=self.dx*self.dt
def plot_line_pic(x,y,label='u-t变化图',ylim1=0,ylim2=0):
plt.plot(x, y, '--r',linewidth=1.6, markeredgecolor='k', markersize=8)
title(label,fontsize=20)
if ylim1!=0 or ylim2!=0:
ylim(ylim1,ylim2)
if __name__ == '__main__':
model = Cable(L=5,all_time=1)
# model.change(type=1)
# model.change(type=2)
model.change(type=3)
model.caculate_u()
subplot(121)
if model.type==1:
for i in range(3,40,6):
plot_line_pic(model.x_list, model.u_list[i], label='u-x图;type='+str(model.type))
else:
for i in range(0,int(model.all_time/model.dt),int(model.all_time/model.dt/7)):
plot_line_pic(model.x_list, model.u_list[i], label='u-x图;type='+str(model.type))
subplot(122)
for i in range(0,int(model.L/model.dx),(int((model.L/model.dx)/7))):
plot_line_pic(model.t_list, [model.u_list[j][i] for j in range(len(model.t_list))], label='u-t图;type='+str(model.type))
plt.savefig("Finite_cable_model_"+str(model.type)+".png")
show()
global fig
fig, ax = plt.subplots()
def update(i):
plt.cla()
plot_line_pic(model.x_list, model.u_list[int(i)], label='u-x图;type='+str(model.type)+';t='+str(round(i*model.dt,2)),ylim2=1)
ani = matplotlib.animation.FuncAnimation(fig, update, frames=50, interval=100, repeat=True)
ani.save("Finite_cable_model_"+str(model.type)+".gif", writer='imgemagick')
二、模型简述
该模型是对有限长度的链型神经结构的模拟:
由cable equation可以推出格林函数,即单位冲击函数对无穷长的链产生的影响为:
G
∞
(
t
,
x
)
=
Θ
(
t
)
4
π
t
e
x
p
(
−
t
−
x
2
4
t
)
G_\infty(t,x)=\frac{\Theta(t)}{\sqrt{4\pi t}}exp(-t-\frac{x^2}{4t})
G∞(t,x)=4πtΘ(t)exp(−t−4tx2)
考虑到限制函数:
∂
∂
x
u
(
t
,
x
=
0
)
=
∂
∂
x
u
(
t
,
x
=
L
)
=
0
\frac{\partial}{\partial x}u(t,x=0)=\frac{\partial}{\partial x}u(t,x=L)=0
∂x∂u(t,x=0)=∂x∂u(t,x=L)=0
由mirror charges技巧,单位冲击函数对有限长为L的链产生的影响为:
G
0
,
L
(
t
0
,
x
0
;
t
,
x
)
=
∑
n
=
−
∞
+
∞
G
∞
(
t
−
t
0
,
x
−
2
n
L
−
x
0
)
+
G
∞
(
t
−
t
0
,
x
−
2
n
L
+
x
0
)
G_{0,L}(t_0,x_0;t,x)=\sum_{n=-\infty}^{+\infty}G_\infty(t-t_0,x-2nL-x_0)+G_\infty(t-t_0,x-2nL+x_0)
G0,L(t0,x0;t,x)=n=−∞∑+∞G∞(t−t0,x−2nL−x0)+G∞(t−t0,x−2nL+x0)
结合输入单位时间和单位长度的外部电流输入iext(t,x)有:
u
(
t
,
x
)
=
∫
−
∞
t
d
t
0
∫
0
L
d
x
0
G
0
,
L
(
t
0
,
x
0
;
t
,
x
)
i
e
x
t
(
t
0
,
x
0
)
u(t,x)=\int_{-\infty}^{t}dt_0\int_0^Ldx_0G_{0,L}(t_0,x_0;t,x)i_{ext}(t_0,x_0)
u(t,x)=∫−∞tdt0∫0Ldx0G0,L(t0,x0;t,x)iext(t0,x0)
三、输出解释
3.1输入为t=0ms时位于x=1处的冲激函数
左图为0.03ms开始每隔0.06ms记录一次的u-x图结果
右图为在X=0与X=L之间均匀分布的七个位置的u随时间t的变化
在冲激函数的作用下,首先x=1处的幅度大幅上升,然后渐渐地下降,离x=1位置越远的地方摆动产生的幅度越小
3.2输入为t=0~1ms持续在x=1处输入脉冲
左图为在t=0与t=1之间均匀分布的七个时间点的u随时间x的变化
右图为在X=0与X=L之间均匀分布的七个位置的u随时间t的变化
在阶跃函数的作用下,各处的u值均持续上升,x=1处的幅度上升最大,离x=1位置越远的地方摆动产生的幅度越小
3.3输入为t=0ms位于x=1一个冲激脉冲,加上t=0~0.1ms的时段内,位于x=4持续输入的一个脉冲
可惜并没有出现和compartmental model一样的如下图的效果,推测是由于compartmental model考虑了分支的存在以及每一段的生物性质的差异,而我们这个模型中默认只有一条链