Store/Restore Stochastic Disruption on Start up

Support
1

Description of problem

I’m using the store/restore system in my model, but it’s not performing as I expected. I run the model for a 50-second spin-up period, and then store the network. I then restore the network and run it for a 300 ms experimental period. However, when I first restore the network, the PoissonInput objects that I’m using for noise in the network seem to spike before settling down again. Is this expected behavior, or is there something that I’m doing incorrectly that is causing the jump?

I’ll also add that I am performing these tests in a Jupyter notebook, and the store/restore are in different cells, but I don’t imagine that to be the cause?

Minimal code to reproduce problem

# network definition 

Thal.I_app = 0*pA
L4e.I_app = 0*pA
L4i.I_app = 0*pA
Ret.I_app = 0*pA
Spike_t = SpikeMonitor(Thal)
Spike_e = SpikeMonitor(L4e)
Spike_i = SpikeMonitor(L4i)
Spike_r = SpikeMonitor(Ret)

run(Init_time, report='stdout')
store()

# a few steps dedicated to establishing directories for output data

restore()

del Spike_t, Spike_e, Spike_i, Spike_r
Spike_t = SpikeMonitor(Thal)
Spike_e = SpikeMonitor(L4e)
Spike_i = SpikeMonitor(L4i)
Spike_r = SpikeMonitor(Ret)
setParam('Istim', str(i) + '*pA', params)
setParam('StimCx', str(s), params)

Thal.I_app = Istim * (1-StimCx)
L4e.I_app = Istim * StimCx
L4i.I_app = Istim * StimCx

run(run_time)

What you have aready tried

I’ve implemented a buffer window, allowing the model to run an extra 1.5 s before I reactivate the spike counters. That way, even though the spike still happens, I get a clean reading from a stable network. I’m more curious about whether or not I’m doing something wrong, or if this is the intended functionality of the store()/restore() commands.

Actual output

Below is a plot of the network activity at the end of the spin-up period and the beginning of the run-time period. The spin-up is on the left and the run-time is on the right, so the two plots can be viewed as one time axis side by side. As shown, at the beginning of the run-time period, a level of disruption occurs when the network first reactivates.

pulsingspinup_connected
pulsingspinup_connected2780×1010 244 KB

2

Hi @cwalsh27 Thanks for the report! Unfortunately, your example is missing the network definition part, so I cannot run it. I tried a simple example with PoissonInput (and PoissonGroup), but did not see anything off ­– I guess you are doing something different, though. Here’s my code for reference:

from brian2 import *

group = NeuronGroup(1, """
                       dx/dt = -x / (10*ms) : 1
                       dy/dt = -y / (10*ms) : 1  
                       """)
poisson_1 = PoissonGroup(100, rates=100*Hz)
syn_poisson = Synapses(poisson_1, group, on_pre="x += 1")
syn_poisson.connect()
poisson_2 = PoissonInput(group, "y", 100, 100*Hz, weight=1)

mon = StateMonitor(group, ["x", "y"], record=0)
run(1*second)
store()
restore()
run(500*ms)
plt.plot(mon.t/ms, mon.x[0])
plt.plot(mon.t/ms, mon.y[0])
plt.axvline(1000, color="darkgray", zorder=-1)
plt.show()

which gives

Figure_1
Figure_1640×480 29.6 KB

3

Hi @mstimberg ,

Thanks for your response, I appreciate you taking the time. Below is the network definition that I’m utilizing. Each block of code is a separate cell in a Jupyter Notebook. After the following three cells, the notebook progresses to the network runs that I linked in my original post. Thank you in advance for any insight you can provide!

