This Article Abstract Full Text (PDF) Supplemental All Versions of this Article: biophysj.105.065771v1 89/6/3790    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Hernjak, N. Articles by Loew, L. M. Search for Related Content PubMed PubMed Citation Articles by Hernjak, N. Articles by Loew, L. M.

Modeling and Analysis of Calcium Signaling Events Leading to Long-Term Depression in Cerebellar Purkinje Cells

Center for Cell Analysis and Modeling, University of Connecticut Health Center, Farmington, Connecticut

Correspondence: Address reprint requests to Leslie M. Loew, Tel.: 860-679-3568; E-mail: les{at}volt.uchc.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Modeling and simulation of the calcium signaling events that precede long-term depression of synaptic activity in cerebellar Purkinje cells are performed using the Virtual Cell biological modeling framework. It is found that the unusually high density and low sensitivity of inositol-1,4,5-trisphosphate receptors (IP₃R) are critical to the ability of the cell to generate and localize a calcium spike in a single dendritic spine. The results also demonstrate the model's capability to simulate the supralinear calcium spike observed experimentally during coincident activation of the parallel and climbing fibers. The sensitivity of the calcium spikes to certain biological and geometrical effects is investigated as well as the mechanisms that underlie the cell's ability to generate the supralinear spike. The sensitivity of calcium release rates from the IP₃R to calcium concentrations, as well as IP₃ concentrations, allows the calcium spike to form. The diffusion barrier caused by the small radius of the spine neck is shown to be important, as a threshold radius is observed above which a spike cannot be formed. Additionally, the calcium buffer capacity and diffusion rates from the spine are found to be important parameters in shaping the calcium spike.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
One cellular basis for learning is the phenomenon of synaptic plasticity that has been observed experimentally in neurons. Synaptic plasticity refers to a temporary or sustained depression or potentiation of the activity level of a synapse. An important form of synaptic plasticity related to motor-learning tasks such as the vestibular-ocular reflex, eye-blink conditioning, posture and locomotion adaptation, motor coordination, and hand/arm movement adaptation (1) is observed in cerebellar Purkinje cells (see Fig. 1 for a micrograph of a Purkinje cell). This particular form of synaptic plasticity, known as long-term depression (LTD), is a lasting decrease in the activity of the synapses between spines on the Purkinje cell dendrites and axons of neighboring granule cells, often referred to as parallel fibers (PF).

View larger version (85K): [in this window] [in a new window]   FIGURE 1  Micrograph of a Purkinje neuron (left) and, at high magnification, a dendritic branchlet with four spines (right). (These images were taken from Palay and Chan-Palay, 1974 (59), with kind permission of Springer Science and Business Media.) The complex structure of the dendritic arbor is visible in the figure on the left. The figure on the right shows the attachment of the spines to the dendritic branchlet via the thin spine necks.

View larger version (36K): [in this window] [in a new window]   FIGURE 2  Cytosolic calcium handling after IP3-dependent release from the endoplasmic reticulum. The cascade of events is triggered when an agonist (glutamate (Glu)) binds to its receptor (mGluR) in the plasma membrane. This sets off a G-protein cascade (the G-protein complex is shown as Gq,ß) that activates phospholipase C (PLC) that, in turn, hydrolyzes the glycerolphosphate bond in phosphatidylinositol-4,5-bisphosphate (PIP2). One product of this hydrolysis goes on to activate protein kinase C (PKC), whereas the other, IP3, is released from the membrane. The IP3 is free to diffuse through the cytosol, be degraded by phosphatases and kinases, and bind to its receptor (IP3R) in the endoplasmic reticulum (ER) membrane. The IP3R is a calcium channel that is triggered to open when IP3 is bound and when calcium itself binds to an activation site; a slower inactivation process also pertains to this channel. The calcium that is thus released binds to calcium buffers (B) in the cytosol including any fluorescent calcium indicators. Calcium influx through plasma membrane channels (CaCh) and release from the ER through the ryanodine receptor calcium channels (RyR) may also contribute to the cytosolic calcium changes. Note that RyR are not present in the spines of Purkinje cells but can be found in the dendrites. Finally, calcium is pumped back into the ER via a calcium ATPase (SERCA) and can also be extruded from the cell by a number of exchange and pump mechanisms (Ca ext).

It has been found experimentally that coincident activation of the PF and CF inputs results in a supralinear increase in [Ca²⁺]. In other words, the change in [Ca²⁺] that is observed is significantly more than the sum of the Ca²⁺ responses obtained by exciting the PF and CF separately (3). It is hypothesized that this supralinear calcium response is the mechanism by which the cell detects the coincident activation of the PF and CF and is the first step in the mechanism leading to LTD. Under normal coincident activation conditions, these supralinear spikes may be confined to single spines.

