Description of problem
I am running a standard AdEx model, with a fluctuating (timed array) input. Each run, I increase the amplitude of this frozen noise input. If the amplitude becomes too large, at a certain point the membrane potential goes to infinity (see image below). Now this is a well-known problem for neuron models and not unique to Brian. However, I would like to do a hard reset when this happens, i.e. force the variables back to within their range. Unfortunately I cannot decrease the time step, due to the limitations of the timed array. Is there a way to force the membrane potential, or other variables, back into their normal range?
Minimal code to reproduce problem
Sorry, I would have to upload the timed array, I’m not sure I can. Please let me know how.
def define_Sim_Parameters(ex_inh, inputfolder):
if ex_inh == 'excitatory':
inputvariable ="input_theory_e"
inputfile = 'input_' + ex_inh + '_tau_250_sampling_40_kHz.mat'
elif ex_inh == 'inhibitory':
inputfile = 'input_' + ex_inh + '_tau_50_sampling_40_kHz.mat'
inputvariable ="input_theory_i"
data = inputfolder + inputfile
hidden_state = np.squeeze(sio.loadmat(data)["hidden_state_" + ex_inh[0]])
# max tmax 360000*ms
Sim_Parameters = {
# simulator settings
"dt" : 0.025*ms,
"tmax" : 5000*ms,
"integrator" : euler,
# input location
"inputfolder" : inputfolder,
"inputfile" : inputfile,
"inputvariable": inputvariable}
return Sim_Parameters, hidden_state
Parameters = {
"gL" : 10*nS,
"Cm" : 50*pF,
"EL" : -70*mV,
"dT" : 1*mV,
# spike frequency adaptation
"tauw" : 144*ms, # time constant adaptation
"a" : 0*nS, # sets the amount of adaptation
"b" : 0*nA, # jumpsize after spike
# threshold adaptation
"tauth" : 5*ms, # time constant threshold adaptation
"Ka" : 0*mV,
"Vt" : -63*mV, # unadapted threshold
"Vi" : -67*mV,
"p" : 0, # sets the Vm dependency
"ki": 5*mV,
"refractory" : 0.8*ms,
# input scaling
"I_baseline" : 0,
"I_scale" : 1500,
# initial values
'v_init' : -70*mV}
Parameters["Vcut"]=Parameters['Vt'] + 5 * Parameters['dT']
Parameters["Vr"]= Parameters['EL']
Parameters["w_init"]= Parameters['a']*(Parameters['v_init'] - Parameters['EL'])
Parameters["th_init"] = Parameters['Vt']
### AdEx model with additional threshold adaptation
def run_simulation(Parameters, Sim_Parameters):
start_scope()
Cm = Parameters['Cm']
gL = Parameters['gL']
EL = Parameters['EL']
dT = Parameters['dT']
Vcut = Parameters['Vcut']
tauw = Parameters['tauw']
a = Parameters['a']
b = Parameters['b']
p = Parameters['p']
ki = Parameters["ki"]
Vr = Parameters['Vr']
tauth = Parameters["tauth"]
Ka = Parameters["Ka"]
Vt = Parameters["Vt"]
Vi = Parameters["Vi"]
inputfolder = Sim_Parameters['inputfolder']
inputfile = Sim_Parameters['inputfile']
# time settings
defaultclock.dt = Sim_Parameters["dt"]
# initial values: sets the values at the start of the simulation
v_init = Parameters['v_init']
w_init = Parameters['w_init']
th_init = Parameters['th_init']
# loading the input and scaling it to make it suitable for the model
I_theory = np.squeeze(sio.loadmat(inputfolder+inputfile)[Sim_Parameters['inputvariable']])
I_baseline = Parameters["I_baseline"]
I_scale = Parameters["I_scale"]
stim = TimedArray(((I_scale*I_theory+I_baseline)*pamp), dt = defaultclock.dt)
# the model's equations
eqs = '''
dVm/dt = (gL*(EL - Vm) + gL*dT*exp((Vm - th)/dT) + I - w)/Cm : volt
dw/dt = (a*(Vm - EL) - w)/tauw : amp
dth/dt = (th_inf-th)/tauth : volt
th_inf = p*(Vm-Vi) + Vt + Ka*log(1+exp((Vm-Vi)/ki)) : volt
I = stim(t) : amp
'''
# creates a neuron and specifies its properties
G = NeuronGroup(1, model = eqs, threshold='Vm>Vcut',
reset="Vm=Vr; w+=b", refractory = Parameters["refractory"], method = Sim_Parameters["integrator"])
initial_values = {'Vm': v_init, 'w': w_init, 'th': th_init}
G.set_states(initial_values)
# set monitors for the variables
S = SpikeMonitor(G)
M = StateMonitor(G, 'Vm', record = True)
if (adaptation_model == 'sf_adaptation') or (adaptation_model == 'combined_adaptation'):
W = StateMonitor(G, 'w', record = True)
if (adaptation_model == 'threshold_adaptation') or (adaptation_model == 'combined_adaptation'):
Th = StateMonitor(G, 'th', record = True)
run(Sim_Parameters["tmax"])
if (adaptation_model == 'no_adaptation'):
return M, S
elif (adaptation_model == 'sf_adaptation'):
return M, S, W
elif (adaptation_model == 'threshold_adaptation'):
return M, S, Th
elif (adaptation_model == 'combined_adaptation'):
return M, S, W, Th
M, S = run_simulation(Parameters, Sim_Parameters)
