Stochastic reaction-diffusion modeling of calcium dynamics in 3D dendritic spines of Purkinje cells.

Calcium (Ca2+) is a second messenger assumed to control changes in synaptic strength in the form of both long-term depression and long-term potentiation at Purkinje cell dendritic spine synapses via inositol trisphosphate (IP3)-induced Ca2+ release. These Ca2+ transients happen in response to stimuli from parallel fibers (PFs) from granule cells and climbing fibers (CFs) from the inferior olivary nucleus. These events occur at low numbers of free Ca2+, requiring stochastic single-particle methods when modeling them. We use the stochastic particle simulation program MCell to simulate Ca2+ transients within a three-dimensional Purkinje cell dendritic spine. The model spine includes the endoplasmic reticulum, several Ca2+ transporters, and endogenous buffer molecules. Our simulations successfully reproduce properties of Ca2+ transients in different dynamical situations. We test two different models of the IP3 receptor (IP3R). The model with nonlinear concentration response of binding of activating Ca2+ reproduces experimental results better than the model with linear response because of the filtering of noise. Our results also suggest that Ca2+-dependent inhibition of the IP3R needs to be slow to reproduce experimental results. Simulations suggest the experimentally observed optimal timing window of CF stimuli arises from the relative timing of CF influx of Ca2+ and IP3 production sensitizing IP3R for Ca2+-induced Ca2+ release. We also model ataxia, a loss of fine motor control assumed to be the result of malfunctioning information transmission at the granule to Purkinje cell synapse, resulting in a decrease or loss of Ca2+ transients. Finally, we propose possible ways of recovering Ca2+ transients under ataxia. Published by Elsevier Inc.

Significance

Ca2+ is a second messenger that can trigger synaptic plasticity in dendritic spines of Purkinje cells, associated among other things with motor learning and motor fine control. Disrupted Ca2+ signals in those spines could lead to pathological conditions such as cerebelar ataxia, a lack of coordination of muscle movements. The molecular reaction mechanisms defining the spatiotemporal aspects of such Ca2+ signals in the noisy environment of dendritic spines in health and disease are not fully understood. Here, we develop a stochastic reaction-diffusion model in MCell to study how synaptic inputs from parallel fibers and climbing fibers reaching spines of Purkinje cells are shaping Ca2+ signals in healthy and pathological conditions and propose a way to recover those Ca2+ signals in pathological conditions.

Introduction

Ca2+ is a second messenger involved in many processes in eukaryotic cells. Ca2+ signals activate several enzymatic targets involved in the induction of synaptic plasticity in dendritic spines of Purkinje cells in the cerebellum and cause an increase or decrease of synaptic strength. Glutamate released at parallel fiber (PF) synapses binds to metabotropic glutamate receptors (mgluRs) on the Purkinje dendritic spines that activate signaling pathways associated with Ca2+ release from endoplasmatic reticulum intracellular stores through inositol trisphosphate receptors (IP3Rs). This response can be further enhanced by a well-timed additional Ca2+ influx triggered by climbing fiber (CF) coactivation (Fig. 1; (6,7)).

Figure 1

Illustration of a spine segment of a Purkinje cell showing the spine head at the top, neck in the middle, and beginning of the dendrite at the bottom. Visible are the signaling pathways of parallel and climbing fiber stimulation (1, 2, 3, 4, 5), which can trigger a cytosolic Ca2+ transient because of an opening of IP3Rs on the ER (6,7). To see this figure in color, go online.

Frequently, the detection and discrimination of transient Ca2+ signals by molecular targets in dendritic spines happen outside mass-action equilibrium, at low particle concentrations of Ca2+ with endogenous buffers involved, and within the rather small three-dimensional structure of the spine which entails a very noisy signaling environment. Therefore, the spatiotemporal characteristics of Ca2+ signals can play an important role in the induction of synaptic plasticity (8, 9, 10, 11, 12), determining whether long-term potentiation (LTP) or depression (LTD) occurs.

Computational models of the three-dimensional structure of dendritic spines and their kinetic mechanisms in their noisy environment can be very helpful in understanding the biophysical constraints shaping the Ca2+ dynamics that control synaptic plasticity in neurons (13). We model the system using stochastic single-particle simulations to account for the noise and low concentration properties of particles involved in the signaling processes (14, 15, 16).

Spines are dynamic extensions of neuronal dendrites and play an important role in cell signaling, neuronal excitability, information processing at the cellular level, and synaptic plasticity (17). They receive synaptic inputs from axons that release neurotransmitters, which bind to postsynaptic receptors on the spines (18,19). Because spines can react to outer and inner stimuli by changes in synaptic efficiency and in their morphological structure, their global topological arrangement becomes a mirror of sensory history and experience. Studying the dynamics of spine behavior is crucial for understanding learning, memory, motor functions, and other large-scale cognitive processes (20, 21, 22, 23, 24).

