本文使用AdEx神经元搭建一个完整的神经网络来进行生物神经脉冲现象的仿真。主要的目的是为了验证数学原理,因此只调用的numpy函数包。对应的代码例程如下:
1.导入所需的Python函数库
import numpy as np
import matplotlib.pyplot as plt
import re
import os
2.定义均值函数以及一些常用函数
def bin_data(data):
try:
return np.mean(data[len(data)-bin_num:len(data)])
except:
return data[len(data)-1]
def extract_number_from_string(s):
match = re.search(r'_(\d+)$', s)
if match:
return int(match.group(1))
else:
return None
3.AdEx神经元定义
class AdExNeuron:
def __init__(self, name, V_Neuron, w_adaptive, G_Synapsis_Excitatory, G_Synapsis_Inhibitory,
E_Excitatory, E_Inhibitory, E_local, G_local, V_disturb, V_Excitatory_Threshold,C_Membrane,
a_w_adaptive, tau_w_adaptive,
tau_Synapsis,
V_Reset_Threshold, V_Reset, b_w_adaptive,
I_Synapsis, T_refractory, T_rest,
Connecting_Neuron, Q_Synapsis, Probability_Connecting):
# variable parameters
self.name = name
self.V_Neuron = V_Neuron
self.w_adaptive = w_adaptive
self.G_Synapsis_Excitatory = G_Synapsis_Excitatory
self.G_Synapsis_Inhibitory = G_Synapsis_Inhibitory
# fixed parameters
self.E_Excitatory = E_Excitatory
self.E_Inhibitory = E_Inhibitory
self.E_local = E_local
self.G_local = G_local
self.V_disturb = V_disturb
self.V_Excitatory_Threshold = V_Excitatory_Threshold
self.C_Membrane = C_Membrane
self.T_refractory = T_refractory
# adaptive parameters
self.a_w_adaptive = a_w_adaptive
self.tau_w_adaptive = tau_w_adaptive
self.tau_Synapsis = tau_Synapsis
# reset parameters
self.V_Reset_Threshold = V_Reset_Threshold
self.V_Reset = V_Reset
self.b_w_adaptive = b_w_adaptive
self.I_Synapsis = I_Synapsis
self.T_rest = T_rest
# connecting neurons
self.Connecting_Neuron = Connecting_Neuron
self.Q_Synapsis = Q_Synapsis
self.Probability_Connecting = Probability_Connecting
def refresh_w_adaptive(self):
# if self.T_rest<=0:
self.w_adaptive = self.w_adaptive+dt*(self.a_w_adaptive*(self.V_Neuron-self.E_local)-self.w_adaptive)/self.tau_w_adaptive
def refresh_G_Synapsis_Excitatory(self):
# if self.T_rest<=0:
self.G_Synapsis_Excitatory = self.G_Synapsis_Excitatory-dt*self.G_Synapsis_Excitatory/self.tau_Synapsis
def refresh_G_Synapsis_Inhibitory(self):
# if self.T_rest<=0:
self.G_Synapsis_Inhibitory = self.G_Synapsis_Inhibitory-dt*self.G_Synapsis_Inhibitory/self.tau_Synapsis
def refresh_membrane_potential(self):
if self.T_rest<=0:
self.V_Neuron =self.V_Neuron+dt*(self.G_Synapsis_Excitatory*(self.E_Excitatory-self.V_Neuron)+
self.G_Synapsis_Inhibitory*(self.E_Inhibitory-self.V_Neuron)+
self.G_local*(self.E_local-self.V_Neuron)+
self.G_local*self.V_disturb*np.exp((self.V_Neuron-self.V_Excitatory_Threshold)/self.V_disturb)-
self.w_adaptive+self.I_Synapsis
)/self.C_Membrane
else:
self.T_rest=self.T_rest-dt
self.refresh_w_adaptive()
self.refresh_G_Synapsis_Excitatory()
self.refresh_G_Synapsis_Inhibitory()
def fire(self, num1, num2):
global fire_matrix
# refresh self parameter
self.V_Neuron = self.V_Reset
self.w_adaptive = self.w_adaptive+self.b_w_adaptive
self.T_rest=self.T_refractory
# refresh the G_Synapsis
# print(self.name)
if self.name[1]=='1':
num1=num1+1
fire_matrix1[extract_number_from_string(self.name)-1,test_input_index]=2
for neuron1 in self.Connecting_Neuron:
neuron1.G_Synapsis_Inhibitory=neuron1.G_Synapsis_Inhibitory+self.Q_Synapsis
if self.name[1]=='2':
num2=num2+1
fire_matrix2[extract_number_from_string(self.name)-1,test_input_index]=2
for neuron2 in self.Connecting_Neuron:
neuron2.G_Synapsis_Excitatory=neuron2.G_Synapsis_Excitatory+self.Q_Synapsis
return num1, num2
def judge_fire(self, num1, num2):
if self.V_Neuron>self.V_Reset_Threshold:
num1, num2=self.fire(num1, num2)
else:
pass
return num1, num2
def Add_Synapsis(self, Synapsis):
self.Connecting_Neuron.append(Synapsis)
4.连接神经元
for neuron_front in G_Group:
for neuron_back in G_Group:
if neuron_front !=neuron_back:
if np.random.rand()<neuron_front.Probability_Connecting:
neuron_front.Connecting_Neuron.append(neuron_back)
5.将神经元接入输入
for neuron_front in P2_Group:
for neuron_back in G_Group:
if neuron_front !=neuron_back:
if np.random.rand()<neuron_front.Probability_Connecting:
neuron_front.Connecting_Neuron.append(neuron_back)
6.定义输入函数脉冲
def heaviside(x):
return 0.5 * (1 + np.sign(x))
def input_rate(t, t1_exc, tau1_exc, tau2_exc, ampl_exc, plateau):
# t1_exc=10. # time of the maximum of external stimulation
# tau1_exc=20. # first time constant of perturbation = rising time
# tau2_exc=50. # decaying time
# ampl_exc=20. # amplitude of excitation
inp = ampl_exc * (np.exp(-(t - t1_exc) ** 2 / (2. * tau1_exc ** 2)) * heaviside(-(t - t1_exc)) +
heaviside(-(t - (t1_exc+plateau))) * heaviside(t - t1_exc) +
np.exp(-(t - (t1_exc+plateau)) ** 2 / (2. * tau2_exc ** 2)) * heaviside(t - (t1_exc+plateau)))
return inp
7.开始仿真
for tick_time in np.arange(0, TotTime, dt):
fire_probability=dt*test_input[test_input_index]
if test_input_index%5000==0:
print(test_input_index)
print("fire_probability:"+str(fire_probability))
for neuron in P2_Group:
if np.random.rand()<fire_probability:
neuron.fire(0,0)
fire1_num=0
fire2_num=0
fire1_frequent=0
fire2_frequent=0
for neuron in G_Group:
neuron.refresh_membrane_potential()
fire1_num, fire2_num=neuron.judge_fire(fire1_num, fire2_num)
fire1_frequent=fire1_num/dt/N1
fire2_frequent=fire2_num/dt/N2
fire1_result.append(fire1_frequent)
fire2_result.append(fire2_frequent)
neuron1_potential_bin.append(G1_1.V_Neuron) # type: ignore
neuron2_potential_bin.append(G2_1.V_Neuron) # type: ignore
fire1_result_bin.append(bin_data(fire1_result))
fire2_result_bin.append(bin_data(fire2_result))
test_input_index=test_input_index+1
8.仿真结果
下面顺序依次是:
兴奋性神经元与抑制性神经元的膜电位
输入点火频率与两类神经元的点火频率
总体点火频率与对应神经元点火节点
抽取局部神经元序列点火时间节点