The properties of the IP₃R involved in this Ca²⁺ signaling pathway are likely to play a key role in the generation of the supralinear Ca²⁺ increase observed during coincident PF and CF activation. In addition to IP₃ binding sites, the IP₃R also contain Ca²⁺ binding sites that further activate Ca²⁺ release. This phenomenon is referred to as calcium-induced calcium release (CICR) (this is also observed in the Ryanodine receptors found in the Purkinje dendritic shafts, but not the spines themselves). Therefore, there is a positive feedback present in the system in which Ca²⁺ release stimulates further Ca²⁺ release. Note that there is an apparent threshold Ca²⁺ concentration above which an IP₃R becomes increasingly deactivated but at a slower timescale (9–12).

The IP₃R in Purkinje cells have certain, unique properties as compared to other neuronal or peripheral cells. Namely, the receptor density is significantly higher (13–16) and the sensitivity to IP₃ is extraordinarily lower in vivo (17,18). It has been shown that the sensitivity of the Purkinje cell IP₃R is similar to those of other cells when analyzed in isolation, leading to the hypothesis that the reduced sensitivity is due to interference by an inhibiting protein found in Purkinje cells (19,20). The physiological significance of these properties in terms of Ca²⁺ release before the onset of LTD is unclear at this time. For example, would similar Ca²⁺ release profiles be observed if the cell had a significantly lower IP₃R density but a much higher sensitivity, i.e., the density and sensitivity common in other cells? It is possible that factors beyond those involved in the LTD-induction mechanism mandate these properties. Additionally, as is supported by the results of this work, it is likely that the CICR phenomenon at the IP₃R is the basis of the coincidence detection mechanism that triggers LTD. Therefore, it is important that the unique properties and influence of the IP₃R be correctly modeled and analyzed as part of a study of the Ca²⁺ dynamics that precede LTD.

When considering the origins of Ca²⁺ release patterns in Purkinje spines, the unique geometry of the spines must also be taken into account. It is hypothesized that the thin neck connecting a spine to its dendrite aids in localizing LTD to a particular spine by creating a significant diffusional resistance (21). This could aid in biochemically decoupling the spine from the parent dendrite (22–25), which is desirable since the Ca²⁺ signaling events can thus be made specific to a single synapse. In summary, the available data in the literature implies that the Purkinje cell's ability to produce the supralinear behavior may rely on a variety of biochemical and geometrical effects.

Given the importance of Ca²⁺ signaling to the induction of the mechanisms leading to LTD, the objective of this work is to use mathematical models of a Purkinje cell that focus on the relevant Ca²⁺ signaling networks to investigate the significance of certain unique characteristics of the Purkinje cell, such as the sensitivity and density of the IP₃R and the geometry of the spines, in terms of LTD induction. The results of this work will aid in identifying those features of the cell that are most critical to the onset of LTD, including consideration of both biochemical and geometrical effects. The wide availability of experimental data on calcium dynamics in Purkinje cells makes such a modeling study feasible.

A recent modeling study of Purkinje spine Ca²⁺ signaling (26) contained analysis of a compartmental model of a single spine that did not consider diffusion out of the spine into the adjacent dendritic shaft or the unique sensitivity and density of the IP₃R. These are key differences between the models proposed here and those in the other work. Indeed, it is these features of Ca²⁺ signaling in Purkinje cell spines that set it apart from other systems. As will be shown, without taking into account the unique in vivo sensitivity of the IP₃R, significant Ca²⁺ release would be produced at low [IP₃] in direct disagreement with available experimental data (e.g., (17)). This combination of biochemical and geometric features also assures that the signal is confined to the activated spine. Additionally, as shown experimentally by Noguchi et al. (25) in pyramidal neurons, diffusion through the spine neck is an important factor in spine Ca²⁺ signaling, suggesting that dynamic spine morphology can be a powerful modulator of synaptic plasticity. The modeling results shown here demonstrate that the spine neck diameter can behave as a sensitive switch for the Ca²⁺ signal.

