Nsingle

[porio]

Hi,

I'm trying to model stochastic single channels instead of continuous conductances. I found the singhh.hoc demo, but I couldn't figure out how to use the Nsingle function with my mechanisms created with NMODL. I guess that I have to point to them with some object name, isn't it?
Also, what does Nsingle do, exactly? What is the reason of the nsing() procedure that is found inside the demo?

Regards

[ted]

Don't use NMODL. In NEURON, stochastic channel gating mechanisms are implemented
with the KSChan class, which you can read about here
http://www.neuron.yale.edu/neuron/stati ... schan.html
You will note that singhh.hoc starts with
load_file("singhhchan.hoc")
The statements in singhhchan.hoc specify a kinetic scheme representations of the
HH gna and gk. If you read that file, you'll see a lot of stuff that may look rather baffling,
but don't worry about it--use the Channel Builder, which is the GUI interface for
managing the KSChan class.

Here's what you need to do for each channel that you want to simulate with stochastic
gating:
1. Devise a kinetic scheme representation. For any HH-style set of differential equations
there is a corresponding kinetic scheme.
2. Start with
NEURON Main Menu / Build / Channel Builder / Point
and then use the channel builder to specify the channel's properties. Be sure to click on
Properties / Allow Single Channels
Then set up channel name, ionic selectivity, ohmic vs. GHK, etc..
Before starting, you might want to work through the
Channel Builder tutorial
http://www.neuron.yale.edu/neuron/stati ... /main.html
if you haven't already done so. The second part of the tutorial involves setting up
a three-state kinetic scheme potassium channel. You might want to try turning it into
a stochastic channel.
3. Once you have a stochastic channel, you can use it like any other Point Process.
It will appear in the Point Process Manager's list of mechanisms, and Nsingle will be the
first parameter in that tool's parameter panel. If you set Nsingle to 0, the channel
conductance will be a continuous function of time. If Nsingle is > 0, you'll see discrete
gating of that number of channels, and gmax specifies the open conductance of a
single channel. With adaptive integration, each channel opening or closing transition
is handled as an independent event. With fixed dt, all transitions are forced to occur
at time step boundaries.

--Ted

[porio]

Thanks for your help, Ted. I followed your instructions and everything worked.
However, there is one thing that I need to do and I don't know how to do it with ChannelBuilder (and that's why I went with NMODL at the beginning): for some channels (if not all) I need to account for the changes in the kinetic constants due to temperature changes. The stochastic mechanism I have just built with ChannenBuilder works exactly as the NMODL mechanism... but only at 22ºC. I need the kinetic constants to change accordingly with temperature.
any suggestions?

Thanks.

[ted]

Export the channel specification to a hoc file
ChannelBuilder / Properties / Clone channel type / Hoc file for KSChan
Then modify the hoc file as needed to assert temerature dependent rate constants.
Refer to the Programmer's Reference entries on KSChan as necessary for details.

Example based on tutorial at
http://www.neuron.yale.edu/neuron/stati ... tline.html
That tutorial constructs a distributed mechanism, so start by
ChannelBuilder / Properties / Clone channel type / As Point Process
and save the output to mykhhpp.ses, then exit NEURON.
Restart NEURON and load mykhhpp.ses.
Convert this to single channel form so it can be used to simulate stochastic channel
gating
ChannelBuilder / Properties / Allow Single Channels
Now export to a hoc file
ChannelBuilder / Properties / Clone channel type / Hoc file for KSChan
and save the output to mykhhpp.hoc
Exit NEURON.
Open mykhhpp.hoc with a text editor.
The statements

Code: Select all

{
  tobj = kstransitions.object(0)
  tobj.type(1)
  tobj.set_f(0, 4, ksvec.c.append(1, -0.164, -48.8))
  tobj.set_f(1, 2, ksvec.c.append(4.4, -0.025, -65))
}
{
  tobj = kstransitions.object(1)
  tobj.type(1)
  tobj.set_f(0, 4, ksvec.c.append(1, -0.036, -22))
  tobj.set_f(1, 2, ksvec.c.append(2.6, -0.007, -65))
}

