How to use hoc to plot a variable vs. time

Particularly useful chunks of hoc and/or NMODL code. May be pedestrian or stunningly brilliant, may make you gasp or bring tears to your eyes, but always makes you think "I wish I had written that; I'm sure going to steal it."
Post Reply
ted
Site Admin
Posts: 6286
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

How to use hoc to plot a variable vs. time

Post by ted »

The problem: how to plot a variable vs. time with hoc code

This problem seems to come up all the time: you want to write some hoc code that will plot certain variables as functions of time. You want to plot these in graphs that are automatically updated in the course of a simulation. Here and there on the WWW, you've seen a bunch of "legacy code" for plotting simulation results, and it's full of odd statements that you don't quite understand. So what do you do?

The easiest way to proceed involves stealing code from a session file. If you use the GUI to create a graph of a variable, then save that graph to a ses file, the ses file will contain hoc statements that specify not only the variable to be plotted, but also axis scaling, and--no less important--properly attach the graph to NEURON's standard run system so that automatic updates occur during a simulation. No more agonizing about beginline, flush, fastflush, doNotify, or any of that stuff.

Stealing code from a session file

For this example, suppose you need to plot membrane potentials, ionic current densities, and gating states.

Make a single compartment model and insert a voltage-gated density mechanism into it, e.g. hh.
Using the NEURONMainMenu toolbar, create:
--a Voltage axis graph of v(0.5) vs. t
--a Current axis graph of soma.ik(0.5) vs. t
--and a State axis graph of soma.m_hh(0.5) vs. t.
Also, use each Graph's "Change text" to write an appropriate label on each graph ("voltage axis", "current axis", or "state axis").
Then use the Print & File Window Manager (PFWM) to select those three graphs and save them to a session file (read how in the FAQ list).
Finally, open that ses file with a text editor.

The ses file will look something like this:

Code: Select all

objectvar save_window_, rvp_
objectvar scene_vector_[8]
objectvar ocbox_, ocbox_list_, scene_, scene_list_
{ocbox_list_ = new List()  scene_list_ = new List()}
{
save_window_ = new Graph(0)
save_window_.size(0,5,-80,40)
scene_vector_[5] = save_window_
{save_window_.view(0, -80, 5, 120, 1071, 66, 300.48, 200.32)}
graphList[0].append(save_window_)
save_window_.save_name("graphList[0].")
save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)
save_window_.label(0.629393, 0.392971, "voltage axis", 2, 1, 0, 0, 1)
}
{
save_window_ = new Graph(0)
save_window_.size(0,5,-1,1)
scene_vector_[6] = save_window_
{save_window_.view(0, -1, 5, 2, 1069, 333, 300.48, 200.32)}
graphList[1].append(save_window_)
save_window_.save_name("graphList[1].")
save_window_.addvar("soma.ik( 0.5 )", 1, 1, 0.8, 0.9, 2)
save_window_.label(0.261981, 0.819489, "current axis", 2, 1, 0, 0, 1)
}
{
save_window_ = new Graph(0)
save_window_.size(0,5,0,1)
scene_vector_[7] = save_window_
{save_window_.view(0, 0, 5, 1, 1064, 599, 300.48, 200.32)}
graphList[2].append(save_window_)
save_window_.save_name("graphList[2].")
save_window_.addvar("soma.m_hh( 0.5 )", 1, 1, 0.8, 0.9, 2)
save_window_.label(0.332268, 0.479233, "state axis", 2, 1, 0, 0, 1)
}
objectvar scene_vector_[1]
{doNotify()}
Let's analyze the code for the voltage axis graph to discover its underlying pattern,
and figure out how you might reuse it for your own purposes.

This

Code: Select all

{
save_window_ = new Graph(0)  // creates but does not display a new Graph
  // so you can control its size and scaling before it is even drawn to the screen
save_window_.size(0,5,-80,40)  // axis scaling
scene_vector_[5] = save_window_
{save_window_.view(0, -80, 5, 120, 1071, 66, 300.48, 200.32)}  // draws it on the screen
  // in a window with user-specified location (5th and 6th args) and size (last 2 args)
graphList[0].append(save_window_)  // graphList[0] is for all graphs whose y values 
  // are to be plotted vs. t
save_window_.save_name("graphList[0].")
save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)  // adds v(0.5) to this graph's plot list
  // also specifies line color and style, location where "v(.5)" will be printed
save_window_.label(0.629393, 0.392971, "voltage axis", 2, 1, 0, 0, 1)  // just prints the 
  // string "voltage axis" on the canvas at specified location, color etc.)
}
is the code that sets up the voltage axis graph, with explanatory comments by me. This is actually a "reusable code pattern" that you can copy, paste, and modify to suit your own needs.
  • Be sure to read the Graph class's documentation for a more detailed description of the various graph methods that are being called. Also, I'm doing this from memory and may have made a mistake.
Here's how you would use this code pattern to create a graph that plots membrane potential at the 0.1 location of apical[3] and the (0.4) location of basilar[21]. You want the x axis to run from 0 to tstop, and the y axis to run from -10 to +5.

Code: Select all

objref vg
vg = new Graph(0)  // creates but does not display a new Graph
vg.size(0,tstop,-10,5)  // axis scaling
// forget about the scene_vector stuff
vg.view(0, -10, tstop, 15, 200, 200, 300.48, 200.32)  // draws it on the screen
  // in a window with user-specified location (5th and 6th args) and size (last 2 args)
graphList[0].append(vg)  // graphList[0] is for all objects that are to be 
  // updated at integer multiples of dt
// forget about the save_name stuff
vg.addexpr("apical[3].v(.1)", 1, 1, 0.8, 0.9, 2)
vg.addexpr("basilar[21].v(.4)", 1, 1, 0.8, 0.7, 2)  // put variable name a bit lower on the canvas
vg.label(0.4, 0.9, "two at once!", 2, 1, 0, 0, 1)
This will put the graph on your screen at (200, 200). If you don't like where it is, or want to change the position of the labels and/or color/linestyle of the traces, just drag it to a better location and use the graph's own menu to make whatever changes you like.
Then use the PFWM to save it to a session file and steal the new coordinates, line style args, etc. from that file.

Now you're ready to analyze the original ses file and extract the code patterns for the current axis and state axis graphs. graphList[1] is for graphs whose y values are to be plotted vs. t - 0.5*dt, and graphList[2] is for graphs whose y values are to be plotted vs. t + 0.5*dt. The reason for these strange time offsets is that, if you are using NEURON's "staggered Crank-Nicholson method" (secondorder == 2), ionic currents and states will be second order correct at t - 0.5*dt and t + 0.5*dt, respectively (membrane potential will be second order correct at t). If you are using the adaptive integrators (cvode), all variables are computed at the same t, and plots generated with graphList[1] and graphList[2] are identical to plots of the same variables that are generated with graphList[0].

An exercise for the reader

Write an obfunc that spawns a new graph that plots a user-specified variable vs. t. It should accept these arguments:
--xy location on screen
--a strdef that is the name of the variable to be plotted
--y axis scaling (i.e. ymin and ymax)

Acknowledgments

"This programming pearl has been brought to you by NEURON's standard run system."

"Aren't you glad you use NEURON's standard run system?
Don't you wish everyone did?"
Post Reply