Spines are targets of signaling and contain molecular signaling mechanisms that regulate and are regulated by intracellular Ca2+ transients. Rapid Ca2+ release is achieved by opening of IP3Rs, which reside on the smooth endoplasmic reticulum (ER), a Ca2+ store inside the spine and dendrite (6,7). IP3Rs require inositol-1,4,5-trisphosphate (IP3) and Ca2+ to open. Both IP3 production and Ca2+ influx are controlled by PF and CF activity (25, 26, 27, 28). These interacting signaling pathways give the IP3Rs the capacity to translate fast PF and CF inputs into longer-lasting slow-output Ca2+ signals (29).

PF stimulation activates two signaling pathways. It triggers glutamate release at the synapse, which leads to an activation of type-1 metabotropic glutamate receptors (mGluR1) located at the postsynaptic density (PSD) at the top of the spine head. The activated mGluRs activate G-protein-coupled receptors (G) that cause the activation of phospholipase Cβ (PLCβ), which synthesizes IP3 from PIP2. IP3 is free to diffuse from the PSD into the dendrite’s cytosol (30, 31, 32) before it vanishes because of degradation by IP3 3-kinase and IP3 5-phosphatase on the timescale of a few seconds (33,34). The second pathway consists of a membrane depolarization, causing Ca2+ influx through voltage-gated Ca2+ channels (Ca2.1 type P/Q-voltage-gated Ca2+ channels), which are highly expressed in Purkinje dendrites (35,36).

CF stimulation also generates a membrane depolarization opening the same type of P/Q-voltage-gated Ca2+ channels, which leads to Ca2+ influx not only into the spine but also into the dendrite (27), summarized in Fig. 1 (1, 2, 3, 4, 5).

Whereas PF synapses are located at the head of a dendritic spine coupled to the PSD, CFs attach to the dendrite itself, creating synapses at the dendritic arbor in ∼2–3 μm intervals (37). One Purkinje cell is connected to up to 105 parallel fibers but to only a single climbing fiber (17).

It is assumed that the activity patterns of PFs attached to Purkinje cells in the cerebellum mediate fine control of movement and promote an increase in synaptic strength (LTP), whereas the activity patterns of CFs encode information about failure of such movement and can trigger a decrease in synaptic strength (LTD) when succeeding PF stimuli (38,39).

The endogenous Ca2+ binding proteins (buffers) calbindin D28k (Cb), parvalbumin (Pv), and calmodulin (CaM) are highly expressed in Purkinje cells (17,40,41). Their role is to shape Ca2+ transients occurring in the cytosol by setting their spatiotemporal parameters such as amplitude and decay time, which are crucial for successful information transmission on cellular level via signaling pathways (42, 43, 44, 45, 46).

Various brain disorders are associated with malfunctioning neuronal information processing which can be related to atypically functioning dendritic spines and IP3Rs (46, 47, 48, 49, 50). Among them, cerebellar ataxia is a serious and heterogeneous neurological condition involving a loss of coordination of muscle movement (51). Most forms of cerebellar ataxia have no cure to this day (52). Thus, it is important to develop computational models to study the role of dendritic spines with respect to ataxia (53). To model ataxia in our approach, we look at Ca2+ transients under reduced IP3 binding rates of IP3Rs and then suggest ways to restore previous Ca2+ transients.

Methods

Model description

The model studies the Ca2+ response to outer stimuli from active PFs and CFs. Simulation methods are explained in the Supporting materials and methods. Stochastic reaction-diffusion, particle-based simulations use MCell (54, 55, 56), and deterministic simulations use Copasi (57). Whereas MCell offers a biophysically realistic approach to a biological problem, accounting for low particle concentrations, inherent stochasticity, and complex three-dimensional geometries (58), Copasi describes the kinetic reactions in a well-mixed volume efficiently, without the influence of diffusion or complex geometry. In this way, the dynamics of the model can be tested in a computationally fast environment before going into more expensive reaction-diffusion simulations in complex geometries. Copasi was, for instance, used to approximately find concentrations of each species at equilibrium, i.e., the initial state then used for simulations in MCell. We focus on explaining model components here. All parameter values not mentioned explicitly in the text are listed in the Supporting materials and methods.

Geometry

We created the simple three-dimensional dendrite-neck-head volume shown in Fig. 2 A for our MCell simulations. The head consists of a sphere with the top and bottom being cut off. The top area models the PSD, and the bottom of the sphere connects to the spine’s cylindrical neck. A cylindrical dendrite is attached to the other side of the neck. The head volume is V = 0.1 μm3. The total volume dendrite, neck and head, is ∼5 times as large, V = 0.512 μm3 (59).

Figure 2

(A) Model geometry. The endoplasmic reticulum (ER) is visible in the head and neck. V = 0.512 μm3, V = 0.100 μm3, and V = 0.020 μm3. Release sites of Ca2+ and IP3 for PF and Ca2+ for CF activation are marked by dots. Exact sizes of the geometry can be found in Table S4. (B) Interaction scheme of particle species used in the simulations from a cytosolic perspective. SERCAs, leak channels, and IP3Rs are located on the ER membrane, and PMCAs, NCXs, and more leak channels are located on the outer plasma membrane. Ca2+, IP3, and the buffers are free to diffuse in the cytosol, the volume within the plasma membrane, and outside the ER.