clearly specify the properties of the two transitions, and the 2nd arguments to
ksvec.c.append() are the rate constants.

Replace each rate constant with a symbol, e.g. kf1, kb1, kf2, kb2.

Write a proc setrates() that specifies how these are affected by the variable
celsius
and a custom proc init() that calls setrates()

Code: Select all

proc init() {
  setrates()
  finitialize(v_init)
}

and put both of these procs in a file called setrates.hoc

Change your init.hoc to read
load_file("nrngui.hoc")
load_file("mykhhpp.hoc")
load_file("setrates.hoc")
. . . other statements that load files that specify the model,
instrumentation, and user interface . . .

Now double clicking on init.hoc (or executing nrngui init.hoc) will pull all the bits together.

--Ted

[porio]

Great, sounds simple and fun to do.
but before I start doing the trials, let's say that I change the temperature while the simulation runs (as I'm doing, indeed). Will this work for me, updating the constants accordingly?

[porio]

Hi Again, Ted.
I'm afraid it's not working. How can I attach the files to a message so that you can see that I'm really doing as you told?

Regards.

[ted]

Thanks for a tough question, porio. Working out the answer taught me a lot. First of all,
I discovered that my initial suggestion wasn't going to work. The hoc file that is emitted by
the Channel Builder uses a bunch of objrefs to set up the internal data structures that
NEURON uses to simulate the channel. However, the very last statement in that hoc file
is {objref ks, ksvec, ksgate, ksstates, kstransitions, tobj}, which destroys the links
between these objrefs and the internal data structures. Result: you have a new object
class that you can treat like any other point process class, but hoc can only see the
public members of this class--which are Nsingle, g, gmax, i, e, and the gating states
themselves. There is no obvious way to get at the parameters for the equations that
govern the state transitions.

A little digging reveals that the parameters are actually stored as elements of Vectors,
so in theory one could access them from hoc using statements of the form
print Vector.x[j]
Vector.x[j] = somevalue
But this means using the "unique name" of an object, which is a bad idea. Quoting from
chapter 13 of The NEURON Book:

Object references vs. object names

Up to this point we have been using object references to refer to objects, emphasizing
the difference between an object itself and what we call it. Actually, each object does
have a unique name that can be used anywhere a reference to the object is used.
However, these unique names are primarily intended for use by the library routines that
construct NEURON's graphical interface. While it may occasionally be useful to employ
these unique names in user-written code (e.g. for diagnostic or didactic purposes),
this should never be done in ordinary programming. Object names are not guaranteed
to be the same between different sessions of NEURON unless the sequence of
creation and destruction of objects of the same type is identical. This is because the
object name is defined as classname[index], where the "index" is automatically
incremented every time a new instance of that class is created. Index numbers are
not reused after objects are deleted except when there are no existing objects of
that type; then the index starts over again at 0.

So to change rate parameters from hoc, we need to use aliases. To make aliases
available, the KSChan must be managed by a Channel Builder. Here's how:

Code: Select all

(this isn't really code; I'm just using this "code block" to preserve the indented 
formatting of the following outline)

Get the channel specification into a Channel Builder
  nrngui sthcn.hoc
  NEURON Main Menu / Build / Channel Builder / import KSChan / sthcn

Aliases exist only when requested by the user, so
  ChannelBuild[0] / Properties / Provide transition aliases
For future use, save this Channel Builder to a session file 
with aliases enabled
  File / save session / sthcn_with_aliases.ses

Use a graph's Plot what? to discover the alias names
  NEURON Main Menu / Graph / Current axis (this does not require a default section)
  Graph's menu box / Plot what?
  Plot what? / Show / Objects
  scroll down list in first column to ChannelBuild_
  double click on ChannelBuild_ -->
    ChannelBuild_ appears in text edit field at top of Plot what?,
    and 0. appears in Plot what?'s middle column, 
  double click on 0. in middle column -->
    edit field changes to ChannelBuild[0].,
    and middle column shows a list of names
  double click on aliases. in middle column -->
    edit field changes to ChannelBuild[0].aliases., 
    and right column shows aC202., aCC2., etc.
  double click on aC202. in right column -->
    edit field changes to ChannelBuild[0].aliases.aC202.,
    and right column shows A, d, k

