Calcium () waves provide a complement to neuronal electrical signaling, forming a key part of a neuron’s second messenger system. We developed a reaction-diffusion model of an apical dendrite with diffusible inositol triphosphate (), diffusible , receptors (s), endoplasmic reticulum (ER) leak, and ER pump (SERCA) on ER. is released from ER stores via s upon binding of and . This results in -induced--release (CICR) and increases spread. At least two modes of wave spread have been suggested: a continuous mode based on presumed relative homogeneity of ER within the cell and a pseudo-saltatory model where regeneration occurs at discrete points with diffusion between them. We compared the effects of three patterns of hypothesized distribution: (1) continuous homogeneous ER, (2) hotspots with increased density ( hotspots), and (3) areas of increased ER density (ER stacks). All three modes produced waves with velocities similar to those measured in vitro (approximately 50–90 m /sec). Continuous ER showed high sensitivity to density increases, with time to onset reduced and speed increased. Increases in SERCA density resulted in opposite effects. The measures were sensitive to changes in density and spacing of hotspots and stacks. Increasing the apparent diffusion coefficient of substantially increased wave speed. An extended electrochemical model, including voltage-gated calcium channels and AMPA synapses, demonstrated that membrane priming via AMPA stimulation enhances subsequent wave amplitude and duration. Our modeling suggests that pharmacological targeting of s and SERCA could allow modulation of wave propagation in diseases where dysregulation has been implicated.
Calcium () is an important second messenger signal in many cell types, with diverse roles, from fertilization (Busa & Nuccitelli, 1985; Kretsinger, 1980) to regulating gene expression (West et al., 2001). is involved in triggering destructive processes including apoptosis (Orrenius, Zhivotovsky, & Nicotera, 2003) and ischemia (Lipton, 1999). Cells, including neurons, therefore regulate cytosolic concentration via buffers (Stern, 1992) and sequestration into mitochondria (Gunter, Yule, Gunter, Eliseev, & Salter, 2004) or endoplasmic reticulum (ER) (Berridge, 1998). In neurons, sequestration is modulated by neuronal activity (Pozzo-Miller et al., 1997), and elevated opens certain ion channels. There is therefore a bidirectional interaction between chemical signaling and electrophysiology (Blackwell, 2013; De Schutter & Smolen, 1998; De Schutter, 2008).
is heavily buffered but travels long distances (more than 100 m) to reach its targets. This poses a temporal problem if relying purely on diffusion. waves increase the rapidity of spread via -induced--release (CICR). CICR occurs in neurons and requires stores of held in the ER. ER is distributed throughout the cytosol in a connected way, through the dendrites and dendritic spines (Harris, 1994; Spacek & Harris, 1997). The ER SERCA pump pushes from the cytosol into the ER. When triggered by and , the ER receptors (s) open and release some of the ER’s into the cytosol. Regions of elevated then spread throughout portions of the dendritic tree (Ross, Nakamura, Watanabe, Larkum, & Lasser-Ross, 2005).
Understanding wave modulation is difficult since cytosolic changes over time and passes between different intracellular compartments and extracellular space (Fall, Wagner, Loew, & Nuccitelli, 2004; Wagner et al., 2004). Adding to this complexity is the nonuniform distribution of s, with clusters forming where local variations in or ER are heightened (Fitzpatrick et al., 2009).
At least two modes of wave spread have been identified: a continuous model that depends on continuous underlying substrate of regenerative potential and a pseudo-saltatory model where regeneration occurs at discrete points with diffusion between them. (We call it pseudo-saltatory here to distinguish it from classical saltatory conduction in myelinated fibers, involving capacitative effects in addition to electrodiffusion.) It has been hypothesized that these two modes produce downstream functional differences in dendrites, where many mechanisms are responsive to the level of , for example (Winograd, Destexhe, & Sanchez-Vives, 2008; Neymotin et al., 2013; Neymotin, McDougal, Hines, & Lytton, 2014) and synaptic plasticity (Kotaleski & Blackwell, 2010).
To investigate waves in a spatiotemporal context relevant for neurons, we developed a model of waves that includes cytosol and ER. Baseline cytosolic and concentration are set to low values. The model includes ER SERCA pumps (pump from cytosol into ER), leak channels ( leaks out of ER), and s. Our model generates waves with realistic physiological properties. We use our model to investigate how density and clustering, alterations in SERCA density, alterations in ER stacking, and diffusibility alter waves. An additional complexity arises when we consider coupling to plasma membrane calcium channels, which will contribute additional calcium flux triggered by rapidly spreading regenerative voltage changes on the membrane. We have therefore also assessed our results in an extended electrochemical model, which included a variety of ion channels (leak, voltage-dependent calcium channels: VGCC, potassium, sodium) with synaptic activation. We used this model to demonstrate how priming due to AMPA-mediated membrane depolarization would enhance subsequent wave amplitude and duration.
2 Materials and Methods
All simulations were run in the NEURON (version 7.3) simulation environment (Carnevale & Hines, 2006). NEURON has traditionally supported electrical modeling but has recently been extended to support reaction-diffusion (RxD) modeling as well (McDougal, Hines, & Lytton, 2013a, 2013b). The full 1D RxD calcium wave model of the neuronal dendrite (depicted in Figure 1) and the data analysis source code is available on ModelDB (Peterson, Healy, Nadkarni, Miller, & Shepherd, 1996; Hines, Morse, Migliore, Carnevale, & Shepherd, 2004) (http://senselab.med.yale.edu/modeldb).
Our dynamics are derived from Wagner et al. (2004), a spatial variant of Li and Rinzel (1994). Parameters are as in Table 1. We modeled a one-dimensional RxD system of intracellular neuronal waves in an unbranched apical dendrite of a hippocampal pyramidal neuron (length of 1000 m and diameter of 1 m). Within the dendrite, we modeled cytosolic and endoplasmic reticulum (ER) compartments by using a fractional volume for each. Suppose that for a given cell volume, fER denotes the fraction occupied by the ER (0.17) and fcyt denotes the fraction occupied by the cytosol (0.83). Necessarily . The inequality is strict if other structures are present, such as mitochondria.
ER-based calcium dynamics were previously modeled (De Young & Keizer, 1992). This previous work showed by a timescale analysis that the qualitative dynamics of receptors (s) could be represented only by considering the slow inactivation binding site state (Li & Rinzel, 1994). This work formed the basis of much subsequent work in intracellular dynamics (Fall & Rinzel, 2006; Fall et al., 2004; Hartsfield, 2005; Peercy, 2008; Wagner et al., 2004). The model presented here is a variant that neglects dynamic production (Wagner et al., 2004).
The ER model involves s and SERCA pumps. As with the plasma membrane model, we combine the net effects of all other channels on the ER into a leak channel. We denote by and the mass flux per unit volume due to the , SERCA pump, and leak channels, respectively. Dividing the mass flux by the volume fraction gives the change in concentration.
2.1 State Variables
2.3 Electrical Dynamics
Electrical dynamics were utilized for the subset of simulations described in the final figure. Electrical dynamics in the dendrite followed the standard parallel conductance model. Ion channels were based on a prior model of a hippocampal pyramidal cell (Safiulina et al., 2010) (http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=126814) and included L-, N-, and T-type channels, which allowed extracellular calcium to enter the dendrite at different voltage levels (McCormick & Huguenard, 1992; Kay & Wong, 1987). The equations for the ion channels follow: leak, Na, and K currents were represented by the conductance approximation: (g conductance; E reversal potential) using ENa = 50 mV; Ek = −77 mV; Eleak = −64 mV; (gleak = 39.4 10−6 S/cm), while the Goldman-Hodgkin-Katz (GHK) flux equation was used for currents: (p permeability). Channel dynamics were corrected for temperature by a Q using the factor of with T = 37C and 25C was taken to be the temperature at which the experiment was done. Conductances and activation curves were not corrected for temperature (Iftinca et al., 2006). Voltage-sensitive channels largely followed variants on the Hodgkin-Huxley formalism, whereby using steady-state value, , the forms were either, : or : , where x is m or n for an activation variable and h for an inactivation variable. F is Faraday’s constant; R is the gas constant. Variations on this scheme are noted below:
with = cm/s
hL (inactivation) was dependent: with ki = 0.001 mM
mL: ; , with , ; ; = 5
mT: , ; with ,
hT: , , with ; ; a0 = 0.04; = 5.
with = cm/s;
mN used ; ; with ;
hN following: ; ; constant = 80; with ; a0 = 0.03, = 5.
The model also contained a transient sodium channel INa and a delayed rectifier channel to allow for action potential generation, a calcium-dependent potassium channel that hyperpolarized the cell after calcium influx, and an A-type potassium channel for rapid inactivation. Equations for these channels follow:
with = 0.11 S/cm
mNa using with:
hNa with special form and using with and ; = 2. sNa:
where = 0.01 S/cm;
following an atypical steady-state: with using the same in with ; a0 = 0.02; = 1
with = 0.009 S/cm; oKCa
with = 0.48 10−3; = 0.13 10−6; using .
where = 0.02 S/cm
nA with and an atypical steady-state value: with ; ; ; a0 = 0.05; = 5. lA: with steady-state: ; ; time constant: ; = 1.
The AMPA receptor was modeled with a double-exponential mechanism: rise 0.05 ms, decay 5.3 ms, E 0 mV. A single AMPA synapse was placed at 500 m, midway on the dendrite; synaptic weight = 0.5 . The metabotropic synapse was not modeled in detail but was given as a local 12.5× increase in (1.25 M) in a 4 micron segment at the same location.
2.5 Simulation Variations
In one set of simulations, we changed the density of and SERCA by scaling and , respectively. hotspots were modeled by scaling at discrete locations, while keeping density in between the hotspots at 0.8× baseline density. ER stacks were simulated by scaling , , and simultaneously, while keeping ER density in between the stacks at 0.8× baseline density. Changes in buffering were simulated by scaling the diffusion coefficients, and , in equal amounts. We also modulated in a set of simulations. The range of parameter values examined was constrained to generate waves with wave features (onset, speed, duration, amplitude) that were around ranges for experimentally observed limits in neurons (Nakamura, Barbara, Nakamura, & Ross, 1999; Ross et al., 2005; Fitzpatrick et al., 2009).
2.6 Electrochemical Simulation Variations
One set of simulations (see Figure 11) was run with the electrical dynamics (see section 2.3) to determine how AMPAergic stimulation’s alteration of membrane potential prior to an stimulus affects the subsequent waves. Simulations were run for 12 s. Background concentration in these simulations was set between 0 and 0.01 mM to prevent spontaneous oscillations. AMPA activations (150 spikes with interspike interval of 25 ms) finished 250 ms prior to stimulus (amplitude: 2.5 mM).
2.7 Data Analysis
and concentrations were recorded in a 3D array using temporal and spatial resolution of 5 ms and 1 m. Cytosolic wave features were extracted by thresholding the array using an amplitude threshold of 2× baseline cytosolic concentration (2.5× for the electrical simulations). Wave onset was defined as the delay after stimulus until the passed threshold. Amplitude was maximum concentration. Speed was calculated using onset time at stimulation location and final onset time at location of wave termination (note that the use of final onset time leads to artifactual lowering of wave speed in pathological condition of multiple waves). Duration was calculated as median at each position from wave onset to offset (going below threshold).
This study involved over 4800 15 to 30 second simulations, testing how variations in levels of density, SERCA density, hotspots, and ER stacks, altered wave initiation, amplitude, speed, and duration. An additional set of 32 12-second simulations were run, testing how the number of AMPA inputs prior to an stimulus affected the calcium waves. Simulations were run using the NEURON simulator on Linux on a 2.27 GHz quad-core Intel XEON CPU. Thirty seconds of simulation ran in 1.5 to 2.0 minutes, depending on simulation type.
We simulated the arrival of from a metabotropic receptor activation cascade as an instantaneous augmentation of concentration to 1.25 M (12.5× background) in a 4 micron segment in mid-dendrite (see Figure 2). then spread gradually along the dendrite, providing a permissive effect for activation of the receptors (s). Activation of s permitted release of from the endoplasmic reticulum (ER) stores into the cytosol (see Figure 2a). The elevation in cytosolic was mirrored by a depression in ER as the wave passed (see Figure 2b). The cytosolic wave started once adequate had built up to spill over and activate neighboring s and initiate sustained positive feedback of induced--release (CICR). The wave then spread bidirectionally as and simultaneously diffused laterally. waves were not able to reverse direction and propagate back toward their source due to an refractory period provided by inactivation gating. Parameters including (1.415 m/ms), (0.08 m/ms), (0.08 m/ms), (120400.0 molecules/mM/ms), (1.9565 molecules/ms), and (18.06 molecules/mM/ms) were adjusted from modeling (Wagner et al., 2004) and experiments (Allbritton, Meyer, & Stryer, 1992), to produce a wave speed of 77 m/s, which is comparable to neuronal wave speeds measured experimentally (6822 m/s in Nakamura et al., 1999). These parameters provided our baseline simulation for further parameter explorations. Median elevation lasted about 1 s and peaked at 1.6 M, both also comparable with experimental results (Ross et al., 2005; Fitzpatrick et al., 2009).
3.1 Density Modulates Excitability
was necessary for wave generation (see Figure 3). waves did not begin below density of 92.2% of baseline, when the first wave was obtained (see Figure 3a and left red dots on Figures 3c to 3f). This wave had initiation delays of about 100 ms longer compared to the baseline wave (see Figure 2a and the middle red dot on Figure 3c), with only a slightly lowered amplitude, lower speed, and lower duration. The major effect, time to onset, was strongly dependent on due to the need to gradually source a sufficient amount of to produce enough lateral spillover to initiate a positive feedback cycle, as the lower density produced a lower total flux from ER to cytosol. ER must also be sourced at a rate that exceeds the reuptake governed by SERCA pumping. With increased density, time to onset decreased rapidly to a low of 25 ms.
Speed increased proportionately with increased (see Figure 3d). At the lowest value, wave speed was 72 m/s. Speed depended on both the rate of ER sourcing at each point of wave regeneration and the rate of diffusion required to reach the next set of s laterally. Proportionate increase in speed with density was due to the reduced requirement for additional for positive feedback, since less is needed to trigger the release of . At each location, must source adequate not only for positive feedback and spillover, but also to exceed the local sink back into ER provided by the SERCA pump. wave duration and amplitude both had a small positive association with increases in (see Figures 3e and 3f), due to the higher density of s liberating more for the wave, which then remained elevated longer.
3.2 SERCA Opposes Effects
The SERCA pump returns back into the ER, opposing the release from conductance. Therefore, the parameter dependence of measures was generally opposite to that seen with density (see Figure 4). The least excitable wave was produced at the highest SERCA density tested (see Figure 4b). The most excitable wave shown here, at 0.66× baseline SERCA (see Figure 4a), was caused by diminished reuptake after release. Further reduction in SERCA resulted first in spontaneous waves and then in simultaneous release from all ER sites at once due to regenerative responses starting from baseline levels (not shown).
Patterns of wave measures were similar to, but reversed from, those seen with increase. Onset to the first wave varied from 30 to 330 ms. Shorter onset times were produced by faster accumulation with less reuptake. Similarly, the most excitable wave had a faster propagation speed (84 m/s) and longer duration due to longer-lasting supporting faster spread. Amplitude showed a slight inverse relationship with increasing SERCA density due to the higher SERCA reuptake of into ER stores diminishing availability for the wave (see Figure 4f).
At the maximal SERCA consistent with wave propagation, 1.07× baseline, onset to wave initiation was increased dramatically (330 ms), wave speed decreased to 75 m/s, and duration decreased to 0.82 s (see Figure 4b). waves were no longer able to initiate at higher levels as was sucked back into the ER before a wave could ignite.
3.3 / SERCA Balance Regulates Wave Induction
The balance between and SERCA density determined whether waves could be elicited via stimulation (see Figure 5). Excessive or insufficient would not allow a wave (regions 1,4). Maximum level in region 1 was low compared to regions 2 and 3, since ER had already partially leaked out of the ER prior to the stimulus and then diffused within the dendrite, preventing the large, localized rise in cytosolic , which formed a wave.
Moving from region 1 toward region 2, density decreased and SERCA density increased, both of which reduced cytosolic , which was no longer sustained above threshold throughout the simulation. This allowed the stimulus to elicit distinct waves. However, the high levels of produced multiple waves without return to a subthreshold state. These repetitive waves caused an artifactual lowering of measured wave speed (see section 2 for details of calculation). amplitude was higher than in region 1 due to repeated higher elevations of cytosolic compared to balanced fluxes in region 1.
In region 3, there was a match across and SERCA density, allowing for a single wave to be elicited from the stimulus. Within this region, for any given SERCA density, an increase in density increased flux from ER to cytosol, reducing time to onset and increasing speed, duration, and amplitude. In region 4, SERCA dominated the dynamics, shuttling back into the ER too quickly to allow stimulation to initiate a wave.
3.4 Diffusibility Regulates Wave Onset and Speed
Altering the diffusion constant () had a pronounced effect on wave initiation (see Figure 6). With low (0.1415 m/ms), the stimulus remained localized in the central stimulus region (approximately 500 m), producing a high prolonged local elevation with more rapid activation of local s and shorter wave onset time (40 ms). near the stimulus then remained elevated longer than at other dendritic locations. Since the stimulus spread to neighboring dendritic locations more slowly, wave speed was slightly reduced (73.5 m/s).
High (1.981 m/ms) significantly delayed the onset of the wave (230 ms; see Figure 6b) due to diffusing quickly away, producing less local stimulation and increasing the time required to trigger the wave. However, once the wave was initiated, it spread slightly faster (77.6 m/s), the prearriving giving a slight boost at subsequent locations along the dendrite. Amplitude had a minimal inverse relationship with increasing (see Figure 6f) because local spread out more quickly, reducing local activation. At values larger than 1.981 m/ms, was so diffuse that it could not elicit the wave.
3.5 Pseudo-Saltatory Waves via Hotspots
The prior simulations were performed using uniform density of ER mechanisms at all locations. However, there is evidence of inhomogeneities, for example, at dendritic branch points, where elevations in local density ( hotspots) might contribute to assisted propagation of waves at these sites of potential failure (Fitzpatrick et al., 2009). In neurons, hotspots have average center-to-center spacing of approximately 20 m (edge-to-edge spacing of 10 m) but show considerable variability in the spacing (Fitzpatrick et al., 2009). In the following simulations, we varied hotspot density while keeping the inter-hotspot density at 80% of baseline.
Both hotspot strength and spacing altered wave speeds (see Figure 7). We started with a 20% reduction in density between hotspots (see Figure 7a), with the hotspots having 93% of the density used in Figure 2. This alteration reduced wave speed from 77 to 68 m/s. However, the propagation pattern differed qualitatively, with its saltatory nature readily seen as spots of high concentration, which occur at the locations of hotspots. Local velocity parallels amplitude in showing an increase at the hotspots, where activation produces a local highly varying gradient (peaked high second spatial derivative) leading to a higher immediate diffusion speed. Slower wave progression occurs between hotspots where lower concentration only slightly boosts the wave progression from diffusion. Increased density at hotspots increased wave speed further to 90 m/s (2.0× baseline; see Figure 7b). Further augmentation provided an approximately linear increase in speed due to faster and larger release of from the ER stores at the hotspots (see Figure 7c). Below about 93% of hotspot density, the waves did not initiate, since there was insufficient release to trigger waves (0 speeds at the left of Figure 7c). Wave propagation also did not occur in the absence of between hotspots, with the threshold for wave initiation having inter-hotspot and hotspot densities of approximately 0.66 and 0.93× baseline, respectively.
Hotspot spacing also modulated wave speeds, with larger spacing producing slowing as and had to travel further via mildly boosted diffusion between hotspots before being fully reboosted. At 1.87× baseline density with a spacing of 15 m, the wave had a speed of 100 m/s (see Figure 7d). Increasing center-to-center spacing between hotspots to 100 m resulted in a reduction of wave propagation speed to 66 m/s. At the densities shown, hotspot spacing had a larger impact on wave speed (about 66 to 100 m/s) than density (about 68 to 90 m/s).
Intracellular concentration is heavily regulated via buffering mechanisms, in part presumed to provide careful regulation of -triggered signaling cascades. buffering also modulates diffusion efficacy and apparent diffusion coefficient (). Because we were not modeling buffering directly in these simulations, we altered instead. Changing dramatically altered propagation speed over a wide range (see Figure 8). With diminished , wave speed was substantially reduced to 25 m/s (see Figure 8a). The converse of this was that was relatively immobile, with duration at one location slightly increased (see Figure 8e). This heightened local elevation also allowed for a shorter onset to wave initiation (20 ms; see Figure 8c). Increasing above baseline level had opposite effects: wave speed was augmented (288 m/s) but time to onset was delayed (40 ms), both due to faster spread of . We also note that speed showed much greater alteration than onset.
3.6 Pseudo-Saltatory Waves Via ER Stacks
There are at least two ways that hotspots could occur in dendrites: type 1 (increased density of at particular locations on homogeneous distribution of ER) and type 2 (increased “density” (accumulation) of ER at particular locations). The former case was explored in the prior section. The latter case has been identified as locations of ER lamellar specialization that are sometimes described as ER stacks. We next explored such stacks as an alternative type of hotspot, noting that this increase in local ER provides increased density of SERCA and leak as well as at the hotspot locations. In these simulations, the ER density between stacks was at 0.8× baseline, consistent with continuous ER throughout the dendrite (Martone et al., 1993; Terasaki, Slater, Fein, Schmidek, & Reese, 1994).
Increasing ER stack density from about 0.8× (see Figure 9a) to about 2× baseline (see Figure 9b) increased wave speed from 68 m/s to 86 m/s (see Figure 9c). This is due to more release of from the leak and channels, which was opposed by heightened SERCA activity. The duration of amplitude elevation was reduced at the ER stacks. In addition, as the ER stack density increased, the onset to wave initiation shortened (220 to 30 ms) due to heightened leak and extrusion of into cytosol. However, elevation duration decreased as the ER stack density increased (965 to 795 ms) due to heightened SERCA pumping.
Increasing the center-to-center spacing of ER stacks while maintaining a fixed ER density of 1.86× baseline (see Figures 9d to 9f) tended to decrease wave speed (93 to 71 m/s). This was due to lower overall availability of ER for releasing and spreading the wave. Interestingly, with heightened center-to-center spacing (15 and 100 m) of ER stacks, the duration of heightened elevation was increased (755 to 960 ms) due to lower overall presence of SERCA pumps.
Waves could not initiate below a minimum of 0.8× ER stack density due to insufficient sourcing. The effects of density and spacing of stacks were generally similar to those seen with hotspots. Onset time depended primarily on ER stack density. Variation of with ER stacks (see Figure 10) produced results very similar to those seen with type 1 hotspots (see Figure 8), except that the effects on time to onset were more pronounced. With ER stacks, the onsets were significantly larger (35–70 ms) compared to those for hotspots (20–40 ms). These heightened onset times with ER stacks were due to the higher SERCA activity, which reduces cytosolic availability.