Phase shifting cell models with periodic spontaneous firing

[Helveg]

Is there a generic and straighforward way in NEURON to apply a distribution of the spike times to a population of identical cell models during one period? eg to make their spike times uniformly scattered across the period so that with a period of 1s 10% of the population has fired between 0 and 0.1s

I'd think not because it highly depends on all of the inner workings and states (voltages, currents, conductances, kinetic state and whatnot) of the sections and mechanisms. That's why I want to make a tool that will run a model, detect a period and completely scan the state of the cell along the period so that you get a discrete collection of states to use as initial states.

Now with the context out of the way, here's the issue: Is every piece of state in HOC necesarily exposed to Python? Are there local variables on mechanisms that I maybe can't read or write to? For example when I

>>> dir(s(0.5).pas)
['__class__', '__delattr__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__next__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', 'e', 'g', 'i', 'is_ion', 'name', 'segment']

Could I completely describe the state of pas with current values of e, g and i and could I force it to that state by setting e, g and i? Maybe pas isn't the greatest example because I suppose these values don't change over time and aren't based on other pieces of state. And would this always/mostly be true for anything I might encounter in the wild world of NEURON models?

[ramcdougal]

Yes, the full state is discoverable from both HOC and Python.

In Python, you get almost* all of a Section's information in the sec.psection() dictionary, where sec is a Section of interest. (* rxd module reactions are missing for some reason).

For example, for a single segment soma with both hh and pas inserted:

>>> pprint.pprint(soma.psection()['density_mechs'])
{'hh': {'el': [-54.3],
        'gk': [0.0],
        'gkbar': [0.036],
        'gl': [0.0003],
        'gna': [0.0],
        'gnabar': [0.12],
        'h': [0.0],
        'il': [0.0],
        'm': [0.0],
        'n': [0.0]},
 'pas': {'e': [-70.0], 'g': [0.001], 'i': [0.0]}}

Background comments:

HOC and Python are two semi-coupled interfaces to NEURON. Everything that can be done in HOC can be done directly in Python. Most (but not all) things that can be done in Python can be done in HOC.

The psection method code is written in Python and uses NEURON's lower-level MechanismStandard and MechanismType interfaces, so you could modify that if you wanted to save time. (These latter two classes can be used from both HOC and Python.)

[nrnhines]

I don't think the "entirely different solution" will work. It is not possible in NEURON to change the topology of a model, eg. number of cells or number of NetCon, without calling finitialize() to get all the internal structures setup again. I would consider creating one ARTIFICIAL_CELL for each population and connecting it to each cell of the population with a NetCon with a random delay (that you can change or re-randomize at desired intervals. But I may not completely understand your conceptual model. At least you can have a Python list of NetCon with
None source and connecting to the cell populations and feed spikes to them via NetCon.event(tdeliver)

Generally every variable in a mod file declared in PARAMETER, ASSIGNED, STATE blocks are available to Python (but ASSIGNED variables must also be listed in the NEURON RANGE statement.) You can change any of those between time steps with the fixed step method. But some variables may be computed from scratch during an integration step. More care needs to be taken with the variable step methods as any change to states or parameters is a new initial value problem which requires re-initialization of the variable step integrator (cvode.re_init())

[Helveg]

@ramcdougal Ok, great info, I'll see if I can mimick every step of the states with the psection dictionary.

@nrnhines I have different populations of cells that have a periodic autorhythm and if I start the simulation all cells of the same population all start in the same state and show synchronized activity:

A --|---|---|---|
B --|---|---|---|
C --|---|---|---|

Whereas I'd like their spike times to be randomly scattered across the period. (and I want to be able to give a phase probability distribution)

A -|---|---|---|--
B --|---|---|---|
C |---|---|---|---

I tried injecting some current at the start and it could be an OK solution but I'd have to find out a general solution to turn the amount of current to inject into a cell model into a linear relationship with the phaseshift. Similarly I tried to inject random numbers of spikes which also helps desync them from eachother but again: how many spikes will I need for any given collection of sections with any mechanisms inserted to phase shift them exactly as much as I want? :/

So then I was thinking that to propel a cell forward by a time t I'd just simulate it up to time t read its full state and in the simulation create it at time 0 with all its state variables set to what it would have at time t. That's the concept and why I wrote this question whether it's possible to completely determine and set the state of a cell.