We created another volume inside the head and neck to represent the smooth ER. It is a scaled copy of the head and neck, just smaller in size, with V = 0.02 μm3. Although the surface of the ER is home to IP3Rs, ER Ca2+-ATPases (SERCAs), and leak channels, we did not model the inside of it. For Ca2+ release by IP3Rs, we assume a constant Ca2+ efflux (60) on the timescale of interest not affected by ER depletion. We are aware of this being an approximation because ER depletion is suggested by the results of Okubo et al. (61), with intralumenal diffusion as the major flux of replenishment. However, lumenal concentrations, together with the ER size required to account for lumenal diffusion on the timescale of release, would render our study unfeasible because of particle numbers above 2 × 105.

Additionally, we used no-flux boundary conditions at the dendrite sites for all particle species, modeling a situation in which neighboring spines also receive stimuli (see (62), their Fig. 1 b).

Molecular components of Ca2+ dynamics

General and steady-state particle concentrations, number of particles, and diffusion coefficients can be found in Tables S5 and S6.

Ca2+ transporters

SERCA pumps are described by a three-state model (63,64), Fig. S3 and Table S4, subsequently binding two Ca2+ before decaying back into the rest state, removing two Ca2+ from the cytosol. We put 68 SERCAs onto the ER membrane (64).

We included five sodium-calcium exchangers (NCXs) (63,64) on the outer plasma membrane without modeling sodium dynamics. Our model assumes constant intracellular and extracellular [Na+] as a simplification. We used a simple two-state model (64), i.e., one NCX receptor can bind one Ca2+ and then either release it back into the cytosol or decay back into the rest state, removing one Ca2+.

Plasma membrane Ca2+-ATPase (PMCA) is another Ca2+ pump that helps to maintain a low Ca2+ concentration in the cytosol of all eukaryotic cells. We used 13 PMCAs (63,64) and also a two-state model similar to the NCX model, but with different reaction rates (65,66).

We include 10 leak channels on each the ER and plasma membrane of the dendrite that yield a small constant influx of into the cytosol. The leak fluxes fix free [Ca2+] 50 nM in steady state.

Buffers. The buffer species in the model are Pv, Cb, and CaM (Fig. 3). We describe Pv by a three-state model. It binds either one Ca2+ or one magnesium (Mg2+) (42,67).

Figure 3

Overview of buffer models for (A) Pv, (B) Cb, and (C) CaM. CaM can hold up to four Ca2+. We used a 16-state model with individual binding sites. For ease of reading, we omitted the six states with two bound Ca2+. Reaction rates for all can be found in Table S3.

Cb is a major Ca2+ buffer in spines of hippocampal and also Purkinje neurons. We used a three-state kinetic model with a low- and a high-affinity site (42,68), with the low-affinity site being present only if the high-affinity site is filled. One Cb protein can therefore bind up to two Ca2+.

CaM is a multifunctional calcium-binding protein expressed in all eukaryotic cells. It is not only a buffer but also acts as a messenger target of Ca2+. Once Ca2+ is bound to CaM, it can modify the interactions of Ca2+ with various other targets like phosphatases or kinases. One CaM protein has four Ca2+ binding sites. We modeled the interactions with Ca2+ using a 16-state model with various medium- and high-affinity binding sites and various reaction rates (66).

IP3R models

We used 54 or 56 IP3Rs (69) on the ER membrane, depending on the parameter set. A large variety of IP3R models have been developed in the last three decades (70, 71, 72, 73). We compare here two models. Doi’s model has been chosen because it has been applied to spine dynamics before (63). We picked Moraru’s model for comparison because its open probability depends nonlinearly on [Ca2+], and it is eligible for easy implementation in MCell (74). This choice allows for statements on the role of IP3R inhibition and IP3R-related slow timescales.

Doi’s model. A seven-state model for the IP3R was proposed by Doi et al. (63). The model is able to reproduce the bell-shaped open probability curve P0([Ca2+]) (Fig. 4 C). An IP3R requires one IP3 to bind before one Ca2+ binds to open; otherwise, it will stay closed and can only bind more Ca2+ if it does not previously release the bound Ca2+. Once the receptor opens, it releases Ca2+ into the cytosol with a constant Ca2+ release rate (Figs. 4 A and S2; Table S2).

Figure 4

