What Is Required for Neuronal Calcium Waves? A Numerical Parameter Study.

Neuronal calcium signals propagating by simple diffusion and reaction with mobile and stationary buffers are limited to cellular microdomains. The distance intracellular calcium signals can travel may be significantly increased by means of calcium-induced calcium release from internal calcium stores, notably the endoplasmic reticulum. The organelle, which can be thought of as a cell-within-a-cell, is able to sequester large amounts of cytosolic calcium ions via SERCA pumps and selectively release them into the cytosol through ryanodine receptor channels leading to the formation of calcium waves. In this study, we set out to investigate the basic properties of such dendritic calcium waves and how they depend on the three parameters dendrite radius, ER radius and ryanodine receptor density in the endoplasmic membrane. We demonstrate that there are stable and abortive regimes for calcium waves, depending on the above morphological and physiological parameters. In stable regimes, calcium waves can travel across long dendritic distances, similar to electrical action potentials. We further observe that abortive regimes exist, which could be relevant for spike-timing dependent plasticity, as travel distances and wave velocities vary with changing intracellular architecture. For some of these regimes, analytic functions could be derived that fit the simulation data. In parameter spaces, that are non-trivially influenced by the three-dimensional calcium concentration profile, we were not able to derive such a functional description, demonstrating the mathematical requirement to model and simulate biochemical signaling in three-dimensional space.

Introduction

Intracellular calcium signals define a transition point between electrical signals and biochemical responses in neurons. While basal calcium concentrations in the cytosol are very low, neurons can modulate local cytosolic calcium concentrations to induce microdomain calcium signals, which can integrate to produce longer ranging signal propagation towards the soma. There they reach the nucleus and trigger gene transcription responses relevant for learning [1-4] and neuroprotection [5-8]. Cellular calcium signals are shaped by calcium transport mechanisms embedded in the plasma membrane, across which calcium can be bi-directionally exchanged between cytosol and the extracellular space (see [9] for an overview). In addition, intracellular organelles like mitochondria and the endoplasmic reticulum (ER) function as large calcium stores [10-16]. Organelle membranes are equipped with calcium exchange mechanisms that allow transport of calcium from the cytosol into the organelles or vice versa. Some of these mechanisms—notably the ryanodine receptor channel (RyR) in the ER membrane, which is able to release large amounts of calcium from the ER into the cytosol—have a positive feedback property: Their opening is facilitated by the presence of cytosolic calcium and will trigger the release of even more calcium through surrounding channels. This calcium-induced calcium release mechanism (CICR) overcomes the limited reach of purely diffusive calcium signals in a buffered regime. CICR in neurons has been implicated in various neurodegenerative diseases [17-23] and has been studied with respect to CICR wave-like properties [24-38]. The importance of CICR wave-like dynamics therefore merits thorough investigation. One major question is how the cellular and intracellular architecture can shape and regulate calcium signals that are able to propagate over long distances, for instance in the context of synapse-to-nucleus communication. It was shown in [26] that variations in the activation of Ca2+ store release through inositol 1,4,5-trisphosphate (IP3) receptors were able to significantly change the calcium wave patterns. Additionally, the distances between distinct IP3 receptor clusters and the pump strength were shown to control wave stability and instability [37, 38]. Sequestration properties through SERCA pumps [36] and mitochondrially controlled CICR [27] were further implicated in CICR. Related work has shown that calcium signals can modulate the shape of the cell nucleus [39, 40]. Thus, there seem to exist mechanisms in which calcium signals shape the geometry of organelles and the intracellular architecture shapes calcium signals.

In this study, we therefore chose to investigate how the relationship between ER and dendrite geometry can influence the dynamics of CICR waves. The ER is typically thought of as a continuous extension of the outer nuclear envelope that forms the rough ER surrounding the cell nucleus. As smooth ER, it extends into dendrites, all the way into synaptic spines [41]. Reconstructions of portions of smooth ER in dendrites have shown the ER to define a moderately branched, continuous structure [41, 42]. From an idealized standpoint, the ER has previously been considered as a neuron-within-a-neuron [43-45], approximated in form of a cable that is embedded in a larger cable (the dendrite). We decided to follow this path of morphology simplification and designed perfectly symmetrical ER/dendrite cable-in-cable model domains (see Fig. 1).

Fig. 1

Simulation domain and model components. The domain for all simulations of this study was a cylindrical dendrite 50 μm in length with variable radius containing a centrally positioned cylindrical ER of variable radius. The rotational symmetry of the domain was used to reduce the problem to two dimensions (axial and radial position). The calcium model contains calcium in the cytosol and the ER as well as calbindin (CalB) in the cytosol. The dynamics of both are governed by a diffusive process and a buffering reaction. Calcium can cross the ER membrane through RyR channels and SERCA pumps, the plasma membrane through PMCA and NCX pumps. Each simulation was initiated by a 1 ms calcium influx through the left end of the cytosol