At least you can have a Python list of NetCon with None source and connecting to the cell populations and feed spikes to them via NetCon.event(tdeliver)

If each cell gets spikes injected at different times from eachother they'd still just all be just 1 EPSC ahead of the normal period and their phase shifted by more or less the same amount

[nrnhines]

In your ABC diagrams you show each cell having identical interspike intervals. Is that just illustrative or are all the cells really that much the same and only the first spike obeys a phase probability distribution due to a different initial state for each cell? I likely am entirely missing critical aspects of your conceptual model but I wonder how much old state information survives a spike. I suppose that depends on how complicated the cells are. I would imagine one could always have a generalized stimulus that forces a spike initiation upon delivery of a spike initiation event (eg. sets voltage and certain gating states in the soma or initial segment of the cell)

[Helveg]

In your ABC diagrams you show each cell having identical interspike intervals. Is that just illustrative or are all the cells really that much the same

Yes each population of cells uses the same cell model (e.g. 100 identical clones that load the same morphology file, the same mechanisms are inserted in the same sections, the same parameters are set on these sections), but there are multiple different populations (basket cells, purkinje cells, stellate cells, golgi cells, ...). When they aren't connected to eachother all cells of the same population undergo exactly the same cycle. Because in the network they receive different inputs they'll eventually diverge from eachother, but I want to start my network simulation with cells that aren't synchronized in this way. And I'd like to solve the problem in general, not just by finding some set of stimuli that make specifically these populations diverge.

I would imagine one could always have a generalized stimulus that forces a spike initiation upon delivery of a spike initiation event (eg. sets voltage and certain gating states in the soma or initial segment of the cell)

Yes? Even without knowing anything about the cell model? I'd like a general solution where I can plug in any cell model with periodic activity. I tried for example setting random section.v for every section per cell but because of the initial state of the mechanisms and their time constants the voltage converges to the same point before the states of the cells are sufficiently altered from eachother causing only a slight shift in phase and it certainly doesn't give me the precision and generality that I'm looking for. Improving this technique could be using maybe a voltage clamp and clamping the voltage until all channels and currents have adjusted, but there is no guarantee that there's a linear relation between the soma membrane voltage and the phase (for example the Vm is -20mV on both the rising and falling edge of the action potential).

But if most variables that change over time are available through psection then I will explore this route first, thank you very much for lending your brainpower to the problem ;) I'll let you know if there are more developments

[Helveg]

@ramcdougal dir(section) does not seem to list all attributes of a section, I've discovered v and Ra, are there more? Why are they omitted from dir(section)?

[nrnhines]

dir(section) begins with ['L', 'Ra', but v is missing because it is the property of a nrn.Segment. I was wondering about rallbranch but I see it is there, so I am not aware of any missing names. @ramcdougal I see that v is missing from psection but states are listed if, e.g. hh is inserted.

[Helveg]

My bad, next thing, I'll open a new issue if it's actually an issue: I can only seem to read and not set values on mechanisms?

>>> bc.soma[0](0.5).ca_ion.cai
0.0061856020933762626
>>> bc.soma[0](0.5).ca_ion.cai = 5
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'nrn.Mechanism' object has no attribute 'cai'

[nrnhines]

That's very odd and needs to be fixed. For now one must bc.soma[0](0.5).cai = 5 which harks back to the hoc syntax where ion variable names did not require a suffix in contrast to names such as gnabar_hh Apparently when I updated getattr I forgot to update setattr for the ion special case. Go ahead and create the issue and I'll make a pull request responding to it.

[Helveg]

Ok, is there no danger of overlap between mechanisms this way?

[nrnhines]

To some extent. It is a good idea for mod files not to declare names with an underbar. All density mechanisms with a SUFFIX can use either section(x).suffix.var or section(x).var_suffix . But lookup in the getattro implementation gives priority to
section(x).mechanism.

[neuroman314]

But at least it is possible to write model cell specifications de novo, in NMODL and hoc|Python, in such a way that one can sample the states from a single cell at intervals during a cycle, then use this data to initialize any cell to any of the captured times. Blind reuse of existing model cells, no.

One potential problem with inducing phase shifts by perturbing cell trajectories is that, even if perturbation is effective, many cycles may be needed for the cell to settle. This could be a special problem for biophysical models that involve ion accumulation.

[Helveg]

Good point neuroman, I could as a default strategy keep running the model until the period becomes stable before I measure the states.