(A) Doi’s model: T00 represents an empty IP3R. The IP3R acts as a coincidence detector and only opens if IP3 binds before Ca2+, making T11 the open state. Otherwise, it buffers Ca2+ and becomes inactivated, preventing binding of IP3. (B) Subsection of Moraru’s model. All vertical and horizontal transitions are possible. When four IP3 and two Ca2+ (activating) are bound, T42 (light gray), there is a probability to go into the open state T, whence additional Ca2+ will be released. Further Ca2+ binding of the IP3R will lower the probability of opening, effectively promoting the 10 states T and T to inhibitory Ca2+ states (gray), x ∈ [0, 4]. We set the Ca2+ binding and dissociation rates of the inhibitory states to be slower than the rates of the activating states, expressed by ratio r. Parameter values can be found in Table S1. (C) IP3R open probability at constant [IP3] = 10 μM for Doi’s (dotted) and our version of Moraru’s model with inhibitory Ca2+ binding scaling r = {10−1, 10−2} (dashed and solid, respectively). Data points are results of stochastic computations with Copasi. The rise of the open probability at low Ca2+ causes CICR.

Moraru’s model. We chose a subsection of the IP3R model of Moraru et al. (74), in which we ignored an additional slow IP3 binding site that is irrelevant for the timescales considered in our model. An IP3R can maximally bind four Ca2+ and four IP3 molecules, and any state transition that includes binding or unbinding of one Ca2+ or one IP3 is allowed. There is an additional transition from the state T42 with four bound IP3 (75) and two bound Ca2+ to the open state, T. Moraru’s model also reproduces the bell-shaped open probability curve (Fig. 4 C). The amplitude of the open probability is easily controlled by changing the reaction rates k and k of the transitions between T42 and T (Fig. S1). Because it is known that the binding dynamics of inhibitory Ca2+ is slower than that of excitatory Ca2+ by a factor of up to 100 (70,76), we decreased the inhibitory Ca2+ binding reaction rates into and from states T and T, which represent the inhibitory Ca2+ binding sites in this model, by an additional factor r (see Fig. 4 B; Table S1).

Parallel and climbing fiber stimulation

We studied the Ca2+ response to four different types of stimuli after experimental results (8): a single PF stimulus, PF burst, CF stimulus, and PF burst + CF stimulus. Stimulation by active PFs was simulated by plasma membrane influx of Ca2+ and IP3 production close to the PSD. CF stimulation was modeled by Ca2+ influx close to the PSD and into the dendrite end close to the spine neck (Fig. 2).

We chose different amounts of Ca2+ per PF and CF stimuli (0–1500 Ca2+ and 0–2000 Ca3+, respectively) as part of our parameter scan. A single PF stimulus consists of one instantaneous injection of Ca2+ and production of IP3 (380 close to the PSD), and a PF burst was made up of five single Ca2+ injections at 100 Hz and 1400 caged IP3, yielding an ∼5 times larger IP3 transient. In the latter, more relevant case, free [IP3] usually peaks around 4.5–5 μM. A CF stimulus included an additional Ca2+ release of 200 particles in the dendrite (Fig. 1). The typical time of the CF stimulus is t = 100 ms after the initiation of the PF stimuli but was varied in Optimal timing of CF stimulus.

IP3 dynamics used constant production and decay rates and was able to capture biexponential IP3 concentration behavior (33,34,63,77). We chose an amount of IP3 in agreement with physiological concentrations (6) such that the IP3Rs were saturated with IP3 for the case of a PF burst. IP3 production was delayed by 100 ms compared to the onset of the PF Ca2+ influx to account for the slower process of IP3 synthesis compared to instantaneous Ca2+ influx from PF and CF stimuli. IP3 diffuses freely in our model.

Results

We successfully constructed a three-dimensional stochastic reaction-diffusion model of Purkinje cell dendritic spine Ca2+ dynamics that reproduces many aspects found in experiments (8,62). Because of the nature of computational modeling, we were able to shed light on some aspects of the system’s response to stimuli that are otherwise extremely hard to control experimentally, e.g., removing certain buffer species or changing the amplitude of Ca2+ associated with a PF or CF stimulus.

Snapshots of the spine head including Ca2+, IP3, and IP3R states on the ER from a typical simulation are shown in Fig. 5. At t = 0 ms in the first frame, the red dot is the initially localized collection of 110 Ca2+ of the first PF stimulus. The particles diffuse and get absorbed by buffers immediately. The CF stimulus consists of 1700 Ca2+, which is visible in three frames corresponding to t = 100 ms to t = 100.032 ms, showing how quickly Ca2+ diffuses. IP3 slowly enters the system at the same time at t = 100 ms (see also Fig. 9 A). IP3Rs start to react to IP3 and increased [Ca2+] (note changing colors of IP3Rs). Eventually, a global Ca2+ transient is initiated, which leads to a prolonged increase of [Ca2+], shown from t = 350–650 ms (see also Fig. 6). Using many of these simulations, we studied the IP3-induced Ca2+ responses of our model spine to different stimuli. We chose these stimuli in accordance with available experimental data (8,62). We were able to approximately reproduce properties of Ca2+ transients in the spine head in response to a PF, CF, PF burst, and PF burst + CF stimulus when using Moraru’s IP3R model, as we will see below.

Figure 5