The associated calcium model consists of a diffusion-reaction partial differential equation for the propagation of calcium in the cytosol and the reaction of free calcium with the buffer calbindin. Endoplasmic calcium is modeled by a simple diffusion equation. Membrane exchange mechanisms are included via flux boundary conditions at the plasma membrane and the ER membrane: Ca2+-ATPase pumps, Na+/Ca2+ exchangers, as well as a leakage term on the plasma membrane; sarco-/endoplasmic reticulum Ca2+-ATPase pumps, ryanodine receptor channels and a leakage flux on the ER membrane.

Our study focused on: (1) how calcium dynamics depend on the parameters dendrite radius, ER radius, and ryanodine receptor density, (2) whether thresholds can be found for stable calcium waves and (3) on gathering information about wave velocities. Simulation results show how cellular architecture, even in this simplified scenario, can determine the shape and fate of calcium waves and their stability. These data were further used to derive empirical laws that explain the behavior of the observed calcium dynamics in distinct geometric regimes. As a side note and in agreement with [26], we also show that the higher spatial dimension (as compared to classical 1D cable models) is indeed required to adequately represent even the most basic properties of calcium wave propagation.

Model and Methods

Model Domain

We conducted all simulations on a perfectly cylindrical model dendrite with a fixed length of 50 μm and variable radius, containing a cylindrically shaped ER of variable radius positioned exactly at the center of the dendrite; see Fig. 1. We used the rotational symmetry to reduce the problem to two dimensions. Calcium signals were induced by a calcium influx density through the left boundary of the dendrite with an initial strength of that linearly decreased to zero within 1 ms.

Model Equations

Spatio-temporal calcium dynamics in the intracellular space are modeled by a system of diffusion-reaction equations described in the following. Transmembrane currents through channels and pumps are incorporated into these equations as flux boundary conditions (cf. Sect. 2.3).

Our model includes cytosolic/endoplasmic calcium and, in the cytosol, calbindin-D28k (CalB) as a calcium buffer. In the equations, we will represent cytosolic and endoplasmic calcium concentrations by the symbols , , and the (unbound) buffer concentration by b, respectively. Mobility in the cytosol/ER is described by the diffusion equation where represents both calcium and calbindin as functions in space and time.

The buffering reaction between cytosolic calcium and calbindin is described by The fact that CalB has four distinct high-affinity calcium-binding sites [46] is taken into account by quadrupling its concentration for the simulation (thus implicitly assuming identical and independent binding sites). The full domain equations for cytosolic calcium and calbindin are thus given by in the cytosolic domain, where the concentration of the CalB-Ca2+ compound is expressed by the difference of the total concentration of CalB present in the cytosol () and free CalB, the former of which is assumed to be constant in space and time (this amounts to the assumption that free and calcium-binding CalB have the same diffusive properties). All diffusion and reaction parameters are listed in Table 1.

Table 1

Model parameters and initial values

Initial and equilibrium values
c_c	50 nM	(chosen)
c_e	250 μM	(chosen)
c_o	1 mM	(chosen)
b^tot	40 μM	[47]
Diffusion/reaction
D_c	220 μm^2s^−1	[48]
D_b	20 μm^2s^−1	[49]
κb−	19 s^−1	[47]
κb+	27 μM^−1s^−1	[47]
RyR channel
ka−	28.8 s^−1	[50]
ka+	1500 μM^−4s^−1	[50]
kb−	385.9 s^−1	[50]
kb+	1500 μM^−3s^−1	[50]
kc−	0.1 s^−1	[50]
kc+	1.75 s^−1	[50]
IRref	3.5×10^−18mols^−1	[51] (approx.)
SERCA pumps
I_S	6.5×10^−21molμMs^−1	[52], (adapt.)
K_S	180 nM	[24]
PMCA pumps
I_P	1.7×10^−23mols^−1	[53]
K_P	60 nM	[54]
ρ_P	500 μm^−2	(estim.)
NCX pumps
I_N	2.5×10^−21mols^−1	[53], (adapt.)
K_N	1.8 μM	[53]
ρ_N	15 μm^−2	(estim.)
Leakage
v_l,e	38 nm s^−1	(calc.)
v_l,p	4.5 nm s^−1	(calc.)

Membrane Transport Mechanisms