# parameter dictionary definition
params = {
    'defaultclock.dt' : '0.1*ms',
    'Init_time' : '50000*ms',
    'buffer_time' : '1200*ms',
    'run_time' : '1200*ms',
    'Cm' : '1* uF/cm**2',
    'S' : '20e3 * um**2',
    'gL' : '.05 * mS/cm**2',
    'EL' : '-60 * mV',
    'Vr' : '-60 * mV',
    'Vt' : '-50 * mV',
    'DeltaT' : '2.5 * mV',
    'tauw' : '600 * ms',
    
    'N_i' : '150',
    'N_e' : '1600',
    'N_t' : '200',
    'N_r' : '100',
    
    'g_E_T' : '0.15 *mS/ cm**2',    
    'g_E_E' : '0.2 *mS/ cm**2',
    'g_E_I' : '0.7 *mS/ cm**2',
    'g_I_T' : '0.2  *mS/ cm**2',
    'g_I_E' : '0.6 *mS/ cm**2',
    'g_I_I' : '0.55 *mS/ cm**2',
    
    'g_R_E' : '0.35 *mS/ cm**2',
    'g_T_E' : '0.35 *mS/ cm**2',
    'g_R_T' : '0.085 *mS/ cm**2',
    'g_T_R' : '2.84 *mS/ cm**2',
    'g_R_R' : '2.84 *mS/ cm**2',
    
    'E_Ampa' : '0.0*mV',
    'E_Gaba' : '-80.0*mV',
    'tau_Ampa' : '2.0 *ms',
    'tau_Gaba' : '3.0 *ms',
    
    'Kappa_e_t' : '100.0',
    'Kappa_i_t' : '75.0',
    
    'Kappa_e_e' : '200.0',
    'Kappa_i_e' : '400.0',
    
    'Kappa_e_i' : '25.0',
    'Kappa_i_i' : '25.0',
    
    'Kappa_t_e' : '32.0',
    'Kappa_r_e' : '32.0',
    'Kappa_r_t' : '2.0',
    'Kappa_t_r' : '8.0',
    'Kappa_r_r' : '8.0',    
    
    
    'syn_weight' : '1',
    
    'Lambda' : '0.1 * mm',  
    'R_Cx' : '0.20 * mm',
    'R_Th' : '0.10 * mm',
    
    'Po_N_exc' : '100',
    'Po_N_inh' : '25',
    'Po_Rate': '25*Hz',         
    'Po_Amp_T': '0.2*mV',
    'Po_Amp_E': '0.2*mV',
    'Po_Amp_I': '0.07*mV',
    'Po_Amp_R': '0.2*mV',
}



def setParams(d):
    '''Sets a dictionary of parameters as global python variables.

    Args:
        d (dict): The dictionary of parameters (param_name, param_vals)
    '''
    for param, val in d.items():
        exec(param + '=' + val, globals())
        
def setParam(param, val, d):
    '''Sets a single parameter/value pair as a global python variable and adds it to data export dict.

    Args:
        param (str): The name of the parameter variable
        val (str): The value of the parameter variable
        d (dict): The dictionary of parameters (param_name, param_vals)
    '''
    exec(param + '=' + val, globals())
    d[param] = val

setParams(params)


#AdEx IF Model 
eqs_aeIF = """
dv/dt = (-gL*(v-EL) + gL*DeltaT*exp((v - Vt)/DeltaT) + I - w/S)/Cm : volt (unless refractory)
dw/dt = (a*(v - EL) - w)/tauw : ampere
I = I_inp/S - I_syn: ampere / meter**2
I_inp = I_app * exp(-r/Lambda) : ampere 
r = i*Rad/N : meter
Rad : meter
I_app : ampere
a: siemens
b: ampere
"""

#synapse current equations, same for each neuron group (from Brunel/Wang (ds_Ampa_T/dt) and Gutnisky)
eqs_I_Ampa_T = '''
I_Ampa_T = g_syn_T * (v - E_Ampa) * s_Ampa_T * 1*ms/(sqrt(Kappa_T)*tau_Ampa): ampere / meter**2
ds_Ampa_T / dt = - s_Ampa_T / tau_Ampa : 1
g_syn_T: siemens / meter**2
Kappa_T :1
''' 
eqs_I_Ampa_E = '''
I_Ampa_E = g_syn_E * (v - E_Ampa) * s_Ampa_E * 1*ms/(sqrt(Kappa_E)*tau_Ampa) : ampere / meter**2
ds_Ampa_E / dt = - s_Ampa_E / tau_Ampa : 1
g_syn_E: siemens / meter**2
Kappa_E :1
''' 
eqs_I_Gaba_I = '''
I_Gaba_I = g_syn_I * (v - E_Gaba) * s_Gaba_I * 1*ms/(sqrt(Kappa_I)*tau_Gaba): ampere / meter**2
ds_Gaba_I / dt = - s_Gaba_I / tau_Gaba : 1
g_syn_I: siemens / meter**2
Kappa_I :1
''' 
eqs_I_Gaba_R = '''
I_Gaba_R = g_syn_R * (v - E_Gaba) * s_Gaba_R * 1*ms/(sqrt(Kappa_R)*tau_Gaba): ampere / meter**2
ds_Gaba_R / dt = - s_Gaba_R / tau_Gaba : 1
g_syn_R: siemens / meter**2
Kappa_R :1
'''

