Steady state init with voltage VS with SaveState

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

Post Reply
gabrielggn
Posts: 10
Joined: Fri May 20, 2022 8:07 am

Steady state init with voltage VS with SaveState

Post by gabrielggn »

Dear all,

I am working with pyramidal cell models and my pipeline is pretty simple:
  1. Compute steady state of the cell given channel mechanismes ( voltage of all compartments do not vary more than 1e-6 mV from previous step)
  • Preforming my modeling which involve synaptic inputs, calling h.finitialize() before h.continuerun(tstop)
I considered two way of doing so at the end (1) :
  • save the voltage of each seg of the cell in a list to further assign the voltage at these segments
  • save the state using the SaveState() functionality provided in NEURON

Using SaveState is not very practical if one wants to change something (netcon or others) this is why using voltages seems simple and more flexible.
However, I notice a difference using the two in the resulting voltage traces with the same input.
I also noticed that after using finitialize() and h.fcurrent() in both cases, the potential of the soma increases until another steady potential (from the previous steady state). I suspect a INITIAL block statement causing this but I am not sure.

I know I probably used a toy model first to test my pipeline. But more generally, do you know if there is a difference between initializing a steady state from voltages and SaveState ?
ted
Site Admin
Posts: 6305
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Steady state init with voltage VS with SaveState

Post by ted »

Compute steady state of the cell given channel mechanismes ( voltage of all compartments do not vary more than 1e-6 mV from previous step)
1. Experimentalists typically wait until v near the spike trigger zone is flat for an amount of time on the order of at least several seconds if not minutes. And even then, they often want each recording to start with at least 1 ms during which that v remains flat before applying any stimuli or injecting current pulses.
2. v is the only state variable in your model that has dynamics? There are no ion accumulation mechanisms (whose NEURON blocks contain one or more USEION statement that WRITEs a concentration), and no second messengers or other signals that are governed by ODEs or kinetic schemes?
Preforming my modeling which involve synaptic inputs, calling h.finitialize() before h.continuerun(tstop)
So you're just doing the typical "voltage clamp initialization" in which each state variable is only dependent on local membrane potential, and is initialized to the value it would have if local membrane potential had been held at some specific level? Which in turn means that you're not dealing with spontaneously active cells, or initializing cells to a "dynamic steady state" (as would happen if a cell is embedded in a milieu of ongoing synaptic activation).
I considered two way of doing so at the end (1) :
  • save the voltage of each seg of the cell in a list to further assign the voltage at these segments
  • save the state using the SaveState() functionality provided in NEURON
. . .
However, I notice a difference using the two in the resulting voltage traces with the same input.
I also noticed that after using finitialize() and h.fcurrent() in both cases, the potential of the soma increases until another steady potential (from the previous steady state).
I'm not sure what you're describing, but it doesn't sound good.
I suspect a INITIAL block statement causing this
Maybe, if one or more of your mechanisms involve user-written NMDOL code.
do you know if there is a difference between initializing a steady state from voltages and SaveState?
The question hasn't arisen before in my experience. To say more, and especially to come up with a solution, I'd have to work with the code that reproduces the symptom(s). If you have a small example that demonstrates what you're seeing, zip it up including all .py, .hoc, .ses, .mod and other source code files that are needed, and email it to ted dot carnevale at yale dot edu,
and I'll tell you what I find out.
gabrielggn
Posts: 10
Joined: Fri May 20, 2022 8:07 am

Re: Steady state init with voltage VS with SaveState

Post by gabrielggn »

Dear Ted,

Thank you for your quick feedback.
1. Experimentalists typically wait until v near the spike trigger zone is flat for an amount of time on the order of at least several seconds if not minutes. And even then, they often want each recording to start with at least 1 ms during which that v remains flat before applying any stimuli or injecting current pulses.
2. v is the only state variable in your model that has dynamics? There are no ion accumulation mechanisms (whose NEURON blocks contain one or more USEION statement that WRITEs a concentration), and no second messengers or other signals that are governed by ODEs or kinetic schemes?
My cell indeed contains mechanisms using USION statements. Your first two statements are pointing toward the fact that I should run the simulation without any input for several seconds to get a steady state for voltage but also for all ionics distributions to be at equilibrium, isn't it ? I though that unstable ions concentrations would be reflected by the voltage accross the membrane as stated in the GHK equation. But know that you point this, there is not a reprocity between voltages and ions concentrations, is it what I was missing ?
So you're just doing the typical "voltage clamp initialization" in which each state variable is only dependent on local membrane potential, and is initialized to the value it would have if local membrane potential had been held at some specific level? Which in turn means that you're not dealing with spontaneously active cells, or initializing cells to a "dynamic steady state" (as would happen if a cell is embedded in a milieu of ongoing synaptic activation).
Indeed this is what I'm doing right now. Should it be preferable to have a "dynamic steady state" if I want to see the effect of a perturabation ? I'm studying extracellular stimulation and I performed first the steady state (ss) as previously expained, then 2 minutes of activity with synaptic inputs and finally N sec of extra stimulation.
I do not have any knowledge of dynamic steady state, is it only doing a pre-simulation with normal synaptic input for few seconds before the simulation of interest ?