In order to study calcium waves, we included ryanodine receptors for calcium release from the ER. We also added sarco-/endoplasmic reticulum Ca2+-ATPase pumps (SERCA) for re-uptake as well as a leakage flux on the ER membrane. Moreover, we defined fluxes across the plasma membrane by addition of Ca2+-ATPase pumps (PMCA), Na+/Ca2+ exchangers (NCX) and a plasma membrane leakage term. This amounts to the ER and plasma membrane flux density equations Here, is the RyR flux density, the SERCA flux density and the leakage flux density on the ER membrane. The flux densities of PMCA, NCX and the leakage flux density of the plasma membrane are denoted by , and , respectively. All flux densities are oriented towards the cytosol so that positive signs mean influx into the cytosol. Besides RyR, ER-membrane-bound IP3 receptor (IP3R) channels are another source for intracellular calcium. They are known to require calcium and, more importantly, IP3 in the cytosol to become permeable to luminal calcium. IP3 is produced at the plasma membrane by receptor activation of phospholipase C. Compared to RyR calcium efflux, IP3R calcium dynamics take place on a slower time scale (governed by the diffusion of IP3). Preliminary simulations for the present study have shown that IP3Rs may contribute to the initiation of a traveling wave, when small calcium influxes are not enough to trigger a wave through RyRs alone. The wave dynamics are, however, dominated by RyR. We therefore decided to exclude IP3R-mediated calcium dynamics.

To incorporate existing single-channel or single-pump models in the partial differential equations formulation, we assume channels and pumps to be continuously distributed along the membranes, using their membrane densities to calculate the flux densities and .

RyR Channels

The calcium flux density through ryanodine receptor channels in the ER membrane is given by an expression of the form where is the density of RyR in the ER membrane, is the open state probability of a single channel, and the single-channel Ca2+ current. We describe the single channel ionic current by where the reference current is approximated from data presented in [51].

The open probability for RyR channels is modeled according to [50] and can be calculated as the sum of the two open states and in the system of ordinary differential equations emerging from a four-state Markov model with the kinetic constants , and .

SERCA Pumps

The current from sarco-/endoplasmic reticulum Ca2+-ATPase pumps is described by a model from [24], which was adapted for the three-dimensional case, and gives rise to the Ca2+ flux density The model reflects the dependence of the Ca2+ current not only on the cytosolic concentration, but also on the endoplasmic saturation. The density of SERCA pumps in the ER membrane was adapted to the RyR channel density in each simulation to ensure a zero net flux through the membrane in equilibrium conditions.

An interesting question with regard to the modeling of trans-membrane currents is their equilibration at stationary, i.e., resting concentrations. The usual approach is to calibrate a leakage current such that it exactly counters the combined net current of all other involved mechanisms at equilibrium. At the ER membrane, this is only possible as long as the resulting direction of the leakage current points outward. With varying RyR densities, however, this is not necessarily the case (as leakage and RyR current have the same direction). Since the SERCA current is the only inward current through the ER membrane that we consider in our model, we chose to calibrate their density to ensure a zero net flux at equilibrium. This results in variable SERCA densities (larger when the RyR density is larger), which appears plausible from an energy consumption point of view: An increase in leakage currents to counter a reduction in RyR density would only lead to an increase in “futile cycling” through leakage and SERCA pumps. Another approach to avoid unneccessary consumption of energy is to use a bi-directional SERCA model (cf. [55], also for a more detailed view on “futile cycling”). Our simulations showed that the exact equilibration of the membrane currents has only a minuscule impact on the results presented here—the threshold values presented in the Results section only change by about 1% with equilibration by leakage flux and constant SERCA density (data not shown). The reason is that the propagation of the wave front takes place on a much faster time scale than the re-uptake of calcium through SERCA pumps. The same applies to the extrusion by PMCA and NCX pumps, which are therefore modeled in a basic way.

NCX and PMCA Pumps

For the Na+/Ca2+ exchanger current, our model assumes a constant Na+ concentration at the plasma membrane. The current densities for both plasma membrane transport mechanisms are expressed as first- and second-order Hill equations, respectively (cf. [53]):

Leakage

Both ER and plasma membranes allow a leakage flux not accounted for by the above transport mechanisms. The plasma membrane leakage flux is calibrated to ensure zero membrane net flux in the equilibrium state. Leakage flux densities are modeled by where is the extracellular Ca2+ concentration, which is assumed to be constant throughout all simulations. Values for all model parameters are gathered in Table 1.

Numerical Methods

For numerical simulations, the three equations are discretized in space using a finite volumes method. This allows for a natural integration of current densities across ER and plasma membranes into the reaction-diffusion process.

The system of ordinary differential equations (in time) arising from this procedure is nonlinear (due to the nonlinear reaction term and, more importantly, the highly nonlinear transport terms across the membranes). As these equations, notably the RyR channel dynamics, require a precise calibration of the time step size used for numerical solution, we employed a linearly-implicit extrapolation (LIMEX) scheme to solve the nonlinear system, since LIMEX offers automated error-estimate-based control of integration order and step size [56, 57].