Snapshots of the spine head from a typical simulation with Ca2+ (red), IP3 (beige), and IP3Rs colored according to their state (Moraru’s model) as explained by the legend. At t = 0 ms, the large red dot shows the 110 Ca2+ coming from the first PF stimulus, and at t ≈ 100 ms, 1700 Ca2+ from the CF stimulus. The highly localized Ca2+ stimuli very rapidly spread by diffusion and are absorbed by buffers. Spatial averages of the Ca2+ concentration are shown in Fig. 6. To see this figure in color, go online.

Figure 9

(A) Time courses of average [IP3](t) and SDs in the spine head for different IP3 decay rates are shown. (B) Data points show averages of Ca2+ peak values and SDs for different values of IP3 decay rates kdecay against the binding rate of IP3. Each color and fit represent one decay rate. Binding rates smaller than the original value possibly represent ataxia. (C) Averages of 12 Ca2+ transients and SDs for cases with modified IP3 decay and binding rates, taken from (B), that approximately peak around the original peak value of 360 Ca2+ (dotted red line in (B)) are shown. To see this figure in color, go online.

Figure 6

(A) Moraru’s model: peak values of Ca2+ transients in the spine head for different parameter sets to PF burst and PF burst + CF (PFb + CF) stimuli, averaged over 12 simulations. Error bars show SDs. Parameters are A, N = 54 IP3Rs, PF = 90 Ca2+, and CF = 1700 Ca2+; B, N = 54 IP3Rs, PF = 150 Ca2+, and CF = 1700 Ca2+; and C, N = 56 IP3Rs, PF = 110 Ca2+, and CF = 1700 Ca2+. Set C is our control parameter set and will be used from now on. (B) Ca2+ in spine head after PF burst (dotted) and PF burst + CF (solid) stimuli with control set. The inset shows a magnification of the first 220 ms. At t = 0.1 s, the highly localized Ca2+ from the CF stimulus is visible as a spike. The injected Ca2+ gets absorbed by buffers immediately explaining the immediate return to lower [Ca2+] shortly after the stimulus. The averages with SDs (greyish areas) of 12 simulations are shown.

We focused on the difference of Ca2+ transients, especially on the peak values, upon a PF burst and a PF burst with a CF stimulus coactivation at t = 100 ms after the onset of the PF burst, as these two cases are assumed to encode the induction of LTP and LTD, respectively (38,39,45). We expect the system to show a clear Ca2+ transient with a spine head peak of ∼2.8 μM under a PF burst stimulus and a 150% increase to ∼7.1 μM with a CF coactivation (8), showing a supralinear response to summation of stimuli, also generally found in other model approaches (29,63,78) and experiments (44,62).

Robustness of IP3R dynamics against Ca2+ concentration noise: a model with linear activation characteristics

Local concentration fluctuations at IP3Rs upon opening or closing are large (79,80), and therefore, channel state noise strongly affects channel state dynamics. In this section, we investigate the noise response of linear Ca2+-dependent channel activation in the IP3-sensitized state as, e.g., Doi’s model uses. We and others (63) were able to reproduce dendritic spine Ca2+ dynamics in well-mixed conditions with Doi’s model as a system of ordinary differential equations (ODEs), generating proper Ca2+ transients to different PF and CF stimuli conditions (see also Figs. S6 and S7).

We found with MCell simulations that the stochastic fluctuations of Ca2+ in the cytosol prevent any possible rise of the Ca2+ transient peak with CF or PF stimuli with linear Ca2+-dependent IP3R channel activation. It takes only one Ca2+ to bind to the receptor to open if sufficient IP3 is already present. Once one or two IP3Rs are in the open state just because of basal Ca2+ fluctuations, they release enough Ca2+ to open more IP3Rs to create a global Ca2+ transient with a peak ∼ 3.0 μM (Fig. S5 A). Adding more Ca2+ because of PF (Fig. S5, B–D) or CF (Fig. S5 E) stimuli did not show any further peak increase because the transients arising from basal fluctuations are saturated already. Even large CF Ca2+ amplitudes left the Ca2+ peak values essentially unchanged. The peaks of the Ca2+ transients are essentially constant for Ca2+ PF amplitudes 0–220 (Fig. S5 F). At large PF amplitudes (>220 Ca2+), the inhibitory action of Ca2+ on the IP3Rs decreased transient amplitudes significantly, as some of the total available IP3Rs bind inhibitory Ca2+ before a global transient can be initiated. Results from adding a CF stimulus with increasing CF Ca2+ amplitude are shown in Fig. S5 F, in which a small peak decrease is visible.

This high sensitivity to noise of IP3Rs in this model has not been observed in stochastic simulations based on molecule numbers only, i.e., in non-spatially-resolved simulations (81). Given the same molecule number amplitude, the local concentration amplitude of fluctuations in our spatially detailed simulations is larger than in the spatially lumped simulations of Koumura et al. (81). This effect of local noise most likely explains the different results with respect to noise sensitivity (14,15,82) and renders spatially resolved simulations necessary (83).