#Synaptic Equations, separate for each synapse group

## Thalamus, Reticular Nrns, L4e,  L4i
Thal = NeuronGroup(N_t, 
                   model = eqs_aeIF +
                   eqs_I_Ampa_E +
                   eqs_I_Gaba_R + 
                   'I_syn = (I_Ampa_E + I_Gaba_R): ampere / meter**2',
                   method = 'euler', threshold = "v>=Vt",reset = 'v=Vr; w+=b', refractory = 2.5*ms)
Thal.a = 0.04 *uS #Destexhe 2009
Thal.b = 0.0 * nA
Thal.v = Vr  
Thal.Rad = R_Th

L4e = NeuronGroup(N_e,
                  model = eqs_aeIF +
                  eqs_I_Ampa_T + 
                  eqs_I_Ampa_E +
                  eqs_I_Gaba_I + 
                  'I_syn = (I_Ampa_T + I_Ampa_E + I_Gaba_I): ampere / meter**2',
                  method = 'euler', threshold = "v>=Vt",reset = 'v=Vr; w+=b',refractory = 2.5*ms)
L4e.a = 0.001 * uS
L4e.b = 0.04 * nA
L4e.v = Vr
L4e.Rad = R_Cx

L4i = NeuronGroup(N_i, 
                  model = eqs_aeIF +
                  eqs_I_Ampa_T + 
                  eqs_I_Ampa_E +
                  eqs_I_Gaba_I + 
                  'I_syn = (I_Ampa_T + I_Ampa_E + I_Gaba_I): ampere / meter**2', 
                  method = 'euler',threshold = "v>=Vt",reset = 'v=Vr; w+=b',refractory = 2.5*ms)
L4i.a = 0.001 * uS
L4i.b = 0.0 * nA
L4i.v = Vr
L4i.Rad = R_Cx

Ret = NeuronGroup(N_r, 
                  model = eqs_aeIF +
                  eqs_I_Ampa_T + 
                  eqs_I_Ampa_E +
                  eqs_I_Gaba_R + 
                  'I_syn = (I_Ampa_T + I_Ampa_E + I_Gaba_R): ampere / meter**2', 
                  method = 'euler',threshold = "v>=Vt",reset = 'v=Vr; w+=b',refractory = 2.5*ms)
Ret.a = 0.08 * uS
Ret.b = 0.03 * nA
Ret.v = Vr
Ret.Rad = R_Th 

Thal.Kappa_R = Kappa_t_r
Thal.Kappa_E = Kappa_t_e
Thal.g_syn_R = g_T_R
Thal.g_syn_E = g_T_E

Ret.Kappa_T = Kappa_r_t
Ret.Kappa_R = Kappa_r_r
Ret.Kappa_E = Kappa_r_e
Ret.g_syn_T = g_R_T
Ret.g_syn_R = g_R_R
Ret.g_syn_E = g_R_E

L4e.Kappa_T  = Kappa_e_t
L4i.Kappa_T  = Kappa_i_t
L4e.Kappa_E  = Kappa_e_e
L4i.Kappa_E  = Kappa_i_e
L4e.Kappa_I  = Kappa_e_i
L4i.Kappa_I  = Kappa_i_i

L4e.g_syn_T = g_E_T
L4i.g_syn_T = g_I_T
L4e.g_syn_E = g_E_E
L4i.g_syn_E = g_I_E
L4e.g_syn_I = g_E_I
L4i.g_syn_I = g_I_I


P_t_exc = PoissonInput(Thal, 'v', Po_N_exc, Po_Rate, weight=Po_Amp_T)
P_i_exc = PoissonInput(L4i,  'v', Po_N_exc, Po_Rate, weight=Po_Amp_I)
P_e_exc = PoissonInput(L4e,  'v', Po_N_exc, Po_Rate, weight=Po_Amp_E)
P_r_exc = PoissonInput(Ret, 'v', Po_N_exc, Po_Rate, weight=Po_Amp_R)

P_t_inh = PoissonInput(Thal, 'v', Po_N_inh, Po_Rate, weight=-Po_Amp_T) 
P_i_inh = PoissonInput(L4i,  'v', Po_N_inh, Po_Rate, weight=-Po_Amp_I)
P_e_inh = PoissonInput(L4e,  'v', Po_N_inh, Po_Rate, weight=-Po_Amp_E)
P_r_inh = PoissonInput(Ret, 'v', Po_N_inh, Po_Rate, weight=-Po_Amp_R)

#%% Connectivity and synapses