So the "complete" names of these three parameters are
ChannelBuild[0].aliases.aC202.A, ChannelBuild[0].aliases.aC202.d, etc.

-------- an aside --------
Well, we don't have to refer to Vector (good because we and the GUI are liable to
create lots of Vectors, so we especially want to avoid referring to Vectors by their
unique names). We still have to refer to a particular ChannelBuild object, but at least
creation and destruction of ChannelBuild objects is something that is completely under
our control, and isn't going to be done automatically by the GUI or runtime library.

If we had a lot of point process ChannelBuild objects, and especially if they controlled
different kinds of channels, we'd probably want to use Lists to manage them (i.e. for
each kind of channel, there would be a List whose elements were each instance of
that channel). That would make it easy to iterate over all instances of a particular
kind of channel, in order to specify their properties.
-------- end of aside --------

Testing what we think we know about aliases: type one of them at the oc> prompt
and see if hoc returns its value--

Code: Select all

oc>ChannelBuild[0].aliases.aC2O2.A
        0.0067

So it works.

For your particular channel model, the forward rate parameter aliases are
ChannelBuild[0].aliases.aCO.A
ChannelBuild[0].aliases.aC2O2.A
ChannelBuild[0].aliases.aCC2.A
ChannelBuild[0].aliases.aOO2.A
and the backward rate parameter aliases are
ChannelBuild[0].aliases.bCO.A
ChannelBuild[0].aliases.bC2O2.A
ChannelBuild[0].aliases.bCC2.A
ChannelBuild[0].aliases.bOO2.A
and we have enough information to set them all from hoc.
Now we need a custom proc init() that will assign the proper values to these, depending
on celsius. A first stab at this, using the temperature dependence that you are
assuming, is

Code: Select all

proc setrates() {
	qv = 4^((celsius - 22)/10)
	ChannelBuild[0].aliases.aCO.A = 0.0015*qv
          . . . plus statements to assign values to the 7 other rate terms . . .
}

proc init() {
	setrates()
	finitialize(v_init)
}

[porio]

Thanks a lot for your answer Ted!!
I haven't had the time to work on it, but as soon as I have I will tel you know.
Regards.

[porio]

Hi Ted
I didn't have the time to work on this until yesterday. And it works!... almost.
You see, the only thing I miss (and I really need it) is that the kinetic constants are updated during the simulation. This is because I'm working with temperature ramps.
The thing would be very straight-through if we had either:
a) the ability to enter user-defined expressions for the transitions in Channel Builder, or
b) a way of specify stochastic single channels in NMODL. And this can not be very difficult, is it? I have some thoughts about it but I would like to have some hints.

[hines]