Therefore, we turn to a model with nonlinear Ca2+-dependent activation characteristics in the following.

Moraru’s IP3R model

A model with nonlinear Ca2+-activation characteristics like Moraru’s model exhibited better robustness against basal fluctuations, and we use it from now on.

Ca2+ transient peak response to a PF and CF stimulus

Piochon et al. (8) estimated the peak of the Ca2+ transient in the spine head after a single PF stimulus paired with a CF stimulus to be ∼0.4 μM, whereas the response to a single PF stimulus was lost in noise. More interestingly, a PF burst stimulus triggered a Ca2+ response with a peak value of ∼2.8 μM, and a peak value of ∼7.1 μM was reached for a PF burst stimulus with CF coactivation, an increase of ∼150%. Ca2+ peak increase with CF coactivation is crucial for the current understanding of initiation of synaptic plasticity in the form of long-term depression (LTD) (9,17), even though LTD has also been observed after very strong PF stimulation alone (8,26,84,85).

We were able to reproduce Ca2+ transients with peak values in agreement with experimental data for the cases of single PF with additional CF coactivation, PF burst, and PF burst with CF coactivation. Summarized Ca2+ results from our simulations for some example parameter sets are shown in Fig. 6 A, where the peak of the Ca2+ transients in response to a PF burst and a PF burst + CF stimulus are shown. The peak values of Ca2+ transients computed deterministically in Copasi increased clearly with increasing PF and CF Ca2+ amplitudes and showed no saturation for tested parameters. The system was very sensitive to CF coactivation (see Figs. S8 and S9).

Averages and the standard deviation (SD) of actual transients of Ca2+ are shown in Fig. 6 B. The SD due to the inherent randomness is large but does not blur the difference between a single PF burst and combined PF burst + CF stimulus.

We simulated a single PF stimulus with CF coactivation and obtained an average peak value of ∼20 Ca2+ = 0.40 μM from 12 simulations, in agreement with Piochon et al. (8) (Fig. S12).

Closing of IP3Rs was caused by a mixture of reaching the inhibitory states T and T with three or four Ca2+ bound for larger values of r due to increasing [Ca2+] during a Ca2+ transient (see bell-shaped open probability curve, Fig. 4) and IP3 becoming less available during IP3 degradation; see IP3R state occupation videos (Videos S1, S2, S3, and S4) with different values of r in the Supporting material. During a typical Ca2+ transient, [Ca2+] stays approximately constant in the spine head, decreases linearly down the neck and becomes constant again in the dendrite segment of the volume, see Video S5.

Results of Ca2+ transients to different stimuli conditions and associated total Ca2+ release from IP3Rs with single buffer species removed, showing how single buffer species shape Ca2+ transients and influence IP3R dynamics, are presented in Fig. S11.

Ca2+ transient peak scaling with PF and CF Ca2+ amplitude

We find that the relation between the peak value of Ca2+ transients and the CF Ca2+ amplitude is strongly affected by the rate of Ca2+-dependent inhibition of the IP3Rs (Table S1). The rate value suggested in the original Moraru model entails saturation at a CF Ca2+ amplitude of 500 already (diamonds and fit with r = 1.0 in Fig. 7 A). Peak values increase with CF amplitude over a large range with a smaller rate of inhibitory Ca2+ binding to the IP3R, r. We find measured responses of the transient peak to a CF stimulus coactivation and beyond with even larger CF amplitudes with r = 0.01 (circles, Fig. 7 A). We used this value throughout the study. The same applies to the response to PF stimuli with various Ca2+ amplitudes, showing larger peak values for smaller r (Fig. 7 B). Additionally, increasing PF Ca2+ amplitudes (without CF coactivation), mimicking a situation of intense PF stimulation, increases peak values even further, reaching the same or even higher peak values than with CF coactivation and smaller PF Ca2+ amplitudes in Fig. 7 B. This resembles the situation in which a strong PF stimulus alone, rather than PF stimulation with CF coactivation, can trigger LTD (8,26,84,85).

Figure 7

Peak scaling of Ca2+ transients against (A) CF Ca2+ amplitude for the PF burst + CF scenario and (B) against PF Ca2+ amplitude for the PF burst scenario for three different inhibitory Ca2+ binding timescale ratios r = {1.0, 0.1, 0.01}; see Table S1. The slower the rate of inhibition, the larger the resulting Ca2+ transients for a given CF or PF Ca2+ amplitude. The value of r sets the largest possible peak value when varying PF and CF Ca2+ amplitudes. Large enough transients triggered by CF coactivation with intermediate Ca2+ amplitudes are only reached when r is on the order of 10−2. Data points are averages of 12 stochastic simulations and were fitted to Hill curves, error bars show SDs.

Optimal timing of CF stimulus

The size of the Ca2+ transient elicited by the CF stimulus, and with it the induction of LTP and LTD, responds optimally to a certain timing of the CF stimulus relative to the PF stimulus as Wang et al. have shown (62). They measured the Ca2+ response to different timing windows between CF and PF stimulus and used a Gaussian to fit the Ca2+ transient’s total integrated response to their results, which peaked around 92 ± 37 ms and had a halfwidth of 212 ± 85 ms.