I was thiking that maybe I could write easily a script to save the density_mechs and ions parameters given by each sections calling psection() to insert back each according section. I did not find the source code of SaveState but I guess it is a part of what is done when it is called.
It is also something that I could use next for example if I want to compare activity off and on stimualtion using the same inputs. Indeed, the pre-activity block of few sec takes hours to run and could be squeezed by restoring the state at the time the stimulation is turn on. However, using SaveState.restore() with all my attemps, I was unable to simulate the extracellular stimulation after restoring the state and it is not very flexible with new events.
I'm sorry if my post is heavy but the final point is: Do you think it would be helpful for me to implement a function to set back all mechanisms and ions parameters (and probably parameters of synapses) using values given by psection and should it give the same final results as SaveState.restore() with more flexibility ?

Thanks for your time and help
gabrielggn
Posts: 10
Joined: Fri May 20, 2022 8:07 am

Re: Steady state init with voltage VS with SaveState

Post by gabrielggn »

Finally, I did as explain in the last post using psection() function.
The code to save the membrane/mechanisms+ions state was the following (using python) :

Code: Select all

    def save_state(self, fn = None):
        l = []

        for sec in self.cell.all:
            dic = sec.psection()
            v = np.array([ seg.v for seg in sec])
            l.append( {"v":v, 
                       "density_mechs": dic["density_mechs"], 
                       "ions": dic["ions"]} )
                       
        with open(fn, "wb") as f:
            pkl.dump(l, f)

    ## END save_state()
After saving the desired state, here the steady state after 5 seconds of simulation with any input, I restore the state with :

Code: Select all

    def restore_fromdict(self, d):
        if type(d) == str:
            with open(d, "rb") as f:
                d = pkl.load(f)
        try:
            for i,sec in enumerate(self.cell.all):
                for j, seg in enumerate(sec):
                    seg.v = d[i]["v"][j]

                    ## setting density mechanisms with taking car to extracellular as it is calling differently
                    mechanisms = d[i]["density_mechs"]
                    for m in mechanisms:
                        if m != "extracellular":
                            for sub_m in mechanisms[m]:
                                exec( f"seg.{m}.{sub_m} = {mechanisms[m][sub_m][j]}"  )
                        else:
                            exec( f"seg.vext[0] = {mechanisms[m]['vext'][j][0]}"  )
                            exec( f"seg.vext[1] = {mechanisms[m]['vext'][j][1]}"  )
                            exec( f"seg.i_membrane = {mechanisms[m]['i_membrane'][j]}"  )
                            
                    ## setting ions concentrations
                    ions = d[i]["ions"]
                    for ion in ions:
                        for ci in ions[ion]:
                            exec( f"seg.{ion}_ion.{ci} = {ions[ion][ci][j]}"  )

        except:
            print("Could not affect density mechs and ions to each seg")
            pass