#read as: connect thalamus to l4e (no connectivity spec = fully connected)
S_e_t = Synapses(Thal, L4e, on_pre='s_Ampa_T+=1', delay = 1*ms)
S_i_t = Synapses(Thal, L4i, on_pre='s_Ampa_T+=1', delay = 1*ms)
S_e_e = Synapses(L4e, L4e,  on_pre='s_Ampa_E+=1', delay = 1*ms)
S_i_e = Synapses(L4e, L4i,  on_pre='s_Ampa_E+=1', delay = 1*ms)
S_e_i = Synapses(L4i, L4e,  on_pre='s_Gaba_I+=1', delay = .85*ms)
S_i_i = Synapses(L4i, L4i,  on_pre='s_Gaba_I+=1', delay = 0.5*ms)

S_t_e = Synapses(L4e, Thal, on_pre='s_Ampa_E+=1', delay = 1*ms)
S_r_e = Synapses(L4e,  Ret, on_pre='s_Ampa_E+=1', delay = 1*ms)
S_r_t = Synapses(Thal, Ret, on_pre='s_Ampa_T+=1', delay = 1*ms)
S_t_r = Synapses(Ret, Thal, on_pre='s_Gaba_R+=1', delay = 0.5*ms)
S_r_r = Synapses(Ret,  Ret, on_pre='s_Gaba_R+=1', delay = 0.5*ms)


S_e_t.connect(p=Kappa_e_t/N_t)  # 200 T -> 1600 E (25%)
S_i_t.connect(p=Kappa_i_t/N_t)
S_e_e.connect(p=Kappa_e_e/N_e)
S_i_e.connect(p=Kappa_i_e/N_e)
S_e_i.connect(p=Kappa_e_i/N_i)
S_i_i.connect(p=Kappa_i_i/N_i)

S_t_e.connect(p=Kappa_t_e/N_e)  # 1600 E -> 200 T (2%)
S_r_e.connect(p=Kappa_r_e/N_e)
S_r_t.connect(p=Kappa_r_t/N_t)
S_t_r.connect(p=Kappa_t_r/N_r)
S_r_r.connect(p=Kappa_r_r/N_r)
4

Hi @cwalsh27, thanks for the additional code. I’m afraid I still cannot run it though, the code you posted first refers to variables i and s that are not set here – I guess you are running things in a loop? Please provide some example parameters. If we do not need to run the same loop as you and a single parameter setting reproduces the issue, even better!

5

Hi @mstimberg-- apologies for the confusion, here is the full loop for the model runs. In this context, i and s are dictating whether stimulation is applied to the thalamic or cortical groups (s) and how much stimulation is applied through Istim (i). That said, only one iteration of the loop is necessary to reproduce the disconnect in network noise that I demonstrated in my original post’s figures.

# Run up to and including this cell to run the loop
StartTime = time.perf_counter()
for i, s, trial in product(Iapp, range(2), range(trials)):
    if trial == 0:
        print(time.perf_counter()-StartTime,
              '\n StimCx: ', s, 
              '\n Istim: ', i,'\n')
    restore()

    del Spike_t, Spike_e, Spike_i, Spike_r
    Spike_t = SpikeMonitor(Thal)
    Spike_e = SpikeMonitor(L4e)
    Spike_i = SpikeMonitor(L4i)
    Spike_r = SpikeMonitor(Ret)
    setParam('Istim', str(i) + '*pA', params)
    setParam('StimCx', str(s), params)

    Thal.I_app = Istim * (1-StimCx)
    L4e.I_app = Istim * StimCx
    L4i.I_app = Istim * StimCx

    run(run_time)

    spikes_T = Spike_t.get_states(['t', 'i'], units = False, format = 'pandas')
    spikes_E = Spike_e.get_states(['t', 'i'], units = False, format = 'pandas')
    spikes_I = Spike_i.get_states(['t', 'i'], units = False, format = 'pandas')
    spikes_R = Spike_r.get_states(['t', 'i'], units = False, format = 'pandas')
6

Hi again. This is still missing the variable Iapp to actually run, but I now realized that your initial plot had parameter values in the title, and that you showed your issue with Istim=0*pA and StimCx=0, so this is what I tried. My full code is here: brian_bug_discussion.py · GitHub

When I run it, I do not see anything particular at the beginning of the run-time period, but also the general activity looks quite different from what you have (more oscillations, less “noise”):

Figure_1
Figure_11000×400 21.4 KB

Do you get the same when you run this code? If yes, what is the difference to the code where you are seeing the issue?