For the results we present here, the linear problems emerging in the LIMEX method were solved using a Bi-CGSTAB [58] linear solver preconditioned by a geometric multigrid method using Gauss–Seidel smoothing and SuperLU [59] as base solver. Computations were facilitated by a domain decomposition parallelization approach and carried out using the UG 4 framework [60] on the JURECA computer system at the Jülich Supercomputing Centre [61].

Implementation

All model components were implemented in a NeuroBox project. NeuroBox [62] is a neuroscientific toolbox that combines 1D, 2D and 3D modeling and simulation of electrical and biochemical signaling in a visual workflow environment. Visual workflows are created with VRL-Studio [63] and the general-purpose numerical framework UG 4 [60] is used to solve the set of coupled nonlinear partial differential equations.

Results

Range of Calcium Waves

In a first series of simulations, we examined how calcium dynamics depend on the ER radius and the density of RyR channels in the ER membrane. To that end, we created five model dendrite domains with a fixed dendrite radius of 0.2 μm each and an ER radius of 40 nm, 50 nm, 60 nm, 70 nm, and 80 nm, respectively. We then simulated calcium dynamics with RyR densities ranging from 0.5  to 4.0  on all five domains and measured how far the initial calcium signal was able to travel through the dendrite.

The propagation distance was defined as the distance between the left end of the dendrite (the activation point) and the rightmost position where the open probability of the RyR channels exceeded 0.1 at any point in time. In our study we observe both stable and abortive calcium waves. The qualitative differences are illustrated in Fig. 2. For a stable wave (depicted for dendrite radius 0.4 μm, ER radius 0.15 μm, RyR density 2.5  in Fig. 2A), the traveling wavefront is clearly visible and maintains a constant shape. Reduction of the ER radius to a value of 0.11 μm makes the wave regime abortive (Fig. 2B): The wave front does not propagate further than 20 μm and the concentration at the wave front diminishes over time. While the velocity of the abortive wave decays to zero, it remains constant in the stable case (Fig. 2C)—justifying the assumption that the stable wave can propagate arbitrarily far. Snapshots of a sample simulation of a stable wave are depicted in Fig. 2D).

Fig. 2

Propagation of stable and abortive calcium waves. (A) Axial calcium profiles of a stable calcium wave at different points in time (dendrite radius 0.4 μm, ER radius 0.15 μm, RyR density 2.5 ). The depicted cytosolic calcium profiles are recorded directly at the ER membrane and over the whole length of the model dendrite. The shape of the wave remains constant and travels from left to right in a convective manner, although driven by a reaction-diffusion process with nonlinear calcium exchange across the ER. (B) Axial calcium wave profiles of an abortive calcium wave (same setup as in A, but with a smaller ER radius of 0.11 μm). While a wavefront traveling from left to right is clearly visible, it has a smaller amplitude than in A and breaks down before reaching the far end of the dendrite. (C) Velocity of the calcium wave fronts in A and B as a function of time. The stable wave quickly reaches a constant velocity of about 1.06  and travels through the whole dendrite at this speed, while the velocity of the abortive wave peaks at about 0.98  and then declines to zero while the wave breaks down. (D) Sample simulation of a 1 μm dendrite with elicitation of a stable calcium wave. The radial coordinate is scaled by a factor of 8 to enhance visibility

We plotted the propagation distance as a function of RyR density in Fig. 3. The traces for ER radii of 50 nm, 60 nm, 70 nm and 80 nm rise slowly at first, but then very rapidly approach a singularity at a RyR density that is specific to each domain. This threshold RyR density decreases with increasing ER radius and also exists for an ER radius of 40 nm, but is outside the plotted RyR density range in Fig. 3.

Fig. 3

Range of calcium waves. Simulations on a cylindrical model dendrite (radius 0.2 μm) containing a cylindrical ER (radius 40 nm, 50 nm, 60 nm, 70 nm, 80 nm) reveal that there is a threshold RyR density (at the singularities of the plotted curves) above which calcium waves are stable. Decreasing the RyR density below the threshold induces abortive waves, whose travel distances rapidly decline with RyR density reduction. Larger ER domains lead to a lowering of the stabilizing RyR density threshold (singularity moves to the left), indicating that the ability to sustain a calcium wave depends on the rate at which calcium can be released from the ER

We draw the following conclusions from these observations for our 0.2 μm model dendrite: In summary, the interplay between ER radius and RyR density enables dendrites to modulate the stability and reach of calcium signals, which may be critical in synaptic cross-communication and synapse-to-nucleus signaling.