We simulated these experiments by varying the CF stimulus time from 0 to 400 ms after the initiation of the PF stimulus. We find that the system is sensitive to the timing of the CF stimulus, as shown when determining the Ca2+ transient peaks under such CF Ca2+ timing variation and also exhibits optimal time windows with a maximal peak (Fig. 8). We find optimal responses for different strengths of PF stimuli also including the parameter value set A in Fig. 6.

Figure 8

(A) Peak of Ca2+ transients against timing of CF stimulus for three parameter sets from Fig. 6 and an additional one with r = 0.03. 12 stochastic simulations per data point were used, and points were fitted to Gaussians with the condition to converge against peak values of PF burst stimulation alone for very large and small t. Error Bars indicate SDs. The maxima of the Gaussian fits occur at 91 ms (red), 59 ms (blue), 16 ms (green), and 80 ms (purple), in qualitative agreement with experimental data (62) (+92 ± 37 ms). The halfwidths of the Gaussian fits are 257 ms (red), 294 ms (blue), 443 ms (green), and 260 ms (purple), in qualitative agreement with experimental data (62) (212 ± 85 ms). (B) Illustrating PF and CF stimulus timing by showing five Ca2+ spikes in spine head from a PF burst stimulus at 100 Hz from 0 to 40 ms and a Ca2+ spike from a CF stimulus, here t = 100 ms. To see this figure in color, go online.

We compare our results in Fig. S10 to the experimental data from Wang et al. (62).

The rising phase of the Ca2+ peaks in Fig. 8 is due to Ca2+-induced Ca2 release (CICR) (see Fig. 4 C). The CF stimulus causes Ca2+ influx. It takes more than 200 ms for this Ca2+ rise to decay back close to prestimulus levels (Figs. 6 and S12). When IP3 production starts 100 ms after onset of PF stimulation, CICR starts because of the presence of IP3 and the remaining Ca2+ from the CF stimulus. The closer to IP3 production the CF stimulus occurs, the stronger the CICR. Interestingly, this does not necessarily lead to an optimal response at a timing window at the time of onset of IP3 production at 100 ms, as the blue and green curves in Fig. 6 show. We did not observe optimal time windows when we released IP3 at the onset of PF stimulation (data not shown).

The decaying phase of the Ca2+ peaks in Fig. 8 toward large time windows is affected by processes terminating Ca2+ release. One of them is the decay rate of IP3 as the simulations in Fig. 9 show (which we will discuss in more detail below). The slower IP3 is removed from the spine head—either by degradation or by diffusing out of the spine—the longer the Ca2+ transient.

Ca2+-dependent inhibition is another process contributing to the termination of Ca2+ release and affecting the peak dependency on the CF time window as the comparison between the purple curve and all other ones in Fig. 8 shows. Additionally, the peak scaling data in Fig. 7 illustrate the role of Ca2+-dependent inhibition in setting peak height. Increasing the rate of Ca2+-dependent inhibition leads to a longer time window providing optimal response (purple curve in comparison to the blue one).

However, increasing the inhibition rate by a factor of three does not shorten the Ca2+ transient by the same factor and does not abolish optimal time windows because Ca2+-dependent inhibition is only one of several factors shaping the transient.

Ataxia

It has been shown that spinocerebellar ataxia type 29 (SCA29), characterized by early-onset motor delay, hypotonia, and gait ataxia, can be caused by malfunctioning type 1 IP3Rs (50). Mutations associated with SCA29 were identified within or near the IP3-binding domain. These mutations interfere with the binding of IP3 and cause IP3Rs of type 1 to lose any channel activity, reducing or removing IP3-induced Ca2+ transients.

On that basis, we decrease the IP3 binding rate k of the IP3R model to mimic ataxia, which results in lower or vanishing Ca2+ transients. We search to rescue the system from this pathological condition by trying to recover the original Ca2+ peak in two ways. First, we increase the amount of IP3 that enters the system, representing increased activity of the PLCβ pathway, which synthesizes IP3. In a second approach, we decrease the degradation rate of IP3, thus increasing the IP3 that is available to the IP3Rs in absolute number as well as in duration.

The first method of increasing IP3 was only able to recover the Ca2+ transients if we increased the amount of IP3 like 1/k (see Fig. S13). Because the decrease of k might be substantial (50), the [IP3] values compensating it are likely beyond the saturation values of the PLCβ pathway.

Prolonging IP3 exposure

We decrease the IP3R’s binding rate of IP3 to values possibly representing ataxia and then reduce kdecay trying to recover the original Ca2+ transient. We start from their standard values kon = 83.3 (μM s)−1 and kdecay = 15 s−1.