You can change the vector values at any time, even every time step, during a run.
However, in the single channel (Nsingle > 0)
case there is an internal data structure based
on the rate values that you will need to
tickle to get it updated if any of the rate
vector values change. These normally
get updated whenever a ligand concentration
is used or the voltage changes from the previous step by more than a specified resolution, KSChan.vres(value) ( don't remember if this is settable from the ChannelBuilder) Anyway, probably the easiest thing to do is add a ligand sensitive transition between any two states and set
the rate values for it equal to 0. This will force an update of the internal structure every time step and your continuous rate changes will take effect. This is not necessary for continuous channels except that you do not want to select the use of tables for the fixed step method. Anyway, I'm curious if you see any kind of performance hit with that otherwise useless ligand sensitive transition.

[porio]

Mmmmm... I'm not sure if I'm getting it. In my scheme there is a ligand-dependent transition (it's a for states schemes, with two v-dependent and two ligand-dependent transitions). So what you're saying whould be happening already, isn't it? Well, it is not :(

[hines]

That's right. Since you have a ligand sensitive transition, the internal data structure is updated every time step. Since you are still having problems, you should send me (michael.hines@yale.edu) the hoc/ses/mod files in a zip file that are needed to run the model.

[porio]

Thanks for the answer. To make it available for everyone, the only thing that was missing is an advance procedure, very similar to the proc init(). So the setrates.hoc now looks like this:

Code: Select all

proc setrates() { 
   qv = 4^((celsius - 22)/10) 
   ChannelBuild[0].aliases.aCO.A = 0.0015*qv 
          . . . plus statements to assign values to the 7 other rate terms . . . 
} 

proc init() { 
   setrates() 
   finitialize(v_init) 
} 

proc advance() {
   setrates()
   fadvance()
}

Now the proc setrates() is called on every time step.
However, I wonder what is the impact of this on the calculation speed. Is this similar to the NMODL mechanisms where the rates are usually called on every timestep? (unless there is a table statement). Will this be slower than that?

[ted]

I'm glad you posted that code.

These are good questions:

Now the proc setrates() is called on every time step.
However, I wonder what is the impact of this on the calculation speed. Is this similar to the NMODL mechanisms where the rates are usually called on every timestep?

Similar in effect.

(unless there is a table statement)

If temperature is really changing throughout a simulation, TABLE serves
no useful purpose. Worse than not having a TABLE statement, in fact, since
the whole table has to be recomputed at each time step.

If temperature changes only at specific times t0, t1, . . . , an
FInitializeHandler and the CVode class's event() method can be used
to call a proc that changes the rates at just those times. The use of
FInitializeHandler and CVode's event() method has been discussed
elsewhere in this Forum in connection with making other kinds of
parameter changes.

General comment regarding "speed"--
When considering how different programming strategies may affect
simulation speed, remember also to ask if the effect is meaningful.
Programmer time is always more valuable than machine time. Some
mechanisms are quite difficult to specify with NMODL, e.g. stochastic
channel gating.

[hines]

...I wonder what is the impact of this on the calculation speed...

The question may be moot since nmodl is not capable (without VERBATIM block c programming) of defining single channel models. However, from the perspective of deterministic channels, KSChan should be slightly (though probably not measurably) faster than NMODL. This is because they are computing with the same algorithms but NMODL channels compute the jacobian contributions di/dv using two calls to the code in a BREAKPOINT block, di/dv = (i(v+.001) - i(v))/.001whereas KSChan computes it analytically, di/dv = g(v). The major performance difference will be due to calling the "setrates()" procedure every time step which is a hoc procedure and therefore, as a rule of thumb, 10 to 100 times slower than equivalent c code. Anyway, one can easily determine performance differences with the idiom,

Code: Select all

runtime = startsw()
run()
runtime = startsw() - runtime
print "runtime = ", runtime

Turning attention to the performance of your single channel simulations it will help to understand a bit of what is going on during the simulation. The example you sent me had Nsingle = 3000, 4 states and 8 transitions. Since at least one of your transitions is ligand dependent, the rates are considered to change every time step and this forces the re-computation of an internal table of length 8 where the values are equal to the 8 transition rates (for specific time step voltage and ligand concentrations). During the time step we compute a sequence of random number pairs and a sequence of new transition probability tables normalized to 1 where the interval between the elements depends on the integer population of each state and the transition rate away from the state. (that is definitely the performance limiting step). The first of the random pair is drawn from a negexp distribution and defines the time of the next transition, ie. one of the 3000 particles changes its state. The second of the pair is drawn from a uniform distribution 0-1 and defines (from the population dependent table) which transition is used and thus the source state population of that transition is reduced by one and the target state population is increased by one. The random pairs and population dependent transition probability table keep getting recomputed until the next transition time is greater than t+dt. Note that this is a very fastidious and statistically precise way to compute stochastic channels. Any performance improvements would have to be analyzed mathematically to verify that the improvement did not come at the cost of too much statistical crappiness.
Note that ONE single channel, Nsingle=1, is computed by quite a different algorithm than when Nsingle > 1.