A minimal RyR density is required in all five domains to elicit a calcium wave. The sharp gradient in the curves depicted in Fig. 3 show that once a critical RyR density is surpassed, calcium waves become stable. Otherwise the initial calcium influx causes only a very localized transient that does not propagate along the dendrite. This result is comparable to minimal distance criteria for the IP3 receptor clusters studied in [37].

Dendrites with smaller ER need a higher RyR density to elicit a wave, indicating that the ability to trigger a calcium wave scales with the rate of calcium release from the ER. This seems to agree with [37], where wave stability is controlled by IP3R pump strength and, as mentioned above, by IP3R cluster distances. The critical threshold for wave stability decreases with increasing ER radius. Since dendrite diameters typically decrease the further they are away from the soma, ER diameters are forced to decrease as well. This would lead to a higher probability for calcium waves to become abortive. Figure 3 demonstrates how such dendritic tapering can be compensated for by an increase in RyR density.

Conversely, as calcium release from the ER increases both with ER radius and RyR density, there is also a threshold ER radius for any sufficiently large RyR density above which calcium waves are stable. Thus, larger ER compartments in dendrites will enable neurons to more readily induce long-distance calcium signals.

There are regimes in which abortive calcium waves are initiated (closely below threshold density). Such threshold dependency has previously been introduced in [26]. Our results confirm the existence of such thresholded regimes and are later quantified by deriving empirical laws. While the general behavior will be either a micro-domain calcium event or a calcium wave, there is a critical intermediate region in which abortive waves with variable travel distances can be elicited.

Threshold for Stable Calcium Waves

Having found the existence of such a threshold, we examined the prerequisites for stable calcium waves more closely. The guiding questions for this investigation were: To address these questions, we varied the dendrite radius between 0.1 μm and 1.0 μm and the RyR channel density in the ER membrane between 0.5  and 4.0 . For various combinations of these two parameters, we determined the threshold ER radius, above which a stable calcium wave can be elicited, using a bisection technique: We first located an interval in the ER radius domain for which the lower bound did not allow for a stable calcium wave, while the upper bound did. Then we simulated the calcium dynamics with ER radius taken from the middle of the interval and repeated the procedure with the lower half of the interval, if a stable wave was detected, with the upper half otherwise, until an adequate precision was reached. For this purpose, we accepted as a stable wave any signal, which, upon activation by calcium influx through the left-hand side of the model dendrite, propagated all the way (i.e., 50 μm) along the dendrite and was detected, in a measurement zone located at the right-hand side (see also Fig. 2). With regard to the steep, singularity-like gradients of the wave range/RyR density functions for ranges between 20 μm and 50 μm in Fig. 3, we are confident that this is an appropriate choice. Since even if a setting existed where a wave had a range of 50 μm but became abortive beyond that range, a minor increase in RyR density (or ER radius) would then substantially increase its range.

Do such thresholds also exist for larger dendrites?

If so, what is the relationship between dendrite radius and ER radius (and RyR density) at the threshold?

With a look at the results in Fig. 4, the first question can be answered affirmatively. In fact, for all combinations of dendrite radius and RyR density that we checked, we found an ER radius threshold above which stable calcium waves were produced. This threshold ER radius decreases with increasing RyR density (indicating that the ability to sustain a stable calcium wave scales with the rate of calcium release from the ER) and increases with growing dendrite radius (because there is more space that the released calcium can diffuse into, which reduces the effective near-ER-membrane concentration). It is interesting to point out three more aspects:

Fig. 4

Threshold ER radius for eliciting stable calcium waves. For any tested combination of RyR channel density and dendrite radius, there exists a threshold ER radius above which stable calcium waves can be elicited. Overall, the threshold ER radius increases with increasing dendrite radius (since there is more space that the released calcium can diffuse into, which reduces the effective near-ER-membrane concentration) and decreasing RyR density (indicating that the ability to sustain a stable calcium wave scales with the rate of calcium release from the ER). (A) Threshold ER radius as a function of RyR channel density in the ER membrane. The lower four traces, corresponding to the lower dendrite radii, exhibit similar quasi-hyperbolic behavior, scaled by the dendrite radius. They can be fitted by a functional description that assumes instant radial distribution of released calcium in the small radial range from ER to plasma membrane, cf. “Empirical threshold laws”. The upper traces, corresponding to larger dendrites, show two separate regimes: Separated to the left, where the distance between ER and plasma membrane is small, they converge towards a limit trace to the right, where the two membranes are far away from each other. This convergence can be attributed to wave propagation that is faster than radial diffusion to the plasma membrane, cf. again “Empirical threshold laws”. (B) Threshold ER radius as a function of dendrite radius (like A, but projected onto perpendicular plane). While the traces are approximately linear in a small dendrite regime, they reach a limit threshold for larger dendrites. No such limit is reached for the lowest RyR density tested. This is due to calcium buffering in the cytosol. Thus, ER size does not need to be increased proportionally with dendrite radius. Neurons could therefore ensure stable calcium waves, even within larger dendrites, where intracellular space is occupied by other organelles

