Variables stepping out of bounds, membrane potential going to infinity

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):
    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}
    # 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)

    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)

Hi @fleurzeldenrust :wave:
I did not yet have time to look into this in detail, but here are a few quick remarks. First of all: did you consider using a different integration method (e.g. rk2 or rk4)?

The TimedArray can have a different dt from your simulation. E.g. you could have a TimedArray that changes its values every 0.1ms, but your simulation dt is 0.05*ms. The only constraint is that the dt values are multiples of each other.

A brute force way is to use run_regularly after the integration step to force values into your desired range (and then the reset should take care of the rest):

G.run_regularly('v = clip(v, -200*mV, 0*mV)', when='after_groups')

O wow, I wasn’t aware of that function of the timed array, that is making my life so much easier! The array we use is limited by the experimental setup, and I wanted to use the same one as in the experiments…

Thank you!

1 Like