In the next section, the models used in this work are introduced and the sources for all parameter values identified. In the subsequent sections, simulation results are presented for a one-dimensional spatial model, a compartmental model, and a two-dimensional spatial model. These results demonstrate the models' ability to rationalize the high density and low sensitivity of the spine IP₃R, to generate the supralinear Ca²⁺ spike, and to elucidate the mechanisms underlying the supralinear behavior. In the last section, conclusions are drawn and future research directions discussed.

	MODEL DEVELOPMENT

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The models developed in this work are based on an assembly of components from our own previous work and those available in the literature relating to cellular calcium dynamics. Additionally, many parameter values are taken from a variety of pertinent references and are summarized in Table 1. The modeling and simulation are performed using the Virtual Cell (http://vcell.org) biological modeling framework (27,28). All models used in this work are publicly available in the Virtual Cell. To access, view, or copy the models described in this article, log on to the Virtual Cell (free registration is required), and go to File Open BioModel. In the BioModel Database dialog, the two-dimensional and compartmental models can be found in the Shared Models directory under "hernjak." The one-dimensional model is similarly available as a MathModel. In what follows, the features shared by all of the models are introduced and the sources for the parameter values discussed. A concise summary of all of the models is provided in the Appendix.

(1)

Each of the terms on the right-hand side of Eq. 1 will be discussed individually in what follows. Because the methodology for modeling diffusion differs among the models developed in this work, discussion of the different forms of R_diffusion can be found in the corresponding sections below.

The rate of binding of Ca²⁺ to a buffering species (X) is modeled as

(2)

where XB denotes the Ca²⁺-bound form of the buffer. CD28k is found to contain both medium-and high-affinity binding sites (29), accounted for in this model under the assumption that the sites bind independently of each other. The overall R_buffering term in Eq. 1 is the summation of the buffering rates for all of the buffers, each modeled as in Eq. 2. Parvalbumin and calbindin are the only endogenous buffers considered in this work. According to the experimental data of Schmidt et al. (30), the intracellular concentrations of other buffers (e.g., calmodulin) are considerably lower than those of parvalbumin and calbindin and provide small contributions to the overall Ca²⁺ dynamics in the system. Similarly, the role of nonprotein buffering species is also likely to be negligible, given the conditions explored in this work.

The behavior of the IP₃R is modeled using the formulation of Li and Rinzel (31). This model, which is a simplified version of the DeYoung-Keizer (32) model, also accounts for the CICR behavior that the receptor demonstrates. Although a number of other models for calcium release from the ER have been proposed (e.g., the models of (33–40)), the Li-Rinzel model is adequate for this work's analysis as it represents the critical dynamics of the calcium release via the IP₃R in a compact form. The Li-Rinzel model for Ca²⁺ flux across the ER membrane (including the behavior of SERCA pumps at the ER and an allowance for Ca²⁺ leakage from the ER) is defined as

(3)

The variable h is the probability of an inhibition site on the receptor being unoccupied. For the purposes of this work, the most important parameters in the model are a and d_IP3. The a parameter corresponds to the density of IP₃R in the system, with large values of a indicating a high density of IP₃R. The value of d_IP3 controls the sensitivity of the IP₃R to [IP₃], with increasing values of d_IP3 indicating decreased sensitivity of the IP₃R to [IP₃] changes. Physically, d_IP3 is the dissociation constant for IP₃ binding to the channel. The d_Ca parameter performs the same task in determining the sensitivity of the CICR behavior to changes in [Ca²⁺]. Although most of the Li-Rinzel model parameters are taken to be the values obtained from modeling studies of other neurons (such as neuroblastoma cells (41)), the values for a and d_IP3 were both increased to correspond to the increased density and decreased sensitivity of the IP₃R found in Purkinje cells. Implicit in the model described by Eq. 3 is the assumption that the concentration of Ca²⁺ in the ER is constant. Although researchers have adapted the Li-Rinzel model to account for varying [Ca²⁺]_ER (e.g., (42–44)), the assumption of constant [Ca²⁺]_ER is appropriate for this work, since the timescales of the Ca²⁺ dynamics investigated here are very fast relative to the proposed timescales for Ca²⁺ dynamics in the ER and the magnitudes of the Ca²⁺ changes that are observed in the cytosol are more than an order-of-magnitude less than [Ca²⁺]_ER. Also note that the magnitude of the leakage term in Eq. 3 is normally very small relative to the other terms, and likely plays a negligible role; however, it is left in the equation for completeness.

To simulate the effect of depolarization of the Purkinje cell by the CF, expressions of the following type are added to the dendrite and spine compartments:

(4)

(where ₁ < ₂). [Ca²⁺]_ex is the concentration of Ca²⁺ in the extracellular space and the J_ch parameter controls the magnitude of the Ca²⁺ entry rate (see Table 1). The two conditional statements in Eq. 4 are equal to 1 when true and 0 when false. Therefore, this expression allows for a pulse of Ca²⁺ to flow between the extracellular space and the cytosol during the time period ₁ < t < ₂. Experimental data showing [Ca²⁺] changes during CF activation in Purkinje cells (30) were used to find realistic J_ch values for the spine and adjacent dendrite regions.

The production of IP₃ due to activation of the PF is modeled using techniques developed in a previous modeling study of neuroblastoma cells (41). Based on the experimental observations of Fink et al. ((41), Fig. 4), activation of mGluR is modeled as an exponentially decaying flux of IP₃ from the membrane of a spine. Multiple activations produce multiple pulses of IP₃. Activation of a PF does not trigger IP₃ production in the dendrite shaft under normal conditions, but generates IP₃ only in the spine. IP₃ can, however, diffuse into the dendritic shaft through the spine neck, but this geometrical feature will be described in the next section. The IP₃ pulses are modeled as

(5)

where n is the number of pulses, ₃ is the time between pulses, and J_p and K₃ control the pulse magnitudes and decay times, respectively. Based on commonly adopted experimental protocols for induction of LTD (e.g., Wang et al. (3)), n was chosen as 12 with a time between pulses (₃) of 0.012 s, corresponding to a frequency of 80 Hz (10). The second term in Eq. 5 accounts for the degradation of IP₃ to IP₂ and IP₄ where K_deg is the degradation rate and [IP₃]_o is the nominal concentration of IP₃. The value for K_deg was obtained from experimental data included in a modeling study of neuroblastoma cells (41). It is possible that the value for K_deg in Purkinje cells may differ significantly from the neuroblastoma values, but data does not exist to support the use of other values of this parameter.

View larger version (19K): [in this window] [in a new window]   FIGURE 4  Kymograph showing distribution of IP3 in the spines along the dendrite during PF activation of the center spine (top). IP3 profile in the center spine during simulated PF activation (bottom). Data obtained from the one-dimensional model.

	ONE-DIMENSIONAL SPATIAL MODEL

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The first model considered in this work is a one-dimensional spatial model of a dendritic segment of length 500 µm, as shown in Fig. 3. The results discussed in this section will demonstrate, via thorough modeling of spatial effects, that the relevant Ca²⁺ dynamics are localized to the activated spine under normal conditions. The governing equation for diffusion and reaction for each species in the dendrite has the form of the standard reaction-diffusion partial differential equation (PDE),

(6)

where [X]_d represents the concentration of species X in the dendrite, z represents the length dimension, and R_i refers to the rates of the buffering reactions and other pertinent reactions as discussed in the previous section. Species diffusion constants (D_X) were obtained from experimental data (30). The validity of this one-dimensional model depends on the assumption that the dimensions of the spines and the radius of the dendrite are small enough such that no spatial gradients exist in those directions. Thus, only diffusion along the length of the dendrite needs to be explicitly considered. The large length of the geometry was selected to make sure that boundary conditions have no effect on the solution. A standard Fick's law expression is used to account for diffusion from a spine to the adjacent dendritic shaft region through the diffusion barrier posed by the spine neck. Specifically, the change in concentration of species X in a spine due to flux between the spine and a point along the dendrite adjacent to the spine and the various flux rates discussed in the previous section is modeled using an expression such as

(7)

where l_n is the length of the spine neck, A_n is the cross-sectional area of the spine neck V_s is the volume of the spine is the concentration of the species in the proximal dendrite, and [X]_s is the concentration of species X in the spine. Production of IP₃ due to PF activation (as in Eq. 5) is only simulated in the spine at the center of the dendrite so that the Ca²⁺ dynamics resulting from stimulation of a single spine can be studied. This model also accounts for Ca²⁺ extrusion through the plasma membrane using

(8)

where is the surface/volume ratio of either the spine or the dendritic shaft, P is the rate of pumping, and [Ca²⁺]_T is a threshold Ca²⁺ concentration below which the pumps are inactive. Equation 8 was obtained from the work of Fink et al. (41) and is based on the experimental observations of Herrington et al. (46), which indicate that the pumping mechanisms behave linearly at moderate [Ca²⁺] above a threshold value. Because high values of [Ca²⁺] are reached only for very short time-intervals in this work, this approximation is sufficient. This approximation may tend to overestimate pumping rates at some values of [Ca²⁺] since pumps are known to saturate. Likewise, the effect of other species' concentrations (e.g., Mg²⁺) on the capacity of the pumps is ignored for similar reasons. The overall system of equations is solved by the Virtual Cell through use of a finite volume algorithm with computation points evenly distributed at 0.5-µm intervals along the dendrite.

View larger version (12K): [in this window] [in a new window]   FIGURE 3  Visualization of the geometry used to simulate calcium dynamics in a dendritic spine. For the one-dimensional spatial model (see article), the diffusion along the dendrite (J23) is modeled explicitly using the PDE in Eq. 6 and the dynamics in the spines is modeled using Eq. 7. For the compartmental model (see article), the first two of the following three compartments are each modeled by a set of ordinary differential equations: 1, The spine under consideration; 2, The region of the parent dendrite directly adjacent to the spine; and 3, The distal regions of the dendrite. Species concentrations in compartment 3 are held constant so that this compartment acts as a sink for material that enters compartment 2. Compartment 1 is modeled as a sphere with volume where r1 (rs) is the radius of the spine. Compartment 2 is modeled as a cylinder with volume where r2 is the radius of the dendrite and l2 is the length of the segment of the dendrite encompassing compartment 2. The fluxes, Jij, between compartments are modeled using Fick's law as shown in Eq. 7.

View larger version (34K): [in this window] [in a new window]   FIGURE 5  (Top) Kymograph showing the spatial distribution of IP3 in the dendritic shaft during PF activation of the center spine. (Middle) The kymograph from Fig. 4 is shown again but with the color scaling set to correspond to that of the top figure (hence the color saturation for the center spine data) for comparison. (Bottom) Comparison of the [IP3] time courses obtained in the spine and dendrite segment at 5 µm to the right of the center spine. The results show that the neighboring spines reach similar [IP3] as the adjacent sections of the dendrite, meaning that the diffusion barrier caused by the spine neck is not sufficient to block IP3 entry at such slow time scales. Data obtained from the one-dimensional model.

View larger version (18K): [in this window] [in a new window]   FIGURE 6  (Top) Kymograph showing the spatial distribution of Ca2+ in the spines during PF activation of the center spine. (Middle) Kymograph showing the spatial distribution of Ca2+ in the spine during coincident PF and CF activation. Ca2+ spikes rapidly in the center spine and then quickly returns to steady-state values, as shown in the timecourse of [Ca2+] in the center spine obtained during coincident activation of the PF and CF (bottom). The kymographs demonstrate that Ca2+ does not appreciably spread to other spines during PF or coincident PF and CF activation. Data obtained from the one-dimensional model.

View larger version (20K): [in this window] [in a new window]   FIGURE 7  (Top) The percent change in Ca2+ spike magnitude during coincident PF and CF activation as a function of the percent change in the radius of the spine neck. (Bottom) Several Ca2+ timecourses obtained during coincident activation with varying neck radii. As seen in the figure on the left, neck radii increased by 30% or more prevent a substantial Ca2+ spike from being formed. Data obtained from the one-dimensional model.

	COMPARTMENTAL MODEL

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
In this section, simulation results are presented that demonstrate the ability of a compartmental model to generate a supralinear Ca²⁺ spike given physically realistic Ca²⁺ signals originating from the PF and CF inputs. Additional simulation results are then presented that demonstrate the sensitivity of the system to certain biochemical features. Results are also shown that elucidate the mechanisms leading to the supralinear Ca²⁺ spike. In total, the results demonstrate that a wide range of biochemical and geometrical effects act in concert to allow the system to generate the supralinear Ca²⁺ signal that triggers the LTD process.

The model used in this section focuses on events that occur in a single spine and in the region of the parent dendrite directly adjacent to the spine. The modeling approach is summarized in Fig. 3. As seen in the figure, the model treats the region as consisting of three well-mixed compartments: the spine under consideration, the region of the parent dendrite directly adjacent to the spine, and a compartment representing the combined regions found at the distant ends of the dendrite. Material is allowed to diffuse from one region to the next as depicted in the diagram in Fig. 3. Species' concentrations in the distal dendrite regions are held constant, implying that the distal regions are sufficiently far from the proximal dendrite that they act as a sink for material coming from the spine. This approximation is justified by the results of our model of a 500-µm length of dendrite, as described in the previous section.

Flux from one compartment to the next is modeled assuming diffusion occurs according to Fick's law, as represented in Eq. 7. To determine the effective distance from the proximal dendrite to the distal dendrite (l₂₃), the one-dimensional PDE model of diffusion along a dendrite was used to identify the value of the length parameter that yields the best agreement between the compartmental model and the PDE model. A key benefit of the use of the ODE model is the significant reduction in model complexity and computational cost with only a small loss in fidelity. Further details about this model may be found in the Appendix.

In all simulations, 160 µM Calcium Green-1 (CG) is included in the model as an indicator dye to provide better comparison of the results to experimental data. In the analysis, the indicator is treated as a buffering species and, in most cases, free [Ca²⁺] is reported. Using the Virtual Cell modeling framework, the relative fluorescence change, F/F, can also be calculated based on the fluorescence properties of CG versus CG bound to Ca²⁺. In any event, it is important to keep the indicator species in the model because, due to its buffering capability, it may have an effect on the magnitude of the Ca²⁺ signals observed (41).

Inducing the supralinear Ca²⁺ signal
Fig. 8 a shows the Ca²⁺ profiles in the spine and adjacent dendrite compartments during the simulated CF activation with no PF stimulation. These results were obtained by adjusting the J_ch parameters to match the experimental data of Schmidt et al. (30). During the simulated opening of the voltage-sensitive calcium channels, [Ca²⁺] rapidly spikes in both the spine and adjacent dendrite with a larger magnitude spike in the spine. Decay rates are comparable in both compartments. Note that the magnitudes of both spikes are in the submicromolar concentration range.

View larger version (13K): [in this window] [in a new window]   FIGURE 8  (a) Simulated Ca2+ signals in the spine (top) and proximal dendrite (bottom) cytosol during CF activation. Data obtained from the compartmental model. (b) Simulated IP3 (top) and Ca2+ (middle) signals in the spine cytosol during activation of the PF with 12 stimuli at 0.012 s intervals. The percent change in fluorescence (bottom) based on the fraction of bound CG sites is also shown. Data obtained from the compartmental model.

Fig. 9 is a plot of [Ca²⁺] in the spine during coincident activation of the PF and CF inputs. A plot showing the relative timing of the individual Ca²⁺ signals is also included in the figure. This relative timing was found by trial and error to be optimal in the sense that if the timing between inputs strayed too far from the selected timing, the large-magnitude supralinear spike was not observed. The effect of the relative timing of the inputs on the supralinear behavior will be discussed in more detail in the following subsection. The [Ca²⁺] results in Fig. 9 clearly show the ability of the model to demonstrate the critical supralinear behavior of the Purkinje cell Ca²⁺-signaling network. Note that the peak magnitude is on the order of 10 µM, as compared to the individual peak magnitudes in Fig. 8 that are orders-of-magnitude lower. This shows that the nonlinearity inherent in the model causes the calcium signal resulting from coincident activation of the two inputs to be greater than the linear sum of the calcium signals caused by the individual inputs. This is consistent with experimental evidence (e.g., Wang et al. (3)). Additionally, the results show an IP₃-sensitive delay before the supralinear rise in [Ca²⁺], which is also consistent with experimental data (see Supplementary Material).

View larger version (17K): [in this window] [in a new window]   FIGURE 9  Relative timings of the individual Ca2+ signals in the spine cytosol (top) and the resulting supralinear Ca2+ spike in the spine cytosol during coincident activation (bottom). Data obtained from the compartmental model.

Sensitivity of the supralinear behavior
The results in this section will show that the ability of the system to demonstrate the strongly nonlinear behavior evidenced by the supralinear Ca²⁺ spike is strongly dependent on many features of the system, including biochemical effects and geometrical effects (i.e., spine geometry, as discussed previously). Recall that the IP₃R in Purkinje cells are more abundant and less sensitive to IP₃ than those found in other neurons or peripheral cells. Our earlier neuroblastoma cell model, for example, produced robust Ca²⁺ release with a and d_IP3 each an order-of-magnitude lower than the values used thus far in this work. By adjusting the parameters a (abundance of IP₃R) and d_IP3 (sensitivity of IP₃R to IP₃) found in Eq. 3, the effect of adjusting the IP₃R abundance and sensitivity can be investigated in terms of the ability of the system to generate the Ca²⁺ spike that induces LTD. In Fig. 10 a, the results of a simulation involving coincident activation of the PF and CF are shown for parameter values that represent an order-of-magnitude decrease in the IP₃R abundance and an order-of-magnitude increase in the IP₃R sensitivity. As can be seen in the figure, altering those parameters by an order of magnitude effectively extinguishes the supralinear [Ca²⁺] spike under these conditions. It is important to stress that a and d_IP3 were increased in our model to conform to the experimental findings of extraordinarily high IP₃R density (13) and low sensitivity to IP₃ (17). They were not adjusted simply to fit the desired Ca²⁺ signaling behavior.

View larger version (15K): [in this window] [in a new window]   FIGURE 10  (a) Spine cytosol Ca2+ signal during coincident activation when the abundance of IP3R is decreased by a factor of 10 and the sensitivity (dIP3) is increased by a factor of 10. The supralinear behavior is negligible under these conditions. Data obtained from the compartmental model. (b) Ca2+ signals in the spine and proximal dendrite cytosol during coincident activation. The top plot shows the responses given the normal Purkinje IP3R sensitivity, whereas the bottom plot shows the responses with the sensitivity to IP3 (dIP3) increased by a factor of 10. When the sensitivity is increased, a secondary Ca2+ spike is observed in the dendrite that diffuses into the spine. This effect may prevent localization of LTD to the intended spine. Data obtained from the compartmental model.

As alluded to in the previous subsection and shown by Wang et al. (3), the relative timing of the CF and PF inputs is critical to the formation of the supralinear Ca²⁺ spike. Experimental results show that the largest magnitude Ca²⁺ spikes are obtained by activating the CF briefly after the PF input train begins (3). As shown in Fig. 11, simulated results agree qualitatively with the experimental evidence in terms of both peak and integral fluorescence. To obtain results that could be correctly compared to the data of Wang et al. (3), the indicator dye was changed to Magnesium Green (K_D = 19 µM at a concentration of 375 µM) and the IP₃ pulse train length was reduced to include only four pulses. The results show that the maximum [Ca²⁺] is obtained when the CF is activated 0.05 s after the PF input train begins. Given timings of this order of magnitude, the CF [Ca²⁺] spike is occurring just before the time that [IP₃] is reaching its peak magnitude at the end of the PF pulse train, thus providing a large Ca²⁺ influx when the maximum release of Ca²⁺ from the ER is occurring. Given this relative timing, [Ca²⁺] peaks at approximately the same time that [IP₃] peaks. Timing shifts in either direction result in a decrease in maximum [Ca²⁺], with shifts in the CF activation time to earlier times resulting in a sharper decrease in maximum [Ca²⁺]. These results show a strong agreement with the experimental results of Wang et al. (3), in terms of the time when the maximum Ca²⁺ peak occurs. Also in qualitative agreement with the experimental data, the results in Fig. 11 show that the integrated Ca²⁺ data trend is more symmetric than the peak Ca²⁺ data trend. Analysis of the results suggests that proper timing of CF activation relative to PF activation may stimulate the CICR phenomenon at the IP₃R to a degree that is optimal in terms of yielding maximum Ca²⁺ release given the fixed IP₃ profile, as is necessary for induction of LTD.

View larger version (26K): [in this window] [in a new window]   FIGURE 11  (a) Maximum spine [Ca2+] observed during coincident activation as a function of the relative timing of the PF and CF inputs reported in terms of fluorescence. Timing is calculated as time of CF activation – time of initial PF activation. For comparison to data of Wang et al. (3), an IP3 pulse train consisting of four pulses was used with Magnesium Green as an indicator dye. The results show that the maximum fluorescence is reached when the CF is activated 0.05 s after the IP3 input train begins. These results are qualitatively similar to the fluorescence results reported by Wang et al. (3) in Fig. 6 b of their article (shown in b, copyright permission granted by Nature) in terms of the time at which the maximum calcium is observed. Data obtained from the compartmental model. (c) The total integrated Ca2+ signal (computed over the full transient period) also reaches its maximum value at a relative timing near 0.05 s, in agreement with the experimental data of Wang et al. as shown in d (copyright permission granted by Nature).

The influx (CF-activated) contribution to the total Ca²⁺ in the spine is a linear function of the density of channels, the open timing, and the current per channel. These are all held constant in our models. However, the IP₃R Ca²⁺ release is intrinsically nonlinear as it is sensitive to the level of Ca²⁺ in the system, with higher levels of Ca²⁺ inducing increased Ca²⁺ release. This positive feedback is part of the mechanism that leads to the Ca²⁺ spike. As evidence, simulations were performed under coincident activation conditions but with the model parameter that controls the sensitivity of the IP₃R to Ca²⁺ levels, d_Ca, increased to correspond to a lessened sensitivity to Ca²⁺. Given an order-of-magnitude increase in this parameter, the Ca²⁺ spike is extinguished. This result shows that the positive feedback in the CICR mechanism is necessary to generate the supralinear spike. In agreement with experimental evidence, the IP₃R with its concerted IP₃ and Ca²⁺ dependencies is necessary for the formation of the supralinear Ca²⁺ spike.

A second important aspect of the mechanism leading to the Ca²⁺ spike is the role played by Ca²⁺ buffers in the system. The Ca²⁺ buffers act as significant stores of bound Ca²⁺. Fig. 12 a is a plot of the fraction of buffer sites occupied for each of the buffers as a function of time during a normal coincident activation simulation. The model results shown in Fig. 12 a indicate that the buffers approach but do not reach saturation during the Ca²⁺ spike. Of course, as the buffers approach saturation, there must be a concomitant increase in free Ca²⁺ (Figs. 6, 9, and 12 b, bottom). The rate of buffer diffusion from the spine may also be important, since the bound form of the buffers are free to diffuse from the spine to the dendrite and to be replaced by buffers with free sites, thus providing additional Ca²⁺-buffering capacity. A useful estimate of the characteristic time for diffusion is given by L²/D_X, where L is a characteristic diffusion length and D_x is a diffusion coefficient. If L is taken to be the length of the spine neck and D_X is taken to be the fastest buffer diffusion coefficient in the model (D_PV = 43 µm²/s), the fastest characteristic diffusion time is 0.010 s. This timescale is the same order of magnitude as the timing of the Ca²⁺ spike formation observed in the data obtained in this work. This implies that, in the best case, the rate of buffer diffusion may not be fast enough to replenish the spine with unbound buffers, thus aiding in the saturation of the buffers present in the spine. The results in Fig. 12 a also suggest that the finite rate at which Ca²⁺ binds to the buffers may be slower than the rate of formation of the Ca²⁺ spike. This is particularly true in the case of PV, where the results suggest that Mg²⁺ must unbind first to free sites for Ca²⁺ to occupy since Mg²⁺ is preferentially bound. All of these effects suggest that the buffers are unable to bind all of the Ca²⁺ that is released into the cytosol thus resulting in the large, nonlinear increase in Ca²⁺ observed during coincident activation.

View larger version (20K): [in this window] [in a new window]   FIGURE 12  (a) Fraction of calbindin (top) and parvalbumin (bottom) buffer sites occupied during coincident activation in the spine cytosol. Data obtained from the compartmental model. (b) Total spine Ca2+ (bound and free) during coincident activation and independent activation of the PF and CF (top). Fraction of spine Ca2+ that is bound during coincident and independent activation (bottom). Data obtained from the compartmental model.

	TWO-DIMENSIONAL SPATIAL MODEL

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Although the results above demonstrate that the one-dimensional and compartmental models are sufficient to model the localized effects considered in this work, future research may consider effects that occur over longer length scales. This type of study will require a model that represents the geometry of the dendrite more accurately than the somewhat simplified one-dimensional model. In this section, a two-dimensional model of IP₃ diffusion in a segment of the dendrite is introduced and used to demonstrate its ability to yield results similar to those seen in the previous sections.

In the two-dimensional model, the PDE in Eq. 6 is solved across the geometry shown in Fig. 13. No approximations are made that would allow for the use of ODE expressions. Diffusion in spines is treated explicitly as in the dendritic shaft. The geometry is based on the high-magnification micrograph in Fig. 1. The Virtual Cell system allows electronic images such as this to be used in defining system geometries. As can be seen in the figure, the geometry used in this model can account for spines on either side of the dendritic shaft as opposed to the one-dimensional model that can only model a single spine at each point along the shaft. Because the one-dimensional model represented a 500-µm segment of dendrite where the boundaries could be safely assumed to have reached basal IP₃ concentrations at all times, the one-dimensional model was used to establish the time-dependent boundary conditions at the open ends of the two-dimensional geometry.

View larger version (64K): [in this window] [in a new window]   FIGURE 13  [IP3] after simulated activation of a PF obtained from the two-dimensional model. As seen in the one-dimensional model, high IP3 levels are not observed in neighboring spines, including those directly across the dendritic shaft from the activated spine.

Although the two-dimensional model was found to produce the same types of results as the compartmental model for analysis of events caused by stimulation of a single spine, this type of modeling will become necessary when consideration is given to events occurring in multiple spines. In that case, the results of this work will be instrumental in identifying parameter values and simulation conditions to be used in extending the low-dimensional models to higher-dimensional cases. The two-dimensional model discussed and demonstrated in this section is an example of such a model.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
The results obtained in this work demonstrate that a unique interplay of biochemical and morphological specializations contribute to the Purkinje cells' ability to generate localized Ca²⁺ spikes in dendritic spines. These features range from the biochemical effects of buffer capacity and the unique properties of the IP₃R to the diffusional barrier imposed by the spine neck. The simulation results suggest that the system is designed in such a way that a lack of any of the features and mechanisms discussed in this work may prevent the onset of LTD. In addition, the results obtained in this work provide some understanding of the degree of sensitivity of the supralinear Ca²⁺ response to many of the important system features.

In the case of the unusual IP₃R sensitivity and abundance observed in Purkinje cells, the results indicate that these features are necessary for the system to induce and localize the supralinear Ca²⁺ spike that is a prerequisite for LTD. The system geometry is also a critical feature in the Ca²⁺ signaling network due to the role it plays in limiting diffusion from the spine for the majority of the system species. In the case of IP₃, its higher diffusion coefficient allows it to overcome the spine neck's diffusion barrier during repeated PF activation. This is consistent with the lack of evidence of rapid degradation of IP₃ after PF activation ends.

Although the modeling and simulation performed in this work focused on events occurring at the level of a single dendritic spine, the ability of the Virtual Cell to model spatially varying systems using PDE models will allow for further study of events occurring over larger length scales, such as branchlet-wide Ca²⁺ signaling, which is generally dominated by voltage-dependent Ca²⁺ entry. The capability to model and simulate PDE systems efficiently will permit the study of events occurring at multiple spines including identification of particular events and mechanisms that cause the spread of LTD effects from one spine to others. Additionally, the role of the unique geometry of the dendritic arbor could be investigated in this context. The parameters identified in this work using the simplified model and existing experimental data will be very useful in extending the model so that Ca²⁺ dynamics on the scale of branchlet segments, or even full branchlets, will be possible. For example, an obvious next step would be to consider simulated activation of multiple spines in a region of dendritic shaft. Neglecting to consider diffusion of material from the spines (i.e., consideration of only compartmental models) would have served to prevent the applicability of the models to problems such as these.

Additionally, proper spatial modeling is required as one begins to consider accurately the role of the mechanisms leading to the production of IP₃. For example, PIP₂ is bound to the membrane and, therefore, its availability for reaction may be limited by lateral diffusion along the membrane. This would limit the Purkinje cell's ability to produce the large amount of IP₃ required to stimulate the IP₃R. Membrane diffusion capabilities will soon be available in the Virtual Cell, thus allowing this type of analysis to be performed in conjunction with a larger spatial modeling study including explicit consideration of PIP₂ hydrolysis.

	APPENDIX

TOP ABSTRACT INTRODUCTION MODEL DEVELOPMENT ONE-DIMENSIONAL SPATIAL MODEL COMPARTMENTAL MODEL TWO-DIMENSIONAL SPATIAL MODEL CONCLUSIONS APPENDIX SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
This Appendix includes a summary of each of the models used in this work. All parameter values can be found in Table 1.

One-dimensional model
In this model, each species X (X = Ca²⁺, IP₃, PV, PVB_c, PVB_m, CD28k, CD28kB₁, CD28kB₂, and CD28kB₁₂) at time t and location z is represented by two concentrations: [X]_d(z,t) and [X]_s(z,t) in the dendrite and spine, respectively. These concentrations are governed by the equations

where D_X is the diffusion coefficient of the species X, l₁₂ is the spine neck length, and the coefficients k₁ and k₂ are determined by the spine neck radius r_n, the spine radius r₁, the dendrite radius r₂, and the linear spine density s: k₁ = s(r_n/r₂)², k₂ = (3/r₁)(r_n/2r₁)². Th