The traces in Fig. 4A (threshold ER radius plotted against RyR density) appear to reach a lower limit for large RyR densities, meaning the ER needs to have a minimal size to produce stable calcium waves. This makes sense as it has to contain a minimal amount of calcium that can actually be released during a wave event. This result is in agreement with [36], where an increase in SERCA density led to an increase in luminal calcium. This was demonstrated to have a positive feedback on the wave dynamics.

The traces in Fig. 4B (threshold ER radius plotted against dendrite radius) reach an upper limit for large dendrite radii. In other words, there is a limit to the thickness of the ER required to elicit stable calcium waves and the ER does not need to be bigger than that in, e.g., very thick proximal dendrites. Figure 5 depicts the limit ER threshold values as a function of RyR density. The data is fitted by a model function of the form with , and .

Fig. 5

Limit threshold ER radius as function of RyR density. The ER radius r above which stable calcium waves can be elicited for arbitrarily large dendrites (for a given RyR channel density ρ, cf. Fig. 4) can be fitted by a function of the form . The parameter values , and have been determined using a least squares fitting. The value of c indicates that such a limit ER radius exists only for RyR densities larger than approximately 0.64 . This is due to calcium buffering in the cytosol: Simulations with the cytosolic buffer concentration reduced by 90% result in a fitting curve with . Neurons could therefore modify RyR densities according to intracellular buffering properties to then maintain a minimal ER “footprint” for stable calcium waves

The trace corresponding to the lowest tested RyR density (0.5 ) in Fig. 4B does not reach such an upper limit (the trace corresponding to a RyR density of 1.0  does reach a limit, but outside the depicted range). One can conclude that a minimal RyR density is required for stable calcium waves (see also [26]) if the ER is restricted to specific dendritic volumes. Further investigation of this observation showed that this is due to calcium buffering in the cytosol: If the RyR-supported calcium current through the membrane is too small, too much of the released calcium is buffered and there is not enough left to trigger further release through RyR channels. This results is in agreement with [37], where similar conclusions were made for IP3R-induced waves. Reducing the buffer concentration by a factor of ten allowed for a limit ER threshold radius for all tested RyR densities (cf. Fig. 5).

Empirical Threshold Laws

To complete the investigation of our second question, we tried to describe the relationship between threshold ER radius and dendrite radius/RyR density by deriving appropriate empirical laws for the traces in Fig. 4. The fact that we found a threshold ER radius for wave elicitation in all combinations of dendrite radius and RyR density suggests that there exists a function (r for ER radius, R for dendrite radius, ρ for RyR density) describing a threshold ER radius manifold in the three-dimensional parameter space spanned by r, R and ρ.

The spatio-temporal patterns of calcium waves can be roughly separated into three groups by means of the distance between ER and plasma membrane. The first group is defined by morphological configurations in which this distance is small enough for the radial calcium concentration profiles to be approximately constant, the second group contains the cases in which the distance is large enough for it not to influence the wave propagation, and the third group comprises all intermediate cases. In this section, we will derive separate empirical threshold laws for the former two groups. To motivate this distinction, we introduce a little thought experiment (cf. Fig. 6A).

Fig. 6

Schematics for deriving empirical threshold laws of wave stability. (A) In a stable wave regime the axial position x, where a cytosolic threshold calcium concentration is attained near the ER membrane, needs to move by Δx to the right within a given time Δt (defined by the wave velocity). The amount of calcium required to increase the original concentration at needs to be released from the ER. (B) In dendrites where the distance between ER and plasma membrane is small (compared to the diffusive range of calcium for time Δt), the iso-concentration surfaces (indicated by blue-to-red sequence of lines) are nearly radially aligned. This means the released calcium is roughly constant in radial direction. (C) In dendrites where the distance between ER and plasma membrane is large (compared to the diffusive range), the iso-concentration manifolds (blue-to-red sequence of lines) are bent backwards in axial direction. Calcium released from the ER only diffuses to an effective radius within the wave propagation time Δt and the wave front does not “sense” the plasma membrane