And it works for the steady state (membrane potential does not vary after init and calling restore_fromdict().
Since it works that well and that was easy, I did the same for point_processes which are here only synapses (arround 5k) and I was planning to use this method to take back the simulation at time tstop1 to continue the stimulation to tstop2 (In my case to get ride of the pre-stimulation period of 2 minutes).
I try this with a prerun of 50ms and then try to get back the state however the time course (and AP events) are similar but not the same (same pre synaptic event). I was wondering if this could be done because of some hiden variable that I can't access from python to store/set back into the point_processes. For example I could not get all the ASSIGNED variables, do you know if it is possible to access these var ? Cause I saw in SaveState that it can save hiden variable and I was wondering if there is a play arround to get the same.

Best
ted
Site Admin
Posts: 6305
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Steady state init with voltage VS with SaveState

Post by ted »

There are no ion accumulation mechanisms (whose NEURON blocks contain one or more USEION statement that WRITEs a concentration), and no second messengers or other signals that are governed by ODEs or kinetic schemes?
My cell indeed contains mechanisms using USION statements.
But do any of them WRITE a concentration? e.g. a calcium accumulation mechanism would have a NEURON block that contains a WRITE statement of the form
USEION ca READ cai, ica WRITE cai
Your first two statements are pointing toward the fact that I should run the simulation without any input for several seconds to get a steady state for voltage but also for all ionics distributions to be at equilibrium
True, but remember that (1) models can contain many other state variables (e.g. channel gating states, second messengers other than Ca), (2) mechanisms may interact with each other in ways that make steady state initialization difficult or even impossible (consider models of cells that fire spontaneously).
I though that unstable ions concentrations would be reflected by the voltage accross the membrane
Not necessarily. Read enough modeling papers and you will discover that modelers have made many different assumptions about ionic currents, concentrations, and equilibrium potentials. It is very common to see that ionic currents are assumed to have no effect on concentrations or equilibrium potentials. In some models, at least one species of ionic current affects concentration, but not equilibrium potential. In others, some currents affect concentration and equilibrium potential, but others don't. The authors of these models all had reasonable justifications for their assumptions. To maximize the utility of NEURON, it was necessary to be able to accommodate all these different strategies. In most cases, NEURON can automatically determine what should happen based on the details of the mechanisms that are present in a section. However, in some cases, the user must be able to override this, and that's the purpose of the ion_style() function. For more information, read the documentation of the ion_style() function in the Programmer's Reference https://nrn.readthedocs.io/en/latest/py ... #ion_style
So you're just doing the typical "voltage clamp initialization" in which each state variable is only dependent on local membrane potential, and is initialized to the value it would have if local membrane potential had been held at some specific level? Which in turn means that you're not dealing with spontaneously active cells, or initializing cells to a "dynamic steady state" (as would happen if a cell is embedded in a milieu of ongoing synaptic activation).
Indeed this is what I'm doing right now. Should it be preferable to have a "dynamic steady state" if I want to see the effect of a perturabation ? I'm studying extracellular stimulation and I performed first the steady state (ss) as previously expained, then 2 minutes of activity with synaptic inputs and finally N sec of extra stimulation.
So your model cell is subjected to "background" synaptic activity for some time, and then you add the extracellular stimulus. Is the background synaptic activity stochastic, and does it differ from run to run?
I do not have any knowledge of dynamic steady state, is it only doing a pre-simulation with normal synaptic input for few seconds before the simulation of interest ?
Maybe a better term would be "perturbed steady state" which suggests a system that is subjected to perturbations yet remains "close" to an unperturbed steady state.
I was thiking that maybe I could write easily a script to save . . .
That's asking for trouble.
the pre-activity block of few sec takes hours to run and could be squeezed by restoring the state at the time the stimulation is turn on.
I think your best bet is

For the moment, forget about extracellular stimulation and even forget about your existing model cell.
Instead, make a toy model (a ball and stick model would be good), attach a couple of synaptic mechanisms to it, and run a simulation that lasts maybe 100 ms of model time and generates reproducible results (that is, each run produces the same result).

Then make a copy of this program, and add SaveState functionality to the copy. The new program should stop execution at 100 ms, and after the end of a simulation it should execute a procedure that calls SaveState.save() and then writes the saved states to a file via SaveState.fwrite().

Make a copy of your second program, and add a conditional statement that controls whether the program {does a 100 ms warmup run and calls the procedure that saves states to a file} or {calls a new procedure that uses SaveState.fread() to get states from a file, then uses SaveState.restore() to restore the states, and finally runs a simulation that stops when t is 200 ms}

Verify that your third program works properly by comparing its output to what your first program does when your first program stops at t == 200 ms; it should match the last 100 ms of your first program's output.

At this point your third program will be organized something like this:

1. declarations of important symbolic constants, e.g. you might use a variable called "runfromsaved" to control program execution
2. model setup code
3. def of a procedure that runs a simulation for 100 ms, then saves states and writes to a file
4. def of a procedure that reads states from a file, then restores states and runs a simulation for 100 ms of model time
5. a conditional statement that does
if runfromsaved is FALSE then execute the procedure that runs a simulation and saves states to a file
else execute the procedure that reads states from a file, restores them, and runs a simulation for 100 ms

Now you're ready to make a fourth program. The fourth program differs from the third one in one way: when runfromsaved is not FALSE, instead of doing this

else execute the procedure that reads states from a file, restores them, and runs a simulation for 100 ms

it does

else execute the procedure that reads states from a file, restores them, sets up extracellular stimulation, and runs a simulation for 100 ms

The only drawback of this final version of the program is that it can only be used to do a single run with extracellular stimulation. If you need to do multiple runs, maybe with different simulus amplitudes or patterns, you have to exit NEURON and then execute the program again. How can you tell NEURON to something different in each run? One way might be for your Python code read a parameter file. That could be done after states have been restored but before executing a simulation. Or maybe you could pass a value to Python via the command line (is that doable with Python? hoc allows it, but hoc allows all kinds of potentially problematical things . . .)
Post Reply