Patch to allow restore after recreating the geometry

[jedlund]

When running multiple simulations on the same geometry (with different extracellular fields and/or synapse weights, etc) it's useful to be able to run the simulation to steady state, save the state of the sections, and then use this saved state to start new simulations after recreating the exact same geometry (in python). In order to do that I made the following small change to savstate.cpp. I thought that this might also be useful to other users. I can also create a simple python program demonstrating it's use if there is interest.

Code: Select all

*** nrn-7.4-orig/src/nrniv/savstate.cpp 2015-04-23 14:06:34.000000000 -0700
--- nrn-7.4/src/nrniv/savstate.cpp      2015-07-18 11:24:21.848503407 -0700
***************
*** 275,281 ****
                }
                return false;
        }
!       if (nsec_ && ss_[0].sec == NULL) { // got the data from a read
                isec = 0; ForAllSections(sec)
                        ss_[isec].sec = sec;
                        section_ref(ss_[isec].sec);
--- 275,282 ----
                }
                return false;
        }
!       if (nsec_ && ((ss_[0].sec == NULL)||(!ss_[0].sec->prop))) { // got the data from a read (or recreated the geometry)
!         fprintf(stderr, "SaveState: Setting sections.\n");
                isec = 0; ForAllSections(sec)
                        ss_[isec].sec = sec;
                        section_ref(ss_[isec].sec);

[ted]

We're all abuzz with curiosity about this. I for one haven't run into a problem using SaveState via hoc, nor have I heard of any problems using it through Python. A short program that demonstrates the utility of this patch would be very interesting. If the demo is more than roughly 50 lines long, perhaps it would be best to email it to
robert dot mcdougal at yale dot edu

[jedlund]

This patch might not be as useful as I thought. Here's the basic outline of what I was attempting:

1. Create neuron geometry (without extracellular fields, iclamps, or synapses)
2. Run the simulation to a steady state
3. Save the state in memory
4. Destroy the neuron geometry
5. Create neuron geometry (without extracellular fields, iclamps, or synapses)
6. Restore the state saved in step 3
7. Add the synapses, extracellular fields, or iclamps (depending on the simulation)
8. Run the simulation and save results
9. destroy the neuron geometry
10. go to step 5 unless finished.

Neuron would complain at step 6 (the restore step) because the original sections didn't exist and the restore wasn't from a file. I thought this was a silly requirement given the later checks on each section. Unfortunately perhaps it wasn't so silly. The above works for extracellular fields and iclamps, but I can't get it to work for ExpSyn. If I do it as above the synapse never fires (because the vecplay never delivers the NetCon event). If I set it up before the save and again (the same way) before the restore I get a core dump.

Code: Select all

(gdb) bt
#0  0x00007f5f6b4a18d1 in element () from x86_64/.libs/libnrnmech.so
#1  0x00007f5f6b4a17d6 in _net_receive () from x86_64/.libs/libnrnmech.so
#2  0x00007f5f6e34edcb in SelfEvent::call_net_receive(NetCvode*) () from /home/jedlund/GitReps/ComsolSpine/ComsolSpineVirtEnv/x86_64/lib/libnrniv.so.0
#3  0x00007f5f6e34ed28 in SelfEvent::deliver(double, NetCvode*, NrnThread*) () from /home/jedlund/GitReps/ComsolSpine/ComsolSpineVirtEnv/x86_64/lib/libnrniv.so.0
#4  0x00007f5f6e34ad1f in NetCvode::deliver_event(double, NrnThread*) () from /home/jedlund/GitReps/ComsolSpine/ComsolSpineVirtEnv/x86_64/lib/libnrniv.so.0

I'm going to have to figure out a better way to do this. This there any way to setup a synapse (using vecplay on a NetCon object) in the middle of a simulation? I know that vector play needs to be setup before neuron.h.finitialize(self.v_init). Is there a way around that?

The code is about 400 lines long (with tests) (I can post it or email it if desired). Below is just the section where I setup the synapse. Note that there is only a single passive section in the test model.

Code: Select all

    def SetupSynapse(self, syn_weight, syn_x, syn_time_ms):
        print "SetupSynapse"
        self.synInfo = {}

        if syn_weight is not None and syn_x is not None and syn_time_ms is not None:
            secIndex = 0
            synTime_ms = self.init_ms + syn_time_ms
            # Now we know secIndex and segx.
            synObj = neuron.h.ExpSyn(self.allSections[secIndex](syn_x))
            synObj.e = 100
            synObj.tau = 1
            self.synInfo['synObj'] = synObj
            self.synInfo['vecStim'] = neuron.h.VecStim()  # Requires the vecevent.mod file compiled and loaded.
            self.synInfo['spikeTimeVec'] = neuron.h.Vector([synTime_ms])
            self.synInfo['vecStim'].play(self.synInfo['spikeTimeVec'])
            threshold = 10  # The default
            delay = 0  # Instead of 1ms.
            weight = 0  # The default
            self.synInfo['netCon'] = neuron.h.NetCon(self.synInfo['vecStim'], self.synInfo["synObj"], threshold, delay, weight, sec=self.allSections[secIndex])
            self.synInfo['netCon'].weight[0] = syn_weight
            pass
        pass

Thanks.

[jedlund]

After having problems with the above patch and ExpSyn, I went back to the standard neuron 7.4 and tried:
1. create a geometry
2. install a synapse (set to fire at 410 ms)
3. ran the simulation for 400ms
4. save the state to disk
5. finish the simulation (synapse fires fine)
6. Steps 1 and 2 again
7. Restore the state from disk (appears to work)
8. finish the simulation (Results in the error: net_send td-t = -410 SelfEvent target=VecStim[0] 0 flag=1)

FYI: I put in extra print statements to determine that td=0 and t = 410. Any idea what's going on?

Here's the code to recreate it. Note that you need to have vecevent.mod compiled and available.

Code: Select all

#!/usr/bin/env python

import neuron
import nrn
import numpy


class SimpleModel:
    def __init__(self):
        neuron.h.celsius = 37
        self.nSegs = 10
        self.length = 300
        self.diam = 3
        self.init_ms = 400
        self.calc_ms = 100
        self.tStep_ms = 10e-3
        self.v_init = -72  # Off from e_pas to demo InitSimulation 

        self.saved_state_filename = False
        pass

    def Execute(self, init_only=False, skip_load=False, syn_weight=None, syn_time_ms=None, syn_x=None):
        print "Execute"
        self.CreateGeometry()
        self.SetupSignal()
        self.SetupNeuronRecord()
        self.SetupSynapse(syn_weight, syn_x, syn_time_ms)

        self.InitSimulation(skip_load=skip_load)
        self.CollectData()

        self.initOutput = self.output
        self.output = None

        if init_only:
            return

        self.RunSimulation()
        # print "Collecting data"
        self.CollectData()
        pass

    def SetupSignal(self, stim_scale=0):
        print "SetupSignal"
        self.initTime_ms = numpy.arange(0, self.init_ms, self.tStep_ms)
        self.initN = len(self.initTime_ms)
        self.stimTime_ms = numpy.arange(self.init_ms, self.init_ms+self.calc_ms, self.tStep_ms)
        self.calcN = len(self.stimTime_ms)
        pass

    def SetupSynapse(self, syn_weight, syn_x, syn_time_ms):
        print "SetupSynapse"
        self.synInfo = {}

        if syn_weight is not None and syn_x is not None and syn_time_ms is not None:
            secIndex = 0
            synTime_ms = self.init_ms + syn_time_ms
            # Now we know secIndex and segx.
            synObj = neuron.h.ExpSyn(self.allSections[secIndex](syn_x))
            synObj.e = 100
            synObj.tau = 1
            self.synInfo['synObj'] = synObj
            self.synInfo['vecStim'] = neuron.h.VecStim()  # Requires the vecevent.mod file compiled and loaded.
            self.synInfo['spikeTimeVec'] = neuron.h.Vector([synTime_ms])
            self.synInfo['vecStim'].play(self.synInfo['spikeTimeVec'])
            threshold = 10  # The default
            delay = 0  # Instead of 1ms.
            weight = 0  # The default
            self.synInfo['netCon'] = neuron.h.NetCon(self.synInfo['vecStim'], self.synInfo["synObj"], threshold, delay, weight, sec=self.allSections[secIndex])
            self.synInfo['netCon'].weight[0] = syn_weight
            pass
        pass

    def CreateGeometry(self):
        print "CreateGeometry"
        '''Build the NEURON Model'''
        self.allSections = []
        # nodeMyelinLength_um = self.nodeLength_um + self.myelinLength_um
        sec = nrn.Section(name='OnlySection')
        sec.L = self.length
        sec.diam = self.diam
        sec.Ra = 87.0  # ohm cm.                    # From Ostroumov 2007 section 3.6 <- Thurbon 1998
        sec.cm = 2.4  # uF/cm2 ,                # From Ostroumov 2007 section 3.6 <- Thurbon 1998
        sec.insert('pas')
        sec.e_pas = -70 # mV                   # From Ostroumov 2007 section 3.6 <- Thurbon 1998
        sec.g_pas = 1/5300.0

        sec.insert('extracellular')

        sec.nseg = self.nSegs

        neuron.h.area(0.5, sec=sec)
        self.segX = [seg.x for seg in sec]
        self.allSections.append(sec)
        pass

    def SetupNeuronRecord(self):
        print "SetupNeuronRecord"
        # Record Time from NEURON (neuron.h._ref_t)
        self.record_time = neuron.h.Vector()
        self.record_time.record(neuron.h._ref_t)

        # Record Voltage from everywhere along the axon
        self.recordVoltVecs = []

        for iSec, sec in enumerate(self.allSections):
            # neuron.h.psection(sec=sec)

            for seg in sec:
                # Setup record for voltage
                record_voltage = neuron.h.Vector()
                record_voltage.record(sec(seg.x)._ref_v)
                self.recordVoltVecs.append(record_voltage)
                pass
            pass
        pass

    def InitSimulation(self, skip_load=False):
        print "InitSimulation"
        neuron.h('dt = '+str(self.tStep_ms))
        # finitialize sets v to v_init and m,h,n to corresponding
        # steady state values. Also delivers vector play. Can't play
        # more vectors after it.
        neuron.h.finitialize(self.v_init)  # finitialize counts as a step.
        neuron.h.fcurrent()

        if not skip_load and self.saved_state_filename:
            # ipython_keyboard()
            print "Restoring SaveState"
            self.RestoreSaveStateFromFile(self.saved_state_filename)
            print "Restored SaveState"
            #ipython_keyboard()
            pass
        else:
            # TODO: consider using Cython for this loop
            for tIndex in range(1, self.initN+1):  # I believe finitialize counts as a step.
                neuron.h.fadvance()
                assert abs(neuron.h.t - self.record_time[-1]) < self.tStep_ms/2
                assert abs(self.record_time[-1] - tIndex*self.tStep_ms) < self.tStep_ms/2
                pass
            self.SaveStateToFile(filename=self.saved_state_filename)
            pass
        assert abs(neuron.h.t - self.init_ms) < self.tStep_ms/2
        # assert abs(self.record_time[-1] - self.init_ms) < self.tStep_ms/2 # record_time isn't correct until frecord_init
        pass

    def RunSimulation(self):
        print "RunSimulation"
        # neuron.h('dt = '+str(self.tStep_ms))
        neuron.h.frecord_init()
        assert abs(neuron.h.t - self.init_ms) < self.tStep_ms/2
        assert abs(self.record_time[-1] - self.init_ms) < self.tStep_ms/2
        for tIndex in range(1, self.calcN+1):  # I believe finitialize counts as a step.
            #print neuron.h.t
            neuron.h.fadvance()
            assert abs(neuron.h.t - self.record_time[-1]) < self.tStep_ms/2
            assert abs(self.record_time[-1] - (tIndex+self.initN)*self.tStep_ms) < self.tStep_ms/2
            pass
        assert abs(neuron.h.t - (self.init_ms+self.calc_ms)) < self.tStep_ms/2
        assert abs(self.record_time[-1] - (self.init_ms+self.calc_ms)) < self.tStep_ms/2
        pass

    def RestoreSaveState(self, save_state):
        print "RestoreSaveState"
        save_state.restore(1)
        pass

    def RestoreSaveStateFromFile(self, filename):
        print "RestoreSaveState"
        neuron_file = neuron.h.File()
        neuron_file.ropen(filename)
        self.save_state = neuron.h.SaveState()
        self.save_state.fread(neuron_file)
        neuron_file.close()
        self.save_state.restore(1)
        pass

    def SaveStateToFile(self, filename=None):
        if not filename:
            import os
            import tempfile
            (fd, filename) = tempfile.mkstemp(suffix='.dat', prefix='neuronsavestate-', dir=os.curdir, text=False)
            os.close(fd)
            pass

        print "SaveState"
        self.save_state = neuron.h.SaveState()
        self.save_state.save()
        neuron_file = neuron.h.File()
        neuron_file.wopen(filename)
        self.save_state.fwrite(neuron_file)
        neuron_file.close()
        self.saved_state_filename = filename
        return self.save_state

    def SaveState(self):
        print "SaveState"
        self.save_state = neuron.h.SaveState()
        self.save_state.save()
        return self.save_state

    def GetSavedState(self):
        print "GetSavedState"
        return self.save_state

    def CollectData(self):
        print "CollectData"
        self.output = {}
        self.output["time"] = numpy.array(self.record_time)
        self.output["Vm"] = [numpy.array(vec) for vec in self.recordVoltVecs]
        pass


def ipython_keyboard():
    import sys
    import IPython
    # use exception trick to pick up the current frame
    try:
        raise None
    except:
        frame = sys.exc_info()[2].tb_frame.f_back
        pass
    # evaluate commands in current namespace
    namespace = frame.f_globals.copy()
    namespace.update(frame.f_locals)
    IPython.embed(user_ns=namespace)
    pass

if __name__ == "__main__":
    print "Compare results from a continuous run and restoring the init part (For ExpSyn) Requires the vecevent.mod file compiled and loaded."
    results = {}
    model = SimpleModel()
    model.Execute(syn_weight=3, syn_x=0.5, syn_time_ms=10)  # This runs all the way through
    results["no restore"] = model.output
    model.Execute(syn_weight=3, syn_x=0.5, syn_time_ms=10)  # Restoring the save_state and then running
    results["restored"] = model.output

    from matplotlib import pyplot
    for key, output in results.iteritems():
        pyplot.plot(output["time"], output["Vm"][3], label="Vm iSeg=3 ExpSyn key={}".format(key))
        pass
    pyplot.legend()
    pyplot.show()
    pass

Also, what's going on with neuron.h.CVode? A bunch of the forum posts use neuron.h.cvode.print_event_queue() or neuron.h.cvode.active(). neuron.h.cvode doesn't seem to exist in my builds of neuron 7.3 and 7.4. I guess there was a change in case at some point? Trying to call neuron.h.CVode.active() results in the message: AttributeError: 'hoc.HocObject' object has no attribute 'active'

Thanks.

Jeffrey

[jedlund]

I haven't solved the above problem with restoring synapses, but I did discover that you can install them after a restore by calling finitialize() after the restore. I developed 2 examples and emailed them to Robert.

SimpleModelSaveToFile.py demonstrates one way of running a geometry to steady state, saving the state of the neuron to disk, restoring the save state, and continuing the simulation (with optional additions like synapses, iclamps, extracellular fields).

SimpleModelWithPatch.py demonstrates the same thing with out saving/loading to/from disk. The patch is required to restore the savestate on the new geometry without loading from disk.

The basic order of execution:

To generate the save state:
1. Create the geometry
2. Call finitialize(v_init)
3. call fadvance() until the simulation reaches steady state
4. save the state (to memory or disk)

To continue the simulation:
1. Destroy any previous sections, synapses, iclamps (You'll likely have problems if you leave some.)
2. Create the geometry
3. restore the state (from memory or disk)
4. Add synapses, iclamps, extracellular fields
5. save t, call finitialize(), reset t
6. call frecord_init to reset recording vectors
7. call fadvance until the simulation is done.