The control Ca2+ peak value can be recovered because slower IP3 degradation increases the amount and duration of IP3 in the system (Fig. 9 A), making up for the negative effects of slower IP3 binding. Additionally, the decrease of kdecay leads to prolonged activity of the open IP3R. The slowed IP3R dynamics also cause some delay in reaching the Ca2+ peak (Fig. 9 C). Whereas the control parameters yield a Ca2+ transient peak at ∼0.47 s (red, Fig. 9 C), slowing IP3 degradation down to one-sixth kdecay = 2.5 s−1 delays the peak to ∼0.65 s (orange), i.e., it increases the response time by ∼40% and increases the width of the Ca2+ transient.

We provide a more systematic analysis in Fig. 9 B. It shows the peak values in dependence on k for five different IP3 decay rates kdecay. The red curve shows results with the control value of k. The curves with reduced kdecay cross the Ca2+ peak control value 360 (red dotted line) at specific values k, which are smaller than the k control value. They are related to kdecay approximately by kdecay ≈ . Simulations with the parameter value pairs of (kdecay, k) calculated according to this equation provide control of Ca2+ peak values with our control parameter set for all other parameter values. A decay rate reduction calculated according to this equation compensates the pathological reduction of k with respect to the Ca2+ transient peak.

Discussion

Cerebellar learning theories suggest that learning is expressed as a change of neuronal weights, i.e., synaptic strengths, reflecting the topological properties of a neuronal network state. Understanding learning therefore requires knowledge of the molecular mechanisms assumed to encode synaptic plasticity and information transmission at the lowest neuronal level, which are Ca2+ transients and the associated cell responses in synapses of spines, eventually. In Purkinje neurons, the IP3-induced Ca2+ transients are dynamical responses to outer stimuli from PFs or CFs happening at low Ca2+ concentrations. Whereas PFs are assumed to carry information about movement and fine motor control, CFs are assumed to carry error information that gives feedback about the network state that triggered the movement (86, 87, 88). We developed a model that is based on complex single-particle stochastic reaction and diffusion processes within a small three-dimensional geometry to study Ca2+ transients in response to dynamical PF and CF stimuli.

Our use of three-dimensional stochastic simulations demonstrated the necessity for IP3R models to filter out Ca2+ binding noise to a sufficient degree. A linear relation between [Ca2+] and the open probability at small concentrations appears not to provide that filtering and entailed Ca2+ dynamics insensitive to CF and PF stimulus Ca2+ amplitudes. However, an increase of the Ca2+ transient peak due to a CF stimulus provides meaning to this stimulus and is thus a necessary model requirement.

Using Moraru’s IP3R model to provide sufficient noise filtering, we were able to reproduce the dynamic behavior of the Ca2+ transients from experiment with respect to the absolute and relative peak values of Ca2+ transients under stimuli (8) and the behavior of peaks under variation of the timing of the CF stimulus (62).

The signal of the CF stimulus turning LTP into LTD might be binary information or graded information. If simply the presence of a stimulus entails LTD, we face binary signaling. If the strength of the stimulus encodes the strength of depression, we see a graded response. We found this characteristic of the signaling by the CF stimulus to depend on the rate of Ca2+-dependent inhibition of the IP3R. We achieved agreement of Ca2+ transient peak values with experimental results at slow inhibition rates. Although these rates are slower than originally suggested by the authors of the model, they are still compatible with puff data of the IP3R taking the large local [Ca2+] at puff sites into account (79,89,90). In summary, these simulation results suggest a graded response of the Ca2+ transients’ peak value to the CF Ca2+ amplitude.

Using our model also allowed for detailed tests on the effects of endogenous Ca2+ buffer molecules. We find clear indication that buffers do not only passively shape amplitude and decay and rise times of Ca2+ transients but actively modulate state dynamics of IP3Rs, resulting in an increase or decrease of released Ca2+.

It has been experimentally established that the timing window of PF and CF stimuli is critical to the induction of LTP and LTD and thereby also to the properties of the IP3R-induced Ca2+ transients (62). The mechanism our results suggest is that the relative timing of IP3 production elicited by PF stimuli, which sensitizes IP3Rs for CICR to the moment of the influx of Ca2+ due to the CF stimulus causing CICR, sets the optimal time window.

We simulated ataxia on the basis of the assumption that it manifests itself in our model by a substantially reduced rate of binding of IP3 to the IP3R. Reducing the rate of IP3 degradation by IP3 3-kinase and IP3 5-phosphatase turned out to be able to compensate the reduced binding without strong stimulation of the PLC pathway. Whether the recovered peak values are enough to trigger an increase of AMPA receptors in agreement with observations by Piochon et al. (8), even with delayed peak times and decay of the Ca2+ transients, is an open question left for future research. Its outcome decides whether reduction of IP3 degradation offers new ways of addressing ataxia pharmacologically.

Author contributions

G.A. and F.M.S.d.S. designed research. V.N.F. performed research, contributed analytic tools. M.F. supervised simulations. V.N.F., M.F., and F.M.S.d.S. analyzed data and wrote the manuscript.