% % ls ( directory listing ) % cp file1 file2 ( copy file1 to make file2 ) % logout ( logs you off the computer )
For example:
% vi circuit1After you enter the code for a NeuronC script, run the script by entering one of the following commands:
neurc file.n Displays graph on screen.
nc file.n Displays only numbers on output.
Used for "debugging" a model.
nc file.n > file.r Saves numbers in file.
Used for running long simulations.
nc -v file.n | vid (equivalent to "neurc file.n")
plotmod file.r | vid Displays graph on screen from file.
plotmod file.r | 8by11 Draws graph on color plotter.
To build a neural circuit, you need to include at least one neural element statement and some stimulus and plot or print statements:
conn 1 to 2 cable length 100 dia 1 rm 5000;
stim node 1 vclamp -.01 start .05 dur .01;
plot 1,2;
step .1;
This simple NeuronC program builds a cable and
defines a stimulus of a voltage clamp to minus 10 mv,
starting at 50 msec into the experiment and stopping
after 60 msec. The voltage at each end of the cable is
plotted. The plotting stops at 0.1 seconds.
When setting up a neural circuit model, remember that a simple model is easier to understand and possibly more believable than a complex one. This is generally true because the utility of a computational model is often that it quantitatively assists (stretches) our imagination. If a model is so complex that we don't understand it, then it's not helping us. If a model uses many parameters for which numerical values are not known, then it may be a focus for disbelief.
The Euler integration method (euler = 1) normally requires a time step on the order of 1 usec (1e-6 sec) so it is usually slower than the Crank-Nicolson method. It is instructive, however, to experiment using different Euler time steps to see under what conditions instability occurs.
After the anatomy and basic synaptic model is all working properly, you should expand the model to include more cells and complexity. The advantage of this technique is that any problems that occur can usually be localized to the addition of a few neural circuit elements, because everything was developed and verified in successive stages.
Several features of NeuronC facilitate the use of "developmental rules". Node numbers (used to define connections between neural elements) may have up to 3 dimensions, and these dimensions are appropriate to define: (dimension 1) neuron types, (dimension 2) neuron number (dimension 2) and (dimension 3) neural morphology. Three dimensional locations and the "ndist" function are useful for finding the distance between two nodes. Finally, it is possible to define stimulus and recording nodes on the basis of cell type and location, by creating a "smart" function that scans the neural circuit and automatically finds the appropriate node.
-----------------------start of morphology file---------------------
#
# beta cell morphology file.
#
# node parent dia xbio ybio zbio region dendr
#
0 0 21 0 0 -25 SOMA 0
1 0 2.17 21 -2.47 -8 PROX 1
2 1 0.833 22.8 0.0287 -5 DEND 1
3 2 1.33 23.3 0.776 0 DEND 1
4 3 0.5 23.9 6.32 0 DEND 1
5 0 1.5 -3.25 -33.6 -30 HILLOCK 2
6 5 0.5 -6.81 -37.6 -30 AXON_THIN 2
7 6 0.5 -25.5 -66.4 -30 AXON_THIN 2
8 7 0.5 -121.2 -160.2 -30 AXON 2
and so on.
--------------------------end of morphology file----------------------
Each point where 2 or more cables connect is a "node", which
is defined on one line. The line contains the node's position
(xbio,ybio,zbio), its diameter, and its parent node, as well a
label defining what region of the cell to assign it to (for
membrane channels, colors, etc.). The "dendr" column is a dendrite
number useful for analyzing the properties of each separate dendrite.
The list of points is read into an array by the "fread" statement
in the main script file. Comments start with "#". You could readily
modify or add to this format to better suit your needs. Note that
for sequentially-labeled nodes along a dendrite, each node's
parent is the previous node. At a branch point, two or more
nodes have the same parent.
Any value in the morphology file can also be specified by a parameter, which must be defined before the file is read by the simulator. For example:
# node parent dia xbio ybio zbio region dendr
#
0 0 21 0 0 -25 SOMA 0
1 0 d_dia 21 -2.47 -8 PROX 1
2 1 d_dia 22.8 0.0287 -5 DEND 1
Where "d_dia" is the dendritic diameter, defined in the simulation
script before the morphology file is read. Then this parameter can
be set automatically according to the specific simulation.
An advantage of this format is that it can represent any type of
morphology, including tree-structured neurons and also syncytial
networks of neurons connected by gap junctions.
The node numbers should be fairly sparse but dense enough that they capture the essential details of the dendritic tree. You may need to experiment a bit with the number and density of nodes to find the appropriate level of detail you want to preserve. For each node, you label its position with a "text" number with depth 50 (i.e. closer than the cell image). The number labels are arbitrary, but it helps to start by labeling the soma with a zero, then label the first dendrite with sequentially increasing numbers. Each branch starts non-sequentially. After all the points are entered, you can save the information as a .fig file, and convert it to the Neuron-C format using "fig2nc".
As a second step, you must edit the file using a text editor to correct the parent node numbers. Normally for a sequentially labeled dendrite each node's parent starts out correctly set as the previous node number. However you must correct the parent node for the the first node of each new branch. You must also add the dendritic diameters for each node. This can be done by measuring with a ruler (making sure to convert to the correct scale), or by adding variable(s) in the file instead of numbers (as described above for "d_dia"), then define the variable in the simulation script. When the file is read into the simulation, the previously defined value is substituted. You can find the "fig2nc" script in "nc/models/retsim/runconf".
-----------------------start of fig2nc script---------------------
#
# fig2nc
#
# Script to convert .fig file data to standard anatomy format for Neuron-C
# when a photograph of a cell is labeled with "text" numbers to represent
# the node numbers in nc.
#
# The parent node is assumed here to be the previous one, which is not
# always correct but since many nodes follow an unbranched dendrite it
# is best to start. Parent nodes then need to be corrected by hand.
#
# original scale is 19200 per 100 um (scale=192), more recently reduced to scale=40
#
awk 'BEGIN {scale = 40;
par = 0;
dia = "amdia"
reg = "DEND";
dend = 0;
z = 0;
printf ("# node parent dia xbio ybio zbio region dendr\n#\n");
}
NR>10 && $1=="4" {
x = $12/scale; y = $13/scale; n=substr($14,1,length($14)-4);
if (n=="0") { xs=x; ys=y; }
printf (" %4s %4s %4s %-8.4g %-8.4g %g %s %d\n",
n, par, dia, x-xs, y-ys, z, reg, dend);
par = n;
}
' $1
-----------------------end of fig2nc script---------------------
http://rsb.info.nih.gov/ij/download.html
The MorphDig plugin can be built and run using a java tool called "ant", which is a Java-based build tool (like "make" , but portable with Java). Ant can be downloaded from:
IJ-MorphDig can be used with the Retsim retinal simulation package included in the Neuron-C distribution. MorphDig allows the user to trace a cell's morphology from a stack of confocal images by positioning points on the images. Each point consists of an (x,y,z) location, the dendritic diameter, measured by a "caliper tool" set by the user. The data is output in the form of Neuron-C morphology files (see above). IJ-MorphDig was written by Mike Schachter (mike.schachter@gmail.org).
--------------------start of density file---------------------------
densdata = {
/* Densities for the cell's different regions (n/um2). */
/*
dend prox soma hillock thin axon */
50, 60, 80, 100, 100, 50, /* Na */
2, 2, 3, 0, 0, 0, /* Ca */
0, 0, 0, 0, 0, 30, /* Kdr */
3, 5, 6, 0, 0.1, 0.1, /* KCa */
30, 30, 30, 0, 0.1, 0.1 /* Ka */
};
---------------------end of density file----------------------------
Example of a simulation script to read these morphology and density files:
---------------------start of definition file----------------
/* gcdefs.n */
/* Definition of columns for "anatfile" which defines cell anatomy:*/
NODE = 0;
PAR = 1;
DIA = 2;
AX = 3;
AY = 4;
AZ = 5;
REG = 6;
/* Regions: Column numbers */
DEND = 0;
PROX = 1;
SOMA = 2;
HILLOCK = 3;
AXON_THIN = 4;
AXON = 5;
nregions = 6;
/* Channel types: Row numbers */
NA = 0; /* noise type 1 */
CA = 1;
KDR = 2; /* type 0 , noise: type 1 */
KCA = 3; /* type 2 , noise: type 6 */
KA = 4; /* type 3 , noise: type 5 */
nchans = 5;
dim chname[] = { "Na", "Ca", "Kdr", "KCa", "KA" };
dim regname[] = { "dend", "prox", "soma", "hillock", "thin",
"axon" };
dim densdata [nchans][nregions];
proc printdens (label) /* procedure to print out densities */
{
printf ("# \n");
printf ("# Channel densities (%s)\n",label);
printf ("# \n");
printf ("# ");
for (ch=0; ch<nregions; ch++)
printf ("%-8s ",regname[ch]);
printf ("\n");
printf ("# \n");
for (ch=0; ch<nchans; ch++) {
printf ("# ");
for (region=0; region<nregions; region ++) {
printf ("%-8.4g ",densdata[ch][region]);
};
printf (" %s\n",chname[ch]);
};
printf ("# \n");
};
-------------------------end of definition file--------------------------
This scheme allows you to read in different densities by defining
a different density file to replace "stddens.n", and different
anatomies by defining a different anatomy file.
-------------------------start of main script file-----------------------
/* main script file */
/* main.n */
timinc = 1e-4; /* simulation variables */
ploti = 1e-3; /* these are defined by the simulator */
endexp = 1.0;
expt = 2; /* which experiment to run, */
/* defined only in this script */
x = setvar(); /* set values from command line */
/* You can change "expt", or "anatfile", etc. */
/* like this:
/* nc --expt 3 --anatfile alpha2 main.n */
if (notinit(deffile)) deffile = "gcdefs.n";
include (deffile); /* definition file */
if (notinit(densfile)) densfile = "stddens.n"
include (densfile); /* define channel densities for cell regions */
if (notinit(anatfile)) anatfile = "betagc";
fread (anatfile, cabldata, clines, ccols); /* read anatomy file */
/* "fread" sets clines = number of lines in file */
/* ccols = number of lines in file */
printdens("N/um2"); /* print out densities actually used */
na_type = 0;
proc makseg(nod1, nod2, d, region) /* procedure to make soma/dendrites */
{
if (region == SOMA){
at [gang][nod1]
sphere dia d vrest=vstart /* this is just the soma */
chan Na type na_type
ndensity=densdata[NA][region]
else conn [gang][nod1] to [gang][nod2] /* else make dendritic segment */
cable dia d vrest=vstart;
chan Na type na_type
ndensity=densdata[NA][region] /* add more chans if wanted */
};
/* Now, make the cell */
for (i=0; i<clines; i++) { /* locate the nodes */
n = cabldata[i][0];
at [gang][n]
loc(cabldata[i][AX],
cabldata[i][AY],
cabldata[i][AZ]);
};
for (i=0; i<clines; i++) { /* make the dendritic tree */
makseg (cabldata[i][NODE],
cabldata[i][PAR],
cabldata[i][DIA],
cabldata[i][REG]);
};
/* Code to make stimuli, set up electrode recordings,
and run the experiment. */
if (notinit(expt)) expt = 1;
if (expt==1) {
"code to test simple stimulus for calibration"
}
else if (expt==2) {
"code to run standard stimulus"
}
--------------end of main script file ------------------------------
The idea is to have standard values and definitions that have
default values but which the user can modify at will.
There are several ways to deal with edge effects. One way, the simplest, is to "live with" edge effects and try to analyze the simulation data in a way that minimizes their influence. Another way is to construct a larger model than necessary, so that data can be taken exclusively from a central region without edge effects. A third way is to "wrap around" one edge's connections to the opposite edge of the array. The resulting model then simulates a large array with a repeating structure, such that a single stimulus "looks like" a stimulus that is spatially repeated over the biological tissue.
A fourth way to deal with edge effects is to connect to each of the edge neurons a small circuit that simulates the effect of an infinite network. If it is not possible to do this exactly, an approximation is better than nothing. In the case of an interconnected array of passive neurons, the effect of an infinite network can be simulated by an impedance equivalent to the impedance "seen" by one neuron "looking" laterally at its neigbor in one direction. This can be accomplished in several ways, depending on the type of neural interactions.
An easy way to simulate an impedance in a steady-state network is to ignore any capacitance and simply attach a load resistor to the edge neurons. The exact value for the load resistor can be calculated by an interative procedure where the input resistance of a centrally located cell is measured from its neighbor through their common connection (e.g a gap junction). One then connects this resistance to the edge cells, and then repeats the measurement at the central cells. After a few iterations of this procedure, the value of resistance settles down and then represents accurately the resistance of an infinite network. Then a stimulus applied in the center of the array will decay towards the edge through the interconnections in a way that accurately simulates an infinite network.
However, if each cell in a network has an identical but nonzero resting voltage, the passive resistor scheme outlined above won't work because the resistor shunts the signal from edge cells. The solution in this case is to connect a "battery" identical to the resting voltage to the appropriate passive resistor. This can be accomplished conveniently by a "buffer amplifier" or "voltage follower" that causes the voltage in one cell to follow the voltage somewhere else in the array. In this case, one connects a voltage follower from each cell one row inside the edge to its neighboring edge cell. Then the edge cells follow their central neighbors and no edge effect is apparent.
In a more complex neural circuit with many neuron types of various sizes and interconnections, it may be very difficult to eliminate edge effects with the "impedance" method outlined above. There are few rows of larger cells in the array, and their edge effect perturbs more rows of smaller cells. The solution is to use a network large enough to include one extra row of the largest cells at the edge, and then use the voltage buffer idea to make these cells "mirror" their inside neighbors.
statements to create first network... run; statements to save voltages from first network... statements to create second network... statements to reset "time" and "endexp" run;Several other housekeeping chores must also be done to get this scheme to work. After the first run, the first network is left intact, so voltages may be read from any node. Before the second network is run, the variables "time" and "endexp" must be reset to their proper values. The "time" variable represents what value of absolute time to start the next run. The "endexp" value represents the absolute time to stop the run. It is possible to set "time" negative, so that zero time corresponds to a point delayed after the start of the simulation. This is useful, for instance, when a network is run to give some time for equilibrium and then a stimulus is given.
Explanation:
1) The easiest way to make a new channel type is to define the closest type, and then shift the rate functions, either by offsetting the voltage activation curves or by a changing the rate function multiplier. You can change voltage offsets with the "offset, offsetm, offseth" parameters, and the multipliers with the "taua, taub, tauc, ... tauf" parameters. The advantages of this method is that it is simple and is done separately for each channel. This allows you to define several "channel types" each with different but closely related rate functions.
The "tau" parameters define a "time constant" (1/rate) which multiplies the default rate function to make it faster or slower. For HH channels, the "taua" parameter is for the forward activation function (alpham), and "taub" is for the reverse function, "tauc" is for inactivation, and "taud" is reverse inactivation. The "tauf" parameter is for flicker rate (channel noise in the conducting state).
For Markov sequential-state channels, to find the which "tau" parameter to use, look at the source code for the channel. The "rateo" parameter defines which one (rateo=0 -> taua, rateo=1 ->taub, ... rateo=5 -> tauf).
2) A an easy way to change the behavior of all the channels of a channel type is to redefine the rate functions. This provides more precise control of the rate functions than changing only the voltage offset and the multipler (as above). To rewrite the rate functions, look at the existing default rate function in the source code for the channel, and copy it. For the interpreted version of the simulator, rename the new function in the script (".n") file to one of the predefined names "calcna1", "calcna2", "calck1", etc. (see the source code for the channel type to find the "user defined" rate function's name). When the "nc" interpreter sees a function defined with such a name in the script file, that function is called instead of the default one when calculating the rates for the channel. For the compiled version, you call "set_chancalc()" to define a different rate function. This way you can modify the functions during a run if necessary. One drawback of this method is that it modifies the rate function for all the channels of that type. Another drawback is that the rate function in the script file is run by the interpreter so it runs slower than if it were compiled. Therefore this method is best used for rate functions that are precalculated and stored in lookup tables [i.e. those (such as HH functions) that are functions of only one parameter such as voltage].
3) A slightly more involved method is to copy an existing channel type to make a new channel type, and then modify the new type to your requirements. You can change the number of states, their conductances, transitions to other states, their rate functions and multipliers. This is the easiest way to change a channel type's exact behavior:
copychan (source name, source type, new name, new type, number_of_states)
setchan_ntrans (name, type, state, number_of_transitions)
setchan_cond (name, type, state, conductance)
setchan_trans (name, type, state, transition, trans_state)
setchan_trate (name, type, state, transition, rate_function)
setchan_mul (name, type, state, transition, rate_multiplier)
setchan_rateo (name, type, state, transition, rate_mul_number)
These functions (in "chanfunc.cc") allow you to copy a channel type and to set all of the channel type's state information that defines the channel's behavior. You can read more about how to do this below or you can look to see how it's done in the source code. To change the number of states, you can set the number_of_states to be greater than the old channel type. Then you must add all the parameters for the new states. If you set "number_of_states" to 0, the number of states for the new channel is set the same as the old.
4) Finally, you can add a new channel type in the source code. To understand how to do this, the easiest way is to copy the source code for an existing channel type and modify it. Several source code files are relevant, "ncomp.cc", where compartment currents are calculated, "chanfunc.cc" (where function look-up tables are computed), "initchan.c", where some useful functions are defined and all the channel tables are created, and the channel definitionnn files, "channa1.cc", channa2.cc", "channa3.cc", "chank1.cc", etc.
Nc contains two classes of channels: the pure HH type and the "sequential-state" type. The HH type is traditional but the sequential-state type is more flexible and is being used more by physiologists in recent literature to describe new channels once they have been characterized. Both use first-order differential equations to describe the kinetics (i.e. the m,n, and h and also the "states" in the sequential-state model all have alpha and beta "rate constants" which are functions of something; usually just voltage, but can be calcium conc., etc.).
The HH (Type 0) class of channel runs faster because fewer state variables need to be modified, and there are no "state transitions" as there are with the "sequential-state" channels. However, the "sequential-state" channels are much more flexible.
struct chanconst { /* master struct for chan constants */
short int ctype; /* type of channel: NA, K, etc. */
short int stype; /* subtype of channel: 0=HH, 1=SS,etc*/
short int numstate; /* number of states used */
short int numparm; /* number of params,set to 2 for HH */
double qc; /* Q10 factor for conductance */
double qt; /* base temp for conductance */
double qcond; /* Q10 multiplier for conductance */
float vrev; /* default reversal potential */
double *perm; /* rel permeability to ions */
float unitary; /* default unitary conductance */
float trconc; /* default tr conc multiplier */
chanstate *state; /* state definitions */
chanparm *parm; /* rate params */
float *respamp; /* response amplitudes */
};
The chanconst structure is set at the beginning of a NeuronC run
and does not change during run time. It defines which tables
contain the rate functions, the states and a few other important
variables. Every channel has a pointer to its chanconst
structure that defines its type.The "ctype" field contains the type (NA, K, etc), and the "stype" field contains the subtype number. The "numstate" and "numparm" fields contain the number of states and "parameter sets", and the basetemp parameter is the temperature at which the rate functions and channel conductances are defined. You set these fields when you run the "makchanconst(ctype, stype, numstate, numparm, basetemp)" procedure in the "mak???()" routine that makes defines your channel. The "state" and "parm" fields are pointers to the channel's state and parameter set definitions that will be created automatically. "States" define the Markov channel's sequential states. The "vrev" parameter is the default reversal potential and will be used if vrev is not specified individually for the channel.
Several fields are used primarily for synaptic channels. The "trconc" field is a default parameter to set the sensitivity of a synapse to a ligand. Since a NeuronC synapse response is normalized (by the binding and saturation mechanism) to range between 0 and 1, the trconc paramter multiplies this value to set an absolute concentration of neurotransmitter for a synaptic Markov sequential state channel. Normally the "trconc" parameter is a "synaptic" parameter but its value will be taken from the channel's data structure if not specified individually for the synapse.
The "respamp" field provides a way to specify the relative response to a set of neurotransmitters or ligands for a channel. For each channel type that you want to respond differently to different ligands, you need to define the response amplitudes by placing them in the respamp array (see the "Markov synaptic" example below). At runtime, the channel's response will be calculated from a weighted sum of the ligand concentrations. These can be set by the "mesgout" parameter of a synapse, or by the "stim puff" statement.
The "perm" field specifies the relative permeability of the channel to different ions. An Na channel is slightly permeable to K ions, and a K channel slightly permeable to Na ions. Typically NMDA and AMPA channels carry Ca currents, with a fraction typically about 10-20%. These channels are connected to the Ca diffusion mechanism at a node (the "calcium compartment") and you can specify Ca flux and shell options for these channel types as you can with Ca channels.
For an HH-type channel (m3h, n4, ca3, etc), each state variable "m", "h", "c", "d", etc. needs to have an "implicit lookup table". The table is created automatically at run time, but you need define how many to make. This is done by defining a separate "parameter set" and a Q10 value for parameter "m", "h", "n", "d", "c", etc, and defining 2 "implicit lookup-tables" per parameter. You must set the number of implicit tables (see below.) For example, a m3h channel would have 2 "parameter sets", and a n4 channel would have only 1 parameter set.
For a Markov sequential-state channel types you normally create one "parameter set", but you would use more than one if you want to have different Q10 values for the some of the rate functions. Markov channel types don't use implicit lookup-tables so you don't need to set them up.
The "chancalc" and "funcname" fields are the default and user-defined functions that define the rate function. Both must be set in the "chanparm" structure, but if the user- defined function is not defined by the user, the rate comes from the default function. The "qt" and "dq" fields define what Q10 temperature coefficient are used. "qt" defines the Q10, and "dq" defines a starting temperature for the rate function definition (by default set to 22). The "oldv", "oldtiminc", and "oldtempcel" fields are temporary storage for the "v", "timinc", and "tempcel" variables when the table was most recently created. When these variables change, the lookup table for the parmeter may need to be regenerated.
The "dqca", "qcao" and "qcavoff" fields allow the simulator to keep track of how external Ca "cao" affects membrane potential. An increase in external Ca tends to remove the negative surface charge of the membrane and so tends to increase the polarization of the membrane, as if it were hyperpolarized. Thus an increase in "cao" causes Na channels to reduce their open probability. To keep track of this effect, the "dqca" field sets how many mV the membrane voltage is offset by a factor of 10 increase in Ca (normally about 18 mV but varies between channel types). The "qcao" field sets the base cao at which the channel rate functions were originally specified. The "qcavoff" is a temporary storage area for the voltage offset for each rate function ("m", "h", etc.) for a channel type.
struct chanparm { /* structure for channel rate params */
short int nfval; /* number of alpha, beta, etc. values */
short int nival; /* number of implicit values */
double *fval; /* alpha, beta, etc. vals for integ. */
double *ival; /* implicit values for integration */
float *chanftab; /* rate func tables */
float *chanitab; /* implicit tables if required */
double (*chancalc)(double,int); /* default rate function */
char *funcname; /* user-defined rate function */
double qrate; /* Q10 rate calculated from dq, qt */
double dq; /* Q10 factor for rate calc */
double qt; /* factor to normalize temp to 22 degC */
double dqca; /* Q10 for cao offset on chan voffset */
double qcao; /* base cao for chan offset: 2.3 mM */
double qcavoff; /* voltage offset calc from dqca,qcao */
double oldvf; /* Variables that need to be checked */
double oldvi;
double oldtimincf; /* to see if table needs to be */
double oldtiminci;
double oldtempcel; /* regenerated. */
double oldcao;
};
Each parameter set should have 2 values "nfval" for its alpha and beta functionns, and also 2 implicit values "nival" for the implicit integration method (see runcomp() in "ncomp.h"). You set the number of function and implicit values (and tables) for each parameter individually after the parameters are defined:
For a HH-type channel:
makna0()
{
nstate = 2; /* 2 Markov states (for noise) */
nparm = 2; /* make 2 sets of params, m, h */
ch=makchanconst(NA,0,nstate,nparm,dbasetc);/* make chan state const, parm */
ch->unitary = 6e-12; /* conductance measured at BASETC. */
ch->perm[PNA] = 1.0; /* rel permeability to Na+ ions */
ch->perm[PK] = dpkna; /* permeability to K+ ions / perm Na */
ch->vrev = vna; /* default reversal potential */
spnt = ch->state;
parm = ch->parm;
parm = ch->parm;
parm[0].nfval = 2; /* number of func vals (am, bm, etc.) */
parm[0].nival = 2; /* number of implicit vals (m1, m2) */
parm[0].chancalc = calcna0m; /* default rate function */
parm[0].funcname = "calcna0m"; /* user rate function */
parm[0].dq = dqm; /* Q10 for m rate function */
parm[1].nfval = 2; /* number of func vals (ah, bh, etc.) */
parm[1].nival = 2; /* number of implicit vals (h1, h2) */
parm[1].chancalc = calcna0h; /* default rate function */
parm[1].funcname = "calcna0h"; /* user rate function */
parm[1].dq = dqh; /* Q10 for h rate function */
mak2state(spnt,rchanf,rchanr,CHANRATE); /* make 2 Markov states for noise */
return ch;
}
For a Markov channel:
makna1()
{
nstate = 8; /* 8 Markov states */
nparm = 2; /* make 2 sets of params, m, h */
ch=makchanconst(NA,1,nstate,nparm,dbasetc); /* make chan state const, parms */
ch->unitary = 6e-12; /* conductance measured at BASETC. */
spnt = ch->state;
parm = ch->parm;
ch->perm[PNA] = 1.0; /* rel permeability to Na+ ions */
ch->perm[PK] = dpkna; /* permeability to K+ ions / perm Na */
ch->vrev = vna; /* default reversal potential */
parm[0].nfval = 2; /* number of func vals (am, bm, etc.) */
parm[0].nival = 0; /* number of implicit vals (m1, m2) */
parm[0].chancalc = calcna1m; /* default rate function */
parm[0].funcname = "calcna1m"; /* user rate function */
parm[0].dq = dqm; /* Q10 for m rate function */
parm[1].nfval = 2; /* number of func vals (ah, bh, etc.) */
parm[1].nival = 0; /* number of implicit vals (h1, h2) */
parm[1].chancalc = calcna1h; /* default rate function */
parm[1].funcname = "calcna1h"; /* user rate function */
parm[1].dq = dqh; /* Q10 for h rate function */
.
.
.
followed by a description of the states and their transitions
.
.
.
return ch;
}
Note the "mak2state()" procedure that is invoked by the
"makna0()" procedure above. It is called with the state pointer
"spnt", 2 rate functions (forward, reverse) and a rate
multiplier. If you use the "rchanf", "rchanr" functions, you'll
get a channel that is driven by either its own HH conductance
(i.e. for Na and K type 0, 2,3) or the conductance of its
controlling synapse (RSYN2, etc):For a Markov synaptic channel:
makampa1()
{
nstate = 7; /* 7 Markov states */
nparm = 1; /* make 1 set of params */
ch=makchanconst(AMPA,1,nstate,nparm,dbasetsyn);/* make chan const, parm */
ch->unitary = dampau;
ch->trconc = dstr;
ch->perm[PNA] = 2.0/3; /* rel permeability to Na+ ions */
ch->perm[PK] = DPK; /* rel permeability to K+ ions */
/* set to give vrev=0 */
ch->perm[PCA] = .1; /* permeability to Ca++ ions */
ch->vrev = 0; /* default reversal potential */
spnt = ch->state;
parm = ch->parm;
respamp = ch->respamp;
parm[0].nfval = 0; /* num of func vals (am, bm, etc.) */
parm[0].nival = 0; /* don't use implicit vals (yet) */
parm[0].chancalc = calcampa1; /* default rate function (dummy) */
parm[0].funcname = "calcampa1"; /* user rate function */
parm[0].dq = dqsyn; /* Q10 for m rate function */
respamp[GLU-GLU] = 0.8; /* response to glutamate */
respamp[AMPA-GLU] = 1.0; /* response to AMPA */
respamp[NMDA-GLU] = 1e-3; /* response to NMDA */
respamp[CNQX-GLU] = -1.0; /* response to CNQX */
.
.
.
followed by a description of the states and their transitions
.
.
.
return ch;
}
mak2state(spnt,rchanf,rchanr,CHANRATE); /* make 2 Markov states for noise */
struct hhchan { /* structure for membrane channel */
short int ctype; /* type of channel: na, k, etc. */
short int stype; /* subtype of channel: na, na2, etc.*/
hhchan *cnext; /* pointer to next channel */
float conduct; /* conductance of channel */
comp *comp1; /* channel compartment */
comp *comp2; /* possible calcium compartment */
float maxcond; /* maximum total channel conductance*/
float voffsm; /* activation offset voltage */
float voffsh; /* inactivation offset voltage */
float vrev; /* reversal potential */
chanconst *consts; /* channel constants */
char *cstate; /* channel random state */
float cafrac; /* fractional Ca conductance */
int nchan; /* total number of channels */
int cno; /* number of channels actually open */
float arate; /* activation rate */
float brate; /* inactivation rate */
float crate; /* deactivation rate */
float drate; /* re-activation rate */
short int numstate; /* number of states used */
stconc *conc; /* used for markov channels */
double m; /* m,n activation, range 0 - 1 */
double h; /* h inactivation range 0 - 1 */
} hhchan;
To add a new HH channel to "runcomp()", you add a new case
statement for it. You can modify the calculation of m*m*m*h
or n*n*n*n in "ncomp.cc" (look inside "runcomp()"). You can
also modify the calculation of the rate constants alpha and
beta for m,h, or n.Note that the original Hodkin-Huxley rate functions were defined at 6.3 deg Celsius, and ordinarily need to be normalized to the temperature at which a model is run. By default, HH rate functions are normalized to 22 deg assuming a Q10 for m of 2 and a Q10 for h of 3. Defaults for other Q10 values are also set to 3, allowing the default rate function definition to be for 22 deg C.
Each rate function is called with 2 parameters, 1) the membrane voltage and 2) the function subtype (alpha, beta, etc.). The function subtype must be less than or equal to the "nfval" field in the chanparm. To add a new subtype for a rate function, make a new "case" statement and put your new definition of the rate function after the case. Remember to place a "break" after the end of your addition.
Some simulators allow interactive modification of the rate functions. However, it is easy enough in "nc" to write a function and graph it.
The function "calcchanval(chan)" is used at run time to return implicit values for the HH equations. It computes function values by interpolating the implicit tables "chanival" and placing the values in "parm->ival". These interpolated values are generated automatically at run time for the alpha and beta rate constants from predefined tables. All the tables are generated at the start of a run by the "makchanconst(ctype,cnum,nstate,nparm)" procedure. Each channel subtype has its own individual tables generated by individual calls to the rate functions. This allows different subtypes to have different rate functions. The lookup tables are tabulated from -150 to +150 mV, in steps of 1 mV. At run time, they are interpolated with a simple linear appoximation algorithm.
The function "calcchanval(v-voffset, tau, parm)" is used at run time to return implicit values for the HH equations. It computes the value of one state variable ("m") by interpolating the implicit tables "chanival" and placing the values in "parm->ival". It must be called once for each state variable (i.e parameter set).
The sequential-state class calls "dochan()" and "dochan2()", defined in the top of "ncomp.cc". The sequential-state algorithm (in dochan2()) looks complicated but it is really very simple. Each sequential-state channel has a structure that defines it somewhat like the HH channels above but containing a little more information:
(taken from ncomp.h:)
struct chan { /* structure for membrane channel */
short int ctype; /* type of channel: na, k, etc. */
short int stype; /* subtype of channel: na, na2, etc.*/
hhchan *cnext; /* pointer to next channel */
float conduct; /* conductance of channel */
comp *comp1; /* channel compartment */
comp *comp2; /* possible calcium compartment */
float maxcond; /* maximum total channel conductance*/
float voffsm; /* activation offset voltage */
float voffsh; /* inactivation offset voltage */
float vrev; /* reversal potential */
chanconst *consts; /* channel constants */
float cdur; /* duration of channel event */
char *cstate; /* channel random state */
int nchan; /* total number of channels */
int cno; /* number of channels actually open */
float taum; /* activation tau */
float tauh; /* inactivation tau */
float tauc; /* deactivation tau */
float taud; /* re-activation tau */
short int numstate; /* number of states used */
stconc *conc; /* state concentr (variable size arr) */
};
struct stconc { /* structure for state concentration */
double cval; /* state concentr */
double cest; /* state concentr estimate */
double nchan; /* number of channels in this state */
};
Each channel state has a concentration and the
concentrations change according to transitions (defined in
"initchan.c") in proportion to alpha and beta, times a
constant, times the state concentration. The "stconc"
structure holds the concentration and delta for each state.
The number of states must be at least MINSTATE but may be as
large as MAXSTATE. Of course these limits may easily be
altered. Each individual channel has a "chan" structure
that defines it and also the appropriate number of "stconc"
structures that define the individual values of the state
concentrations. The channel also has a pointer to a table
(array) of "chanpar" structures. Each chanpar structure
contains an array of "chanstate" structures which define the
number of transitions, the states to which the transitions
are defined, the function returning the rate (alpha, beta),
a rate multiplier, and which compartment's voltage should be
used to calculate the rate:
(taken from "ncomp.h":)
struct chanstate {
/* structure for chan state constants*/
float cond; /* normalized conductance of state */
short int numtrans; /* number of transitions from state */
char trans[NUMTRANS]; /* four possible state transitions */
double (*trate[NUMTRANS])(); /* function returning rate for trans */
float ratemul[NUMTRANS]; /* rate multiplier for transition */
char rateo[NUMTRANS]; /* comp. for transition volt sens */
}; /* rateo=0 -> taua; rateo=1 -> taub */
To set up a new sequential-state channel, you can modify or
copy the "makna1()" or "makk1()" subroutines. The "makna3()"
subroutine is a copy of the "makna1()" subroutine, except that it
has 6 states instead of 8 and the transitions are a little
different. First you should draw the state diagram, and label the states with numbers starting with 0, and draw the allowable transitions. Then you set up the type, number of states, and a chanstate pointer "spnt".
nstate = 8; /* 8 Markov states */ nparm = 2; /* make 2 sets of params, m, h */ ch=makchanconst(NA,1,nstate,nparm,baset); /* make chan const, parm */ spnt = ch->state; parm = ch->parm;Then you set up the states. For each state, put in the number of allowable transitions, and the conductance. Then for each transition, (labeled starting with zero for each state) put in the state the transition goes to, the function returning the rate constant for the transition (e.g. "alm" or "alh", etc), a possible multiplier for the rate constant, and a number (either 0 or 1) defining which user-defined offset and rate to use (i.e. offsm or offsh) to pass to the rate function). A zero means use offsm.
m = 2;
state trans
spnt[0].numtrans = 1;
spnt[0].cond = 0;
spnt[0].trans [0] = 1;
spnt[0].trate [0] = alm;
spnt[0].ratemul[0] = 3.0 * m;
spnt[0].rateo [0] = 0;
spnt[1].numtrans = 2;
spnt[1].cond = 0;
spnt[1].trans [0] = 2; /* 0 <-> 1 <-> 2 <-> 3 */
spnt[1].trate [0] = alm;
spnt[1].ratemul[0] = 2.0 * m;
spnt[1].rateo [0] = 0;
spnt[1].trans [1] = 0; /* ^ | */
spnt[1].trate [1] = betm;
spnt[1].ratemul[1] = 1.0 * m;
spnt[1].rateo [1] = 0; /* 4 */
etc, etc.
Once you have defined the states and transitions, the
next change is to modify the alpha and beta rate functions.
You can use any of the existing ones but you can also make
up arbitrary new functions. To do this, see the discussion
of HH functions above.
For a different way to set up a sequential-state channel, look in "channa2.cc". The Na type 2 channel has 10 rate functions but some of them are "activation" parameters and others are "inactivation parameters". All the parameters that have the same Q10 are part of one "parameter set":
nstate = 9; /* 9 Markov states */ nparm = 2; /* make 2 sets of seq. state states */ ch = makchanconst(NA,2, nstate, nparm, baset); /* make chan const, parm */ ch->unitary = 6e-12; /* conductance measured at BASETC. */ spnt = ch->state; parm = ch->parm; ch->perm[PNA] = 1.0; /* rel permeability to Na+ ions */ ch->perm[PK] = dpkna; /* permeability to K+ ions / perm Na */ ch->vrev = vna; /* default reversal potential */ parm[0].nfval = 6; /* num of func vals (am, bm, etc.) */ parm[0].nival = 0; /* don't use implicit vals (yet) */ parm[0].chancalc = calcna2m; /* default rate function */ parm[0].funcname = "calcna2h"; /* user rate function */ parm[0].dq = dqm; /* Q10 for m rate function */ parm[1].nfval = 4; /* num of func vals (ah, bh, etc.) */ parm[1].nival = 0; /* don't use implicit vals (yet) */ parm[1].chancalc = calcna2h; /* default rate function */ parm[1].funcname = "calcna2h"; /* user rate function */ parm[0].dq = dqh; /* Q10 for m rate function */If you have designed a new channel type that you believe others would be interested in, please send the source code to allow others to share it.