In a stable wave regime, the axial position x, where a cytosolic threshold calcium concentration is attained near the ER membrane, needs to shift by Δx to the right within a certain time Δt (governed by the wave velocity). The amount of calcium necessary to increase the original concentration at to within this time span needs to be released from the ER. It is proportional to the number of releasing RyR channels (see also [26]), which, in turn, is proportional to: (i) the RyR density ρ and (ii) the radius r of the ER membrane, both assumed constant in time. The single-channel current, however, is time-dependent in that it depends on locally available endoplasmic calcium (see Eq. (8)), which decays rapidly during RyR release. While the amount of released calcium is proportional to r, the endoplasmic volume is proportional to . Assuming rapid radial equilibration within the ER, the decay of a single-channel current is proportional to . Integrating the decreasing current from t to shows that the amount of released calcium is also proportional to (iii) the factor with some constant b, in a first-order approximation. Putting (i) to (iii) together, an amount proportional to is released into the cytosol, This is where we distinguish between regimes with small and large distances between plasma and ER membrane.

Small Distance Between Plasma and ER Membrane

When the distance between ER and plasma membrane is small, then there is only a negligible radial concentration gradient at all times, as diffusion will almost instantaneously produce constant radial concentration profiles (cf. Fig. 6B). With this in mind, the increase in cytosolic calcium concentration, caused by the amount of calcium released from the ER, is inversely proportional to , assuming that a fixed percentage of the released calcium is bound by the mobile buffer, which leaves proportionality intact. As the required increase to bridge the gap between the concentrations at x and is fix, a constant can be introduced to arrive at the relation with the parameters b and c, usable in a least squares fitting, which very closely fits the 0.1 μm dendrite radius trace in Fig. 4A. The one for the 0.2 μm dendrite radius trace is also acceptable, but using R as an additional parameter in the fitting expression produces substantially better fitting results for the 0.3 μm and 0.4 μm dendrite radius traces, indicating that the model assumption of radially constant calcium concentrations is already largely violated in theses cases.

Large Distance Between Plasma and ER Membrane

When the distance between ER and plasma membrane is large, the wave front moves faster in the axial direction than its radial component reaches the plasma membrane, and the wave does not “sense” the presence of the plasma membrane (cf. Fig. 6C). The fact that the upper traces in Fig. 4A converge towards a limit trace and that the lower traces in Fig. 4B attain a limit is due to this circumstance. In this situation, radial calcium diffusion will be limited to a distance a from the ER membrane until the threshold concentration is reached at axial position . It is reasonable to assume that the calcium concentration near the ER membrane after release through RyR channels is inversely proportional to the volume it diffuses into, i.e., . Like above, with a fixed percentage of the released calcium being bound by the buffer, we can introduce a constant (which depends on the buffer concentration) to arrive at the equation The parameters a, b, and c can be determined to fit the limit threshold ER radii attained for large dendrite radii (Fig. 5, blue trace), also for buffer concentrations reduced by a factor of ten, in which case the parameter c is significantly smaller (Fig. 5, orange trace). Note that the buffer concentration puts a constraint on the RyR density: If the density falls below a threshold (the constant c), then there can be no calcium wave, regardless of how big the ER is.

Intermediate Distances

We were not able to find descriptive model functions to fit the threshold values for intermediate distances. In such regimes, the significant radial concentration gradient at the wave front, which varies strongly with the distance between the two membranes, makes it difficult to find expressions similar to the two presented ones.

Wave Velocity

During the described simulations, we observed that the velocity of calcium waves varies with ER radius and RyR density. This led us to finally examine the velocity of stable waves. Following the initial activation, all stable waves reached a constant velocity shortly after their formation (cf. Fig. 2C)). These constant velocities were measured using the rightmost axial coordinate where the RyR open probability exceeded 0.1 as the wave front position in every time step of the simulation. Results (shown for a dendrite radius of 0.2 μm in Fig. 7) indicate that calcium wave velocity increases both with increasing ER radius and RyR density. This can be expected given that the driving force for the movement of the wave is diffusion, which becomes faster when axial concentration gradients become bigger, which, in turn, is the case when the efflux density from the ER increases. The traces for wave velocity as a function of RyR density rise more rapidly in a small range above the threshold for stable waves and then exhibit a quasi-linear behavior with a slope seemingly independent of the ER radius. Velocity traces for a larger dendrite of radius 1.0 μm show the same basic behavior (data not shown), though with a less pronounced dependency on the ER radius, which is in line with our findings for the ER radius threshold in the previous section.

Fig. 7

Stable calcium wave velocity. The velocity of stable calcium waves depends on the RyR channel density in the ER membrane and on the ER radius. The traces for wave velocity as a function of RyR density rise quickly in a small range above the threshold for stable waves and then exhibit a quasi-linear behavior with a slope seemingly independent of the ER radius. In dendrites where multiple synaptic sites are activated in a spatio-temporal manner, this linear dependence could result in linear additivity of multiple calcium wave velocities and thus may be useful in describing multi-wave interactions in more complex dendritic branching patterns

