Synapses cannot be simulated after using `set_states`

[Jasmine9691]

Description of problem

Hi Brian2 developers and users,

I used Network.set_states to restart simulation with the file output by Network.get_states. The neurons can be simulated as expected, but the synaptic activity cannot: when the membrane potential of the presynaptic neuron exceeds the threshold, brian2 did not detect this event and the synaptic conductance did not change.

NOTE: store and restore cannot be used in my case because I’m computing on the server and post-processes results on my PC, which are different platforms. I could ONLY use set_states and get_states here.

Brian2 version: 2.9.0.
System: Both CentOS 7 and Windows 10.

Minimal code to reproduce problem

from brian2.only import *
import brian2.numpy_ as np
import matplotlib.pyplot as plt
import pickle

eqs = '''
dv/dt = (v0 - v) / tau : volt
tau : second
v0 : volt
'''
group = NeuronGroup(2, eqs, threshold='v > 10*mV', reset='v = 0*mV',
                    refractory=5 * ms, method='exact')
group.v = 0 * mV
group.v0 = np.array([20, 0]) * mV
group.tau = np.array([20, 200]) * ms

S = Synapses(group, group, on_pre='v_post += 1 * mV')
S.connect(i=0, j=1)
mon = StateMonitor(group, 'v', record=True)
net = Network(group, S, mon)
# =============================== save states
# net.run(100 * ms)
# with open('state1.pkl', 'wb') as pf:
#     pickle.dump(net.get_states(read_only_variables=True), pf)

# =============================== load states
with open('state1.pkl', 'rb') as pf:
    states = pickle.load(pf)
    net.set_states(states)
net.run(50 * ms)
plt.plot(mon.t / ms, mon.v.T / mV)
plt.xlabel('t (ms)')
plt.ylabel('v (mV)')
plt.show()

Expected output (if relevant)

The synaptic event can be detected.

Actual output (if relevant)

Neuron #1 did not reset its potential when v>10*mV

[mstimberg]

Hi @Jasmine9691. I tried to reproduce this issue with your code, but I wasn’t able to right away. I first run the code with the save_states part and then again with the load_states part and this led to an error since it cannot write read-only variables with set_state (the get_states stored e.g. N of the NeuronGroup). After changing read_only_variables to False it ran, and I can reproduce the problem now.
The problem is that when you store the states, you are also storing the time of the last spike for the refractoriness mechanism. Since the new simulation starts at 0s, the neuron is considered refractory given that lastspike + refractory is in the future. For this specific scenario, you can correct these values with something like:

    # Correct the lastspike... variables to not point into the future
    for objects in states.values():
        for k, v in objects.items():
            if k.startswith("lastspike"):
                v -= 100*ms

before using net.set_states.

But get_states will not store all the internal state of the system, in particular for random numbers and spikes in the queue that have not been delivered yet (particularly relevant if you have synaptic delays). Can you explain in more detail why you cannot use store and restore? It also uses pickle under the hood, and its file should be transferable in the same way as you do for your pickle file that you create from get_states.

[Jasmine9691]

Hi @mstimberg. Thanks for your kind advice.

After changing read_only_variables to False it ran, and I can reproduce the problem now.

Sorry about the typo. This is indeed False.

For this specific scenario, you can correct these values with something like:

The problem has been resolved with your code. For my simulation with more complicated neuron models, I resolved it with v-=1*ksecond (v-=100*ms cannot). However, for neurons that have long refractory periods, this remedy will be inapplicable/inaccurate because a neuron will fire spikes ahead of expected time moments. And I have two questions about it.

1. The neuron is only refractory from lastspike to lastspike + refractory. In the minimal working example, lastspike=97.2*ms, refractory=5*ms. So the neuron is only refractory from 97.2*ms to 102.2*ms. When I restarted the simulation, it is reasonable that the neuron didn’t fire spikes from 100*ms to 102.2*ms. But why couldn’t it fire after 102.2*ms, which has already passed the refractory period?

2. Modifying lastspike does not influence SpikeMonitor, does it?

Can you explain in more detail why you cannot use store and restore ?

Hmm, the reasons are rather complex and multi-faceted:

• set_states are very flexible for users. For example, I can modify some states like lastspike, and I can remove some states if I delete a variable (restore will report errors).
• I’m currently coupling brian2 with another library (LAMMPS) to obtain the stimulus to the neurons. LAMMPS is running using MPICH. Brian2 is running in the main process and LAMMPS is running in 24 processes. I did use store and restore before. However, weeks ago, I found that the file written by Network.store cannot be read normally on my PC (Windows 10) during post-processing (in serial). But the problem would be resolved if I ran the simulation in serial so that the stored file is obtained in the serial mode. Then I guess Brian2 is not compatible with MPICH. The error is as follows:

ERROR      Brian 2 encountered an unexpected error. If you think this is a bug in Brian 2, please report this issue either to the discourse forum at <http://brian.discourse.group/>, or to the issue tracker at <https://github.com/brian-team/brian2/issues>. Please include this file with debug information in your report: C:\Users\13372\AppData\Local\Temp\brian_debug_hend29pl.log  Additionally, you can also include a copy of the script that was run, available at: C:\Users\13372\AppData\Local\Temp\brian_script_t7rc8j5o.py Thanks! [brian2]
Traceback (most recent call last):
  File "E:\EntericNervousSystem\my_work\postprocess\visualize_network.py", line 24, in <module>
    net.restore(filename=f'{path}/{case_name}/store/net_{read_step}.store')
  File "D:\Anaconda3\envs\enteric_new\lib\site-packages\brian2\core\base.py", line 346, in device_override_decorated_function
    return func(*args, **kwds)
  File "D:\Anaconda3\envs\enteric_new\lib\site-packages\brian2\core\network.py", line 705, in restore
    obj._restore_from_full_state(state[obj.name])
  File "D:\Anaconda3\envs\enteric_new\lib\site-packages\brian2\synapses\synapses.py", line 452, in _restore_from_full_state
    self.queue._restore_from_full_state(queue_state)
  File "brian2\\synapses\\cythonspikequeue.pyx", line 63, in brian2.synapses.cythonspikequeue.SpikeQueue._restore_from_full_state
  File "<stringsource>", line 178, in pair.from_py.__pyx_convert_pair_from_py_int__and_std_3a__3a_vector_3c_std_3a__3a_vector_3c_int32_t_3e____3e___
ValueError: too many values to unpack (expected 2)

I did not report this error in this Brian2 community at that time because I think it might be more flexible if I used get_states and set_states. Also, I reckoned that brian2 might give low priority to this issue because most users would choose multiprocessing for parallel simulation. The reason why LAMMPS chooses MPI is that it is suitable for space decomposition, which is not applicable to computational neuroscience. Anyway, I switched to set_states finally.

[mstimberg]

@Jasmine9691

The problem has been resolved with your code. For my simulation with more complicated neuron models, I resolved it with v-=1*ksecond (v-=100*ms cannot). However, for neurons that have long refractory periods, this remedy will be inapplicable/inaccurate because a neuron will fire spikes ahead of expected time moments.

The time you subtract should correspond to the time at the end of the previous run (100ms in your example code). E.g. if a neuron with a 5ms refractory period spikes at 99ms and you save its state after 100ms, you want the lastspike variable to be at -1ms for the second simulation, so that it stays refractory until 4ms.

@Jasmine9691

The neuron is only refractory from lastspike to lastspike + refractory. In the minimal working example, lastspike=97.2*ms, refractory=5*ms. So the neuron is only refractory from 97.2*ms to 102.2*ms. When I restarted the simulation, it is reasonable that the neuron didn’t fire spikes from 100*ms to 102.2*ms. But why couldn’t it fire after 102.2*ms, which has already passed the refractory period?

The time in the second simulation starts at 0ms, and it will not spike until lastspike + refractory, i.e. for the first 102.2ms of the simulation. In your minimal example from above, you only run the second simulation for 50ms, so it never gets to the point where it spikes again. If you run it for e.g. 150ms, you see this effect quite clearly:

@Jasmine9691

Modifying lastspike does not influence SpikeMonitor, does it?

The lastspike variable is only used for the refractoriness mechanism, so it does not affect SpikeMonitor directly (of course, it affects it indirectly by preventing the neuron from spiking…).

@Jasmine9691

set_states are very flexible for users. For example, I can modify some states like lastspike, and I can remove some states if I delete a variable (restore will report errors).

Yes, but note that you can also access and modify variables before/after you restore, giving you the same flexibility in that regard. It is true that store/restore assumes an unchanged network, but you can always use it on a full network and then set objects you do not want to simulate as inactive (by setting their active attribute to False) instead of removing them from the network.

Regarding the restore error: from the message, I think the issue is not with MPI, but I rather suspect that the machine that called Network.storedid not use the optimized Cython version of the spike queue and fell back to the Python implementation (there should have been a warning about this). This has less good performance, but also uses a slightly different format when storing the internal state to disk (it probably shouldn’t!), leading to the error you saw. These days, with Python wheels and conda packages, it is rather unusual to not have the Cython spike queue available. Is it possible that you were running Brian 2 directly from the source directory without installation, or that you did an installation from source without Cython and/or a C++ compiler available?