In every cell type, intracellular \({Ca}^{2+}\) signaling plays a major role in a diverse set of activities, including process homeostasis and genetic expression. In neurons, \({Ca}^{2+}\) wave propagation has been implicated as a major component of intracellular signaling and is complementary to membrane electrical signaling via action potentials. It follows that the cell must contain tightly-controlled intrinsic mechanisms for regulating intracellular \({Ca}^{2+}\) content and waves; these include cytosolic buffers, storage in the ER, and sequestration into mitochondria. In this tutorial, we examine how some of these basic mechanisms interact to produce a \({Ca}^{2+}\) wave.
\({Ca}^{2+}\) is well-buffered by several cytosolic proteins (e.g. calmodulin, calbindin, etc...), but must still traverse a significant distance within the neuron to reach its targets. This distance is long enough that pure diffusion would pose a temporal problem for cell communication. Instead, the cell uses oscillating \({Ca}^{2+}\) waves to harness the power of chemical signaling. These waves can travel up and down the Endoplasmic Reticulum (ER) to effect homeostatic changes. Understanding how these waves are modulated in real time by the cell poses a great difficulty because of complex spatiotemporal patterns and nonuniform distribution of the channels that control \({Ca}^{2+}\) flux across the ER. Inositol triphosphate (\(IP_3\)) receptors (\(IP_3Rs\)), Sarco/endoplasmic reticulum \({Ca}^{2+}\)-ATPase (SERCA pump), and Ryanoine receptors (RyRs). For simplicity’s sake, we will not deal with RyRs in this tutorial.
Diffusible \({Ca}^{2+}\) is released directly from the ER via \(IP_3Rs\) upon cooperative binding of free \({Ca}^{2+}\) and \(IP_3\) to the receptor. When \(IP_3\) levels are ideal, this results in \({Ca}^{2+}\)-induced \({Ca}^{2+}\) release (CICR), yielding a wave that can travel along the ER. As mentioned above, these waves are implicated in the gene expression and protein upregulation associated with neuronal plasticity by way of their interaction with the neuronal soma. The SERCA pumps within ER are normally responsible for maintaining appropriate intracellular \({Ca}^{2+}\) levels. Two models exist for wave propagation: continuous and saltatory. To understand the physical and functional differences between these models, it is ideal to compare them using Reaction-Diffusion (RxD) wtihin the NEURON environment.
The model presented in this tutorial generates \({Ca}^{2+}\) waves with realistic physiological properties, including wave speed, amplitude, duration, and onset.
Python is a rich and user-friendly programming language that offers significant options for systems integration. For this reason, we will run the NEURON and RxD modules within the Python interpreter.
Like other RxD problems within NEURON, the first step is to import and load the appropriate libraries into the Python interpreter. Importing “gui” will open up the GUI for graphical manipulation of the hoc programming environment.
from neuron import rxd, h, gui
First, we define parameters for the section to be used. Like any imported module within python, the imported hoc and rxd libraries can be accessed by the “h” and “rxd” prefixes:
sec = h.Section()
sec.L=100
sec.diam=1
sec.nseg=100
The above code defines a section of a neuron with Length = 100 \({\mu}m\), diameter = 1 \({\mu}m\), and 100 discrete segments. By chunking the neuron into a large number of segments, we increase the resolution of the output at the expense of processing power. For a smaller number of segments, the opposite is true. In larger models, it is useful to represent “nseg” as a variable – this way, if changes need to be made during any number of simulations, the single variable can be altered and will affect all subsequent iterations of “nseg”. In this example, we only reference a single section, so there is no need to elaborate.
Next, we activate the time-step integration algorithm. This is a class of functions which are part of the SUNDIALS package and allow variable-order, variable time-step, and implicit solving of ordinary differential equations (ODE). This is a particularly useful tool for solving a system with multiple species, whose concentrations may change suddenly from stable conditions.
h.CVode().active(1)
The argument of active specifies whether or not the CVode solving system will be active during the program. When active, this line allows finitialize() and fadvance() to call two additional functions respectively: CVode.re_init() and CVode.solve(). To activate CVode, the samp:{arg} = 1 as in the above example. This calls the CVode solver for a single time-step.
CVode.re_init() initializes the function integrator, and CVode.solve() will integrate the differential equations by one time-step. Essentially, this step primes the program for action.
Note
CVode.solve(tout) can specify integration until the step passes the time-step specified by variable tout
Now we can set initial parameters for the cell component variables:
Baseline diffusion constants for \({Ca}^{2+}\) and \(IP_3\) are numbers confirmed by prior biological experiments. These rates represent the diffusive movement after buffering by native cytosol proteins is taken into account. In this tutorial, some numbers have been optimized to produce an output \({Ca}^{2+}\) wave that corresponds well to biological observations. Below is an example of some of the parameters used that are necessary in the tutorial model:
caDiff = 0.016
ip3Diff = 0.283
cac_init = 1.e-4
ip3_init = 0.1
gip3r = 12040
gserca = 0.3913
gleak = 6.020
kserca = 0.1
kip3 = 0.15
kact = 0.4
ip3rtau = 2000
A region is the space within which reactions occur.
{cyt} = rxd.Region({h.allsec()}, nrn_region='{i}', geometry=rxd.FractionalVolume({fc},surface_fraction={1}))
This uses the rxd.Region class to define the properties for the region on the immediate interior of the membrane (nrn_region = "i"). Electrophysiologically, the only regions that matter for this simulation are the spaces directly inside the membranes. The “i” marks that we are concerned with the inner region only. For any region, different geometries can be specified; here we model the region as fractional volumes. The geometry argument tells NEURON what to expect the region to actually look like (i.e. a cylinder has a specific 3D spatial representation and therefore requires specific diffusion modeling for ions within that particular geometry).
For this model, we make the assumption that ER is roughly evenly-mixed throughout the neuron segment. By similar logic, we set surface_fraction=1 because 100% of the cytosol is in contact with the ER membrane (i.e. \({Ca}^{2+}\) must enter the cytosol before it enters the ER – there is no direct connection between the ER and the extracellular fluid). This is in contrast to models with significant ER stacking or ER in close proximity to the plasma membrane, in which cases this fraction would necessarily be less than 1.
er= rxd.Region(h.allsec(), geometry=rxd.FractionalVolume(fe))
We do the same for the ER region, although this time it is not necessary to specify the “surface_fraction” because we are inside the ER membrane and there is no contact with the cytosol:
Where the rxd.ScalableBorder(scale, on_cell_surface=False) takes two arguments. The scale is the scaling factor for the diameter of the region. In this case, it is fine to use “1”, although in other cases, with complicated regional geometry, various scaling factors can be specified.
Note
The older code for this was rxd.FixedPerimeter(num), which was a bit of a hack, because it did not allow for scaling.
Taken together, the above code specifies a perimeter region for the ER membrane geometry that scales proportionally by diameter in all sections. This allows us to simplify many of the calculations that deal with density of \(IP_3\) receptors on the ER membrane. If the area is 1, then it is easy to calculate #units/area.
It is important to remember that the rxd.Region class specifies reactions that occur in a single spatial compartment. Therefore, the variables er and cyt and cyt_er_membrane each specify single regions, as far as the model is concerned. This is somewhat intuitive from looking at the arguments. We will see shortly that there is another class of functions to deal with reactions across compartments, termed rxd.MultiCompartmentReaction.
Definition of a species in RxD follows the general form:
var = rxd.Species(regions, d, name, charge, initial)
Because these are keyword arguments, they can be left out or even entered in any order. If the arguments are given in the order they are specified in __init__.py for the rxd.Species class (i.e. regions, d, name, charge, initial) then they keywords can be omitted. For more on this, see the older Python documentation on keyword arguments.
Below is an example of how rxd.Species might be used:
{ca} = rxd.Species([cyt, er], d={caDiff}, nam={ca}, charge={2}, initial={caCYT_init})
Note that the following two pieces of code are equivalent to the above line:
ca = rxd.Species([cyt, er], caDiff, ca,2,caCYT_init)
ca = rxd.Species([cyt, er], name=ca, d=caDiff, initial=caCYT_init, charge=2)
Notice that the variable ca is defined in two separate regions. Most of the arguments reference parameter values that may be specified at the beginning of the file in a parameter table.
In this tutorial, we will also define the state of the \(IP_3R\) gate and the slow \({Ca}^{2+}\) inactivation gate (the “h” gate). The rxd.State class is functionally similar to rxd.Species for now, but this will likely change in the future:
ip3 = rxd.Species(cyt, d=ip3Diff, initial=ip3_init)`
ip3r_gate_state = rxd.State(cyt_er_membrane, initial=0.8)
h_gate = ip3r_gate_state[cyt_er_membrane]
The rxd.MultiCompartmentReaction uses angle brackets to delineate the basic reaction direction. For instance, the code is interpreted thusly: “>” points from reactant towards product, and “<>” indicates a bidirectional reaction. Bidirectional reactions necessarily must have two rate constants, which must be listed in succession. Therefore:
ip3r = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], k, k, membrane=cyt_er_membrane)
The rxd.MultiCompartmentReaction class requires a membrane as a reaction parameter. In this example, cyt_er_membrane was called from an earlier definition, as part of a region within the reaction. Note that this class of functions specifies the movement of a particular species across compartments that have been previously defined using rxd.Region. Therefore, ca[er] refers to the concentration of \({Ca}^{2+}\) within the ER compartment. You might think that rxd.MultiCompartmentReaction is a general class, with rxd.Region being a simpler version. This is not the case. Reactions across multiple compartments must necessarily scale by membrane area (or other specified region area) to properly accomodate geometries.
Therefore, care must be used in determining which class of functions to use when defining a new species.
The rxd.MultiCompartmentReaction class also contains an argument for custom_dynamics. The SERCA pump and leak channel code examples below exemplify this nicely:
serca = rxd.MultiCompartmentReaction(ca[cyt]>ca[er], gserca/((kserca / (1000. * ca[cyt])) ** 2 + 1), membrane=cyt_er_membrane, custom_dynamics=True)
leak = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], gleak, gleak, membrane=cyt_er_membrane)
If custom_dynamics is set to False then the reaction assumes laws of mass action and will multiply the rate by appropriate powers of species. Otherwise, it will use the equations given to determined rate. In this example, the SERCA pump is modeled with Hill-type dynamics, and the leak channel with Hodgkin-Huxley dynamics. The SERCA pump is active and acts against a gradient. Therefore, it has dynamics that are not strictly mass action and custom_dynamics=True, which accesses another bit of code within the rxd.MultiCompartmentReaction module that enables this environment to be properly computed.
Here are some additional definitions included in the model:
minf = ip3[cyt] * 1000. * ca[cyt] / (ip3[cyt] + kip3) / (1000. * ca[cyt] + kact)
k = gip3r * (minf * h_gate) ** 3
ip3r = rxd.MultiCompartmentReaction(ca[er]<>ca[cyt], k, k, membrane=cyt_er_membrane)
ip3rg = rxd.Rate(h_gate, (1. / (1 + 1000. * ca[cyt] / (0.3)) - h_gate) / ip3rtau)
Where minf refers to \(m_{\infty}\), which is the m-gate open probability at steady state.
The rxd.Rate class is discussed elsewhere <> but will be reviewed here briefly, for clarity. Objects specified by this class take the form:
r = rxd.Rate(var, reaction_rate, regions=None)
Here, we specify the rate dynamics for regions containing var and variables contained in reaction_rate. By setting regions=None (or by ommitting this parameter entirely), the code is interpreted as occurring at all possible regions that have been previously specified.
Within the RxD environment, variable information will be fed into the time-step integration as the concentrations change in real time. Effectively, this sets up the proper framework for running the reaction.
Let’s use the \(IP_3R\) gating as an example. We can model the rate of the \(IP_3R\) activity by:
We define \({h_{gate}}\) as h_gate, which is the slow gate of the \({Ca}^{2+}\) channel. Therefore, we end up with:
ip3rg = rxd.Rate(h_gate, (1./ (1+ 1000.* ca[cyt] / (0.3)) - h_gate)/ ip3rtau)
{ip3rg} = rxd.Rate({h_gate}, (1. / (1 + 1000. * {ca}[cyt] / (0.3)) - {h_gate}) / {ip3rtau})
This line uses the rxd.Rate class to specify what is changing and the equation for the rate of change. Note that this includes a value for ip3rtau, which is a time constant for the opening of the \(IP_3R\) gate.
Warning
Because of unit conversion issues, many of the variables must be mutliplied or divided by factors of 10^3. Don’t be alarmed when you see this in the equations to be solved!
Once we have finished setting up all the parameters (phew!) we can begin the code to run the program:
h.finitialize()
This command (run via the hoc module) initializes all sections and point processes. Whatever intial values have been specified previously will now be set. This is the first line of the code that executes the wave.
cae_init = (0.0017 - cac_init *fc) / fe
ca[er].concentration = cae_init
These lines set the initial \({Ca}^{2+}\) concentration in the ER based on the fractional volumes specified earlier and a value for the average cytosolic \({Ca}^{2+}\) concentration (‘0.0017’).
What about density of \(IP_3\) along the ER membrane? Well, there is code for that as well:
for node in ip3.nodes:
if node.x < 0.2:
node.concentration = 2
First, this loop accesses every node along the ER, where the number of nodes is determined by the value of nseg specified earlier. In NEURON, every section that is specified has a fractional length value that can be specified by .x, where x is a value from 0-1. This preserves directionality for a given connectivity of compartments. Therefore, the loop specifies that for the leftmost 20% of every node, the \(IP_3\) density should be set to 2 (which, relative to the normal value of 1, means a 100% increase in \(IP_3\) density). Because nodes have properties of space, they have discrete locations and can be identified within NEURON. This is particularly useful for modeling saltatory wave conduction; experimental evidence has shown that \(IP_3\) and \(IP_3\) Receptors are not uniformly distributed along the ER membrane within a given neuron.
Now we re-initialize the integrator. We must do this, because state variables change with the variable step methods and keep their own local copies:
h.CVode().re_init()
When using variable timesteps, it may be prudent to manually set an absolute error value when dealing with low concentrations (e.g. \({\mu}m\), or 1e-6). This can be done by using the h.cvode.atol(arg) command. When no argument is specified (i.e. h.cvode.atol()) the program will use the default local absolute error tolerance. The default value can also be set by specifying the value as the argument. Functionally, this command forces the CVode solver to set the local error below a value defined by the function:
rtol*|state| + atol*atolscale_for_state
As might be expected, the local relative error tolerance can be set by a similar command, h.cvode.rtol(arg). See the CVode documentation page on the NEURON website for more on this.
We also need code to save and print the output of our simulation:
def save():
print 't = %g cai = %g' % (h.t, sec.cai)
s.flush()
s.printfile('figs/wave1d_%g.ps' % h.t)
for t in xrange(0, 10000, 100):
h.CVode().event(t, save)
We will be looking at snapshots from 0-10000 msec at time-steps of 100 msec. And now the final piece:
h.continuerun(10000)
This final command tells the program to run until the time-step reaches 10,000 ms. The value of 10000 was chosen for clarity in demonstrating the changes in \({Ca}^{2+}\) concentration for this model. Obviously, this parameter can be modified as needed for other models.
Sample output shows \({Ca}^{2+}\) concentration increasing and then decreasing over time. This output can easily be plotted in a variety of different ways to produce graphs of the wave, or can be stored in arrays for vector analysis when multiple simulations are run:
t = 0, cai = 0.0001
t = 100, cai = 7.69054e-05
t = 200, cai = 6.26621e-05
t = 300, cai = 5.48893e-05
t = 400, cai = 5.10401e-05
t = 500, cai = 4.93336e-05
t = 600, cai = 4.89449e-05
t = 700, cai = 4.92243e-05
t = 800, cai = 4.97612e-05
t = 900, cai = 5.03961e-05
t = 1000, cai = 5.17513e-05
t = 1100, cai = 5.87617e-05
t = 1200, cai = 9.46684e-05
t = 1300, cai = 0.000233663
t = 1400, cai = 0.000564351
t = 1500, cai = 0.000961326
t = 1600, cai = 0.00126194
t = 1700, cai = 0.00144737
t = 1800, cai = 0.001544
t = 1900, cai = 0.0016034
t = 2000, cai = 0.00164564
t = 2100, cai = 0.00166171
t = 2200, cai = 0.00166515
t = 2300, cai = 0.00166192
t = 2400, cai = 0.0016552
t = 2500, cai = 0.00164667
t = 2600, cai = 0.00163655
t = 2700, cai = 0.00162452
t = 2800, cai = 0.00161084
t = 2900, cai = 0.00159581
t = 3000, cai = 0.0015794
t = 3100, cai = 0.00156158
t = 3200, cai = 0.00154239
t = 3300, cai = 0.00152189
t = 3400, cai = 0.00150008
t = 3500, cai = 0.00147698
t = 3600, cai = 0.00145261
t = 3700, cai = 0.00142702
t = 3800, cai = 0.00140023
t = 3900, cai = 0.00137228
t = 4000, cai = 0.0013432
t = 4100, cai = 0.00131304
t = 4200, cai = 0.00128187
t = 4300, cai = 0.00124974
t = 4400, cai = 0.0012167
t = 4500, cai = 0.00118283
t = 4600, cai = 0.00114817
t = 4700, cai = 0.0011128
t = 4800, cai = 0.00107678
t = 4900, cai = 0.00104017
t = 5000, cai = 0.00100305
t = 5100, cai = 0.00096545
t = 5200, cai = 0.00092745
t = 5300, cai = 0.000889132
t = 5400, cai = 0.000850495
t = 5500, cai = 0.000811596
t = 5600, cai = 0.000772436
t = 5700, cai = 0.000733064
t = 5800, cai = 0.000693477
t = 5900, cai = 0.000653636
t = 6000, cai = 0.000613516
t = 6100, cai = 0.00057314
t = 6200, cai = 0.00053243
t = 6300, cai = 0.000491127
t = 6400, cai = 0.000449115
t = 6500, cai = 0.000406105
t = 6600, cai = 0.000362066
t = 6700, cai = 0.000316843
t = 6800, cai = 0.000270977
t = 6900, cai = 0.00022551
t = 7000, cai = 0.000183225
t = 7100, cai = 0.000145728
t = 7200, cai = 0.000114834
t = 7300, cai = 9.03681e-05
t = 7400, cai = 7.39168e-05
t = 7500, cai = 6.2258e-05
t = 7600, cai = 5.1986e-05
t = 7700, cai = 4.68036e-05
t = 7800, cai = 4.66802e-05
t = 7900, cai = 4.28365e-05
t = 8000, cai = 3.89536e-05
t = 8100, cai = 3.86087e-05
t = 8200, cai = 3.93548e-05
t = 8300, cai = 4.03637e-05
t = 8400, cai = 4.18398e-05
t = 8500, cai = 4.36184e-05
t = 8600, cai = 4.44633e-05
t = 8700, cai = 4.52969e-05
t = 8800, cai = 4.69812e-05
t = 8900, cai = 4.85148e-05
t = 9000, cai = 4.96922e-05
t = 9100, cai = 5.03194e-05
t = 9200, cai = 5.06684e-05
t = 9300, cai = 5.13789e-05
t = 9400, cai = 5.19757e-05
t = 9500, cai = 5.22195e-05
t = 9600, cai = 5.27178e-05
t = 9700, cai = 5.34486e-05
t = 9800, cai = 5.40841e-05
t = 9900, cai = 5.47133e-05