Discussion

Our computational study finds that calcium waves can be triggered in an idealized morphological setting. If certain geometric and physiological conditions are met, wave propagation is a stable process, similar to the all-or-nothing properties of electrical action potentials. A minimal ryanodine receptor density is required to guarantee sufficient calcium release from the ER and this density scales inversely with ER surface. We also observed abortive regimes in small parameter regions below thresholds for ER size and ryanodine receptor densities. When considering realistic morphologies, spatial variations in dendritic and ER morphology as well as changes in RyR density are likely to occur, thus, abortive regimes may be important when it comes to synaptic cross-communication with respect to calcium signals.

We used spatially constant RyR densities and membrane radii throughout the study, since we were interested in obtaining functions expressing their relationships with respect to the potential for the elicitation of stable calcium waves. Of course, in real dendrites neither dendrite radius nor RyR density are spatially constant [64], nor can the ER be expected to be a perfectly shaped cylinder that is centrally positioned within the dendrite. Future research, which is outside the scope of this study, could integrate heterogeneous membrane composition and detailed three-dimensional morphologies.

In the context of spike-timing-dependent plasticity, calcium wave velocity may be a controlling parameter. Our study shows that calcium wave velocities vary depending on ER size and RyR density. This merits a closer investigation of the timing of calcium signals in connection to, e.g., back-propagating action potentials.

A more precise understanding of the intracellular architecture may be critical in wave timing and stability. For example, mitochondria are an important sink for excess cytosolic calcium. They can quickly absorb considerable amounts of calcium through mitochondrial calcium uniporter channels [65] and extend through dendrites as filamentous networks [41, 42]. It would therefore be interesting for future work to extend this study to a fully three-dimensional, non-symmetric case, to investigate the interplay between mitochondria and ER in the context of dendritic calcium waves.

In addition to the stability and velocity of dendritic calcium waves, changes in the cellular and intracellular architecture may promote direction selectivity of calcium signals. Since the radius of dendrites is variable, typically large in proximity of the soma and smaller in more distal regions, we wondered whether there could be some kind of direction selectivity in the propagation of calcium waves. Due to the increasing amount of releasable calcium in the growing ER, we suspected that stable waves might be better supported when propagating from thin dendrites to larger ones with larger ER than in the direction of decreasing radii. We tested this hypothesis with three simulations (data not shown), the first on a medium-size dendrite with constant radius, the second on the same dendrite with increasing radius, and the last one on the same dendrite with decreasing radius. In all cases, the ER radius at each location was chosen to be just below the threshold for stable calcium waves at the chosen (constant) RyR density, so that one would expect the waves to terminate at some point. This indeed happened in all cases, however, the distance the waves were able to travel differed: shortest in the shrinking dendrite and longest in the expanding dendrite. The results confirmed our hypothesis. Yet, the effect was relatively small (approx. 11 μm difference in travel distance between small-to-large and large-to-small scenarios) and required very precise calibration of the ER radius (too small and all waves would have terminated even earlier, too large and all waves would have been stable). We therefore believe that this effect is unlikely to be relevant under realistic conditions.

Our study further highlights the question of model dimensionality. While in many cases, in which electrical properties of neurons are studied, one-dimensional multi-compartment models are employed, this may not always be possible for intracellular biochemical modeling. The fact that we could not find a simple descriptive model for regimes with neither very small nor very large distance between ER and plasma membrane reflects the non-trivial dynamics of the calcium signal with respect to space. Even in this relatively simple model setup with rotational symmetry in the geometry, it is not generally possible to reduce dimensionality (by ignoring the radial concentration profile in the cytosol) to a model that is essentially 1D. A fortiori, such a model would not be able to generate the correct stability thresholds and wave velocities in medium-size and large dendrites. For very thin dendrites, however, a simplification of 3D calcium models to 1D is conceivable provided the symmetry prerequisites of this study are met.

Below are the links to the electronic supplementary material.

Video of wave propagation in thin dendrite: Typical wave in a dendrite of radius 0.2 μm, with an ER radius of 50 nm and a ryanodine receptor density in the ER membrane of 3.5 . One second in the video corresponds to 10 ms of simulated time. Calcium concentrations are color-coded blue (, ) to red (, ). Radial coordinate scaled by a factor of 40 to enhance visibility (MKV 74 kB)

Video of wave propagation in thick dendrite: Typical wave in a dendrite of radius 1.0 μm, with an ER radius of 0.2 μm and a ryanodine receptor density in the ER membrane of 2.5 . One second in the video corresponds to 10 ms of simulated time. Calcium concentrations are color-coded blue (, ) to red (, ). Radial coordinate scaled by a factor of 8 to enhance visibility (MKV 88 kB)