Dendritic spine geometry and spine apparatus organization govern the spatiotemporal dynamics of calcium

Dendritic spines, small protein- and actin-rich protrusions located on dendrites of neurons, have emerged as a critical hub for learning, memory, and synaptic plasticity in both short-term and long-term synaptic events (Bourne and Harris, 2008; Rangamani et al., 2016). These subcompartments provide valuable surface area for cell–cell interaction at synapses, and compartmentalization of signaling proteins to control and process incoming signals from the presynaptic terminal (Nishiyama and Yasuda, 2015; Yasuda, 2017). Thus, dendritic spines are hotbeds of electrical and chemical activity. Since calcium is the first incoming signal into the postsynaptic terminal, calcium temporal dynamics have been extensively studied experimentally and more recently computationally (Denk et al., 1996; Augustine et al., 2003; Bloodgood and Sabatini, 2005; Bartol et al., 2015b; Yasuda, 2017). In particular, calcium acts as a vital second messenger, triggering various signaling cascades that can lead to long-term potentiation, long-term depression, actin cytoskeleton rearrangements, and volume expansion, among other events (Holmes, 1990; Bourne and Harris, 2008; Rangamani et al., 2016).

Dendritic spine activity has numerous timescales, with signaling pathways operating on the millisecond to the hour timescale following spine activation (Yuste et al., 2000; Segal, 2005; Rangamani et al., 2016). Calcium dynamics are on the millisecond timescale, since calcium is the second messenger that floods the spine following the release of neurotransmitter from the presynaptic terminal. The temporal dynamics of calcium have provided valuable insight into the signaling dynamics in dendritic spines, and it is quite clear that calcium dynamics are influenced by a large number of factors. Multiple studies have connected the electrical activity of the plasma membrane (PM) voltage to calcium dynamics of N-methyl-D-aspartate receptors (NMDARs; Jahr and Stevens, 1993; Shouval et al., 2002; Rackham et al., 2010). The electrophysiology of dendritic spines influences many signaling dynamics through voltage-sensitive (or voltage-dependent) ion channels (Jaffe et al., 1994), and thus, models of these dynamics can be linked to downstream signaling.

Calcium is critical for almost all the reactions in the brain (Augustine et al., 2003) and is believed to accomplish a vast variety of functions through localization (Clapham, 1995, 2007; Yuste and Denk, 1995). One possible way to achieve localization is by restricting the distance between the calcium source (often a channel) and the sink (calcium sensors and buffers). Thus, the localization of calcium can result from the location and mobility of different buffers and sensors (Schmidt, 2012). Spatial models of calcium dynamics in dendritic spines that consider such effects have been proposed previously (Holcman et al., 2004; Means et al., 2006; Bartol et al., 2015b). Spatial-temporal models of calcium dynamics have highlighted the role of calcium-induced cytosolic flow and calcium influx regulation by Ca²⁺-activated K⁺ channels (SK channels; Holcman et al., 2004; Rackham et al., 2010). Computational studies using stochastic models of calcium transients have revealed that the readout of fluorescent probes can alter experimental readouts (Yuste et al., 2000; Higley and Sabatini, 2012; Bartol et al., 2015b). In particular, fluorescent probes can sequester calcium, effectively acting as calcium buffers and lowering the perceived calcium concentrations. Therefore, computational studies have already provided some insight into the spatiotemporal dynamics of calcium in both dendritic spines (Yuste et al., 2000; Segal, 2005; Rangamani et al., 2016) and whole neurons (Loewenstein and Sompolinsky, 2003). In this study, we specifically focus on the effect of spine geometry on calcium signals in the postsynaptic spine.

Recent advances in imaging and reconstruction techniques have shed light into the complex surface area of a single spine and the ultrastructure of the spine apparatus (SpApp; Kuwajima et al., 2013; Bartol et al., 2015a,b; Wu et al., 2017). Experimental evidence shows that calcium signals remain predominantly localized to single spines (Holmes, 1990; Koch and Zador, 1993; Sabatini et al., 2002). Dendritic spines have a unique set of shapes and recently the size and shape distribution of spines has been linked to their physiological function (Hering and Sheng, 2001; Berry and Nedivi, 2017). Additionally, only ∼14% of spines have a specialized ER known as the SpApp (Basnayake et al., 2018 ,Preprint; Jedlicka et al., 2008), which serves as a calcium store (Wuytack et al., 2002; Jedlicka et al., 2008). Indeed, calcium dynamics in dendritic spines are quite complex.

Given the importance of calcium dynamics in dendritic spines and the complexity of spine ultrastructure (Spacek and Harris, 1997; Holbro et al., 2009; Wu et al., 2017), we sought to use a computational approach to probe the role of spine geometry in modulating the spatiotemporal dynamics of calcium. Specifically, we seek to address the following questions: (i) How does the size and shape of both the spine head and SpApp affect calcium dynamics? (ii) How does the distribution of channels along the synaptic membrane modulate calcium dynamics within the spine? (iii) How do calcium buffers and calcium diffusion rates affect the spatiotemporal dynamics of calcium? To answer these questions, we develop a spatial 3-D reaction-diffusion model with multiple compartments. We chose a computational approach because it is not yet possible to manipulate spine size, shape, or SpApp location with precise control experimentally. However, the insights obtained from computational approaches can lay the groundwork for generating experimentally testable predictions (Kotaleski and Blackwell, 2010).

To interrogate the spatiotemporal dynamics of calcium in dendritic spines, we developed a reaction-diffusion model that accounts for the fluxes through the different sources and sinks shown in Fig. 1. We briefly discuss the assumptions made in developing this model and outline the key equations below. See Table 1 for common notation.

We model a single dendritic spine of a rat hippocampal area CA1 pyramidal neuron as an isolated system, because we focus on the 10- to 100-ms timescale and the ability of the spine neck to act as a diffusion barrier for calcium (Yuste and Denk, 1995; Bloodgood and Sabatini, 2005; Byrne et al., 2011; Gallimore et al., 2016). As a result, we do not consider calcium dynamics due to the mitochondria. This assumption is valid in our model, because even though mitochondria are known to act as calcium stores, their dynamics are on a longer timescale (10–100 s; Wacquier et al., 2016) and mitochondria are located in the dendrite outside of the dendritic spine. Although, it is well known that the SpApp acts as a source for calcium, Ca²⁺ release from the SpApp occurs at longer timescales because of IP₃ receptors, RYR, and CICR (Jedlicka et al., 2008). It should also be noted that not all neurons have RYR and IP₃R on the SpApp (Mattson et al., 2000). Therefore, for the purpose and timescale of this study, we do not focus on CICR and therefore do not include RYR or IP₃R dynamics in this study.

We model membrane voltage as an algebraic equation based on the summation of an excitatory postsynaptic potential (EPSP) and a back-propagating action potential (BPAP) applied to the entire PM (Jahr and Stevens, 1993; Shouval et al., 2002; Griffith et al., 2016). The EPSP arrives 2 ms before the BPAP to stimulate the maximum possible membrane depolarization (Fig. 1 d; Hu et al., 2018). This stimulus triggers the influx of calcium into the spine head.

The average spine head volume is ∼0.03 µm³ (Spacek and Harris, 1997; Bartol et al., 2015b), but a large variation has been observed physiologically. The commonly observed shapes of dendritic spines include filopodia-like, stubby, short, and mushroom-shaped spines (Bourne and Harris, 2008; Berry and Nedivi, 2017). In this work, we consider two idealized geometries: a spherical spine head to represent a younger spine and an ellipsoidal spine head to represent a more mature mushroom spine (Spacek and Harris, 1997). The postsynaptic density (PSD) is modeled as a section of the membrane at the top of the spine head (Fig. 1 c). In simulations where the spine size is varied, we assume that the PSD changes surface area approximately proportionally to the spine head volume (Arellano et al., 2007).

SpApps are found primarily in larger mushroom spines (Spacek and Harris, 1997), which hints at their role in potential regulation of sustained spine volume (Ostroff et al., 2017). We assume that the SpApp acts as a calcium sink within the timescale of interest (Segal et al., 2010). Another assumption is that the SpApp has the same shape as the spine head, a simplification of the more complicated and intricate SpApp geometry (Spacek and Harris, 1997).

To maintain our focus on the short-timescale events associated with calcium dynamics, we include the following PM sources and sinks of calcium: voltage-sensitive calcium channels (VSCC), NMDARs, PM Ca²⁺-ATPase (PMCA) pumps, and sodium–calcium exchanger (NCX) pumps. NMDARs are localized on the postsynaptic membrane adjacent to the PSD, designated at the top of the spine head. VSCC, PMCA, and NCX pumps are uniformly distributed along the PM, including at the base of the spine neck. Therefore, we model the dendritic spine as an isolated system with the spine neck base modeled in the same manner as the rest of the PM rather than explicitly modeling the base with an outward flux into the dendrite (see the following assumption).

We model the spine as an isolated system (Yuste and Denk, 1995; Sabatini et al., 2002; Bloodgood and Sabatini, 2005; Lee et al., 2009). Therefore, for most of our analyses, we use the same boundary conditions at the base of the spine neck as the rest of the PM. In the online supplemental material, we relax this assumption and test different boundary conditions, including a clamped calcium concentration at the base of the neck and an explicit effect of a dendritic shaft attached to the spine neck (Figs. S4 and S5).

We explicitly model calcium in the cytoplasm and in the SpApp. We assume that the calcium concentration in the extracellular space (ECS) is large enough (2 mM; Clapham, 1995; Bartol et al., 2015b) that the calcium influx into the spine has an insignificant effect on ECS calcium concentration. Therefore, ECS calcium concentration is assumed constant. We only solve the volumetric reactions in the cytoplasm. The ER calcium concentration is assumed to only affect the fluxes on the ER membrane.

There are numerous CBPs present in the cytoplasm that act rapidly on any free calcium in the spine head (Yuste and Denk, 1995; Yuste et al., 2000; Sabatini et al., 2001, 2002; Bartol et al., 2015b). These CBPs are modeled as both mobile and fixed buffers in our system. Mobile buffers are modeled as volume components in the cytoplasm (Schmidt and Eilers, 2009; Schmidt, 2012), and they are modeled as a diffusive species with mass-action kinetics in the cytoplasm. We assume that the mobile buffers have a buffering capacity, κ, of 20, where κ = [B_m]/K_d (Sabatini et al., 2002; Matthews and Dietrich, 2015; B_m is the mobile buffer concentration; K_d is the dissociation constant). Dendritic spines also have fixed or immobile buffers, but the molecular identity of fixed buffers remains more elusive. Studies suggest that they are primarily membrane-bound components (Matthews and Dietrich, 2015); therefore, we model fixed buffers as immobile species localized to the PM. As a result, the interactions of calcium with these fixed buffers are treated as flux boundary conditions for the surface reactions with mass-action kinetics. We assume fixed buffers have a K_d = 2 µM (Bartol et al., 2015b) and have a concentration of 78.7 µM (Bartol et al., 2015b). These values are converted to a membrane density by multiplying by the spine volume over PM surface area (n_PMr; see online supplemental material section S1.1 for more details.).

In this model, the SpApp acts as a Ca²⁺ sink in the 10- to 100-ms timescale (Fifková et al., 1983; Fifková, 1985; Jedlicka et al., 2008). The implications of this assumption are discussed later in the paper. We assume that SERCA pumps are located uniformly on the SpApp membrane. SERCA pumps have been observed to buffer a large percentage of calcium within the spine, so we include SERCA pumps with a relatively large influx (Higley and Sabatini, 2012; Hu et al., 2018). We also include a small leak current from the SpApp to the cytoplasm and set this leak current to offset pump dynamics at basal calcium concentrations of 100 nM (Bartol et al., 2015b; Futagi and Kitano, 2015).

Based on these assumptions, we constructed a 3-D spatial model of Ca²⁺ dynamics in dendritic spines. Our control geometry is a medium-sized spine with volume of ∼0.06 µm³ including the spine head and neck, with a SpApp of volume ∼0.003 µm³. We use a spherical spine with spherical SpApp and ellipsoidal spine with ellipsoidal SpApp as our two control spines of interest. Most results are shown as a 2-D cross section for ease of interpretation (see Fig. S7 for examples of the full 3-D solutions).

The spatiotemporal dynamics of calcium are determined by the combination of dynamics within the spine volume and boundary conditions at the PM and SpApp. Calcium dynamics in the spine volume are represented by a single reaction–diffusion equation:

Here, D is the diffusion coefficient of calcium and $f(Cacyto,2+CDP)$ is a function that represents the reaction between cytosolic calcium and mobile buffers in the cytoplasm; ∇² is the Laplacian operator in three dimensions. The stimulus to the system is the depolarization of the membrane based on an ESPS and BPAP separated by 2 ms (Fig. 1 d). The boundary conditions at the PM and the SpApp are given by time-dependent fluxes that represent the kinetics of different channels, pumps, and fixed buffers.

We model the calcium influx through activated NMDARs in response to glutamate release in the synaptic cleft and the calcium influx through VSCCs in response to membrane depolarization (Jahr and Stevens, 1993; Shouval et al., 2002; Bartol et al., 2015b). We should note that the majority of existing models for NMDAR and VSCC calcium influx assume well-mixed conditions. In this model, these species are restricted to the PM, which is the boundary of the geometry. This results in a time-dependent flux at the PM. Both the NMDAR and VSCC-mediated calcium influx depend on the membrane voltage (see Fig. 1 d); we prescribe this voltage as a set of biexponentials to capture a BPAP and EPSP based on previous studies (Jahr and Stevens, 1993; Shouval et al., 2002). On the PM, we also include PMCA and NCX that are activated in response to a change in cytosolic calcium concentration (Calizo et al., 2017 ,Preprint; Maurya and Subramaniam, 2007). We also localize fixed CBPs (fixed buffers) to the PM. Therefore, the binding of cytosolic calcium to fixed buffers (B_f) is modeled as a membrane flux, J_Bf. The flux boundary condition at the PM is then the sum of all these fluxes and is given by

The functions that define the flux terms in Eq. 2 are given in Table S1.

In the cases where we included a SpApp, we included SERCA pumps and a leak term along the SpApp membrane that are functions of the cytosolic calcium concentration. The boundary condition for the flux across the SpApp membrane is given by

The functions that define the flux term in Eq. 3 are given in Table S1.

To briefly summarize, these governing equations are simply the balance equations that keep track of the spatiotemporal dynamics of cytosolic calcium due to calcium diffusion and mobile buffers (Eq. 1), influx and efflux through the PM (Eq. 2), and influx and efflux through the SpApp membrane (Eq. 3). The coupled nature of this system of equations and time-dependent fluxes limits the possibility of obtaining analytical solutions even for simple geometries (Cugno et al., 2018  Preprint). Therefore, we use computational methods to solve these equations.

Given the vast number of parameters in this model, we constrain the parameters in our model as follows: parameter values are chosen from experimental observations or existing computational models or to match overall experimental and computational observations with respect to pump or channel dynamics (section S1). Overall, we predict a high calcium concentration and relatively fast decay dynamics (Sabatini et al., 2002; Higley and Sabatini, 2012; Hu et al., 2018). See Fig. S2 for a comparison of temporal dynamics to existing literature. We conducted a kinetic parameter sensitivity analysis for our model using the open source software COPASI—a COmplex PAthway SImulator (Supplemental text). The sensitivity analysis was performed in COPASI by converting our spatial model into a compartmental ordinary differential equation (ODE) system. This conversion involves transforming boundary flux equations into volumetric reaction rate through the lengthscale factor n, the volume-to-surface area ratio.

We modeled the dendritic spines using idealized geometries of spheres and ellipsoids; dendritic spines consist of a spine head attached to a neck, with a similarly structured SpApp within the spine, see Fig. 1 c for the different model geometries used in this study. These geometries were inspired by reconstructions (Bartol et al., 2015a; Griffith et al., 2016; Paulin et al., 2016) and were idealized for ease of computation. For the variations of the base of the neck, we also include a condition with an explicit dendrite modeled as a cylinder attached to the spine neck, Figs. S4 and S5. The geometric parameters, including volume and surface area, are given in Table S8.

Simulations for calcium dynamics were conducted using commercial finite-element software (COMSOL Multiphysics 5.4). Specifically, the general form and boundary partial differential equations (PDEs) interface were used and time-dependent flux boundary conditions were prescribed. A user-defined tetrahedral mesh with a maximum and minimum element size of 0.0574 µm and 0.00717 µm, respectively, was used. Due to the time-dependent, nonlinear boundary conditions used in this model, we also prescribed four boundary layers (prism mesh elements) on all membranes in COMSOL. A time-dependent solver was used to solve the system, specifically a MUMPS (MUltifrontal Massively Parallel sparse direct Solver) solver with backward differentiation formula time-stepping method with a free time stepper. Results were exported to MATLAB for further analysis. All COMSOL files will be posted on the Rangamani Lab website for public dissemination.

All model equations and parameters and geometric parameters can be found in the online supplemental material. Fig. S1 shows the effects due to using a different lengthscale, n. Fig. S2 shows temporal comparison with other existing models. Fig. S3 is a sensitivity analysis of model parameters. Figs. S4 and S5 show the case when the neck base has different boundary conditions in spherical and ellipsoidal spines, respectively. Fig. S6 shows the effect of spine neck radius. Fig. S7 shows calcium spatial dynamics in 3D for spherical and ellipsoidal spines. Fig. S8 shows temporal dynamics and peak concentrations for spines of various shapes. Fig. S9 shows the effects of different calcium buffers in ellipsoidal spines. Table S1 lists the main equations in the model. Table S2 lists parameters used in the model for the volumes. Table S3 lists parameters for NMDAR. Table S4 lists parameters for membrane voltage. Table S5 lists parameters for VSCC. Table S6 lists parameters for PMCA and NCX pumps. Table S7 lists parameters for SERCA. Table S8 lists the size and shape parameters of the different geometries. Table S9 lists size variations for the spherical spine head with spherical SpApp. Table S10 lists size variations for the ellipsoidal spine head with ellipsoidal SpApp. Table S11 lists size variations for the spherical spine head with spherical SpApp when varying SpApp volume. Table S12 lists size variations for the ellipsoidal spine head with ellipsoidal SpApp when varying SpApp volume.

Using the model developed above (Eqs. 1–3), we investigated how different geometric factors of the spine head and SpApp affect calcium dynamics. Because of the coupling between the volume dynamics (Eq. 1) and the fluxes on the membranes due to biochemical components (Eqs. 2 and 3), the effect of spine geometry on calcium dynamics is quite complex. To parse the coupled effects, we have categorized and organized our simulations as follows and discuss each case in detail. First, we investigated the effect of spine volume to surface area ratio. This parameter can be changed in multiple ways: by changing the shape of the spine head and the shape or presence of the SpApp (Figs. 2 and 3), by changing the size of the spine head alone (Fig. 4), or by changing the size of the SpApp alone (Fig. 5). Next, we investigated the effect of spatial distribution of the fluxes on the PM (Fig. 6). Finally, we demonstrate the effect of buffer types and their location (Figs. 7 and 8). We describe each of these results in detail below.

We first analyzed how the volume to surface area ratio of the spine affects the spatiotemporal dynamics of calcium by simulating calcium dynamics in the different geometries shown in Fig. 1 c. All of these geometries were constructed such that they have the same volume of the spine cytosol but the surface area of the PM and SpApp vary because of the shape (Table S8). We note that in all the geometries, the temporal dynamics of calcium shows a rapid increase in the first 2–3 ms and a decay over ∼50 ms (Fig. S8). This time course is consistent with experimentally observed and computationally modeled calcium dynamics (Ngo-Anh et al., 2005; Matthews et al., 2013; Hu et al., 2018; Figs. 2 d and S2, respectively). The spatial profiles of calcium in spheres and ellipsoids (Fig. 2, a and b) show that spine shape can alter the spatial gradients and decay profiles of calcium. In particular, we observe that while at 5 ms all shapes demonstrate a gradient from the PSD region to the spine neck, because of the localized influx of calcium through the NMDAR in the PSD, at 10 ms, the ellipsoidal spine heads have almost no gradient within the spine head. Instead, they have a more pronounced difference between the spine head and spine neck. We also changed the shape of the SpApp to get different combinations of spine head and SpApp geometries (Fig. 2). Since the apparatus acts as a sink, we find that the reduction in surface area to volume ratio of the spherical apparatus leads to less influx into the SpApp, compared with the ellipsoidal SpApp (Fig. 2 c). This result emphasizes the nonintuitive relationships between organelle shape and spine shape.

In addition to the transient response of calcium, the cumulative calcium (total calcium) also carries information with respect to synaptic plasticity (Bourne and Harris, 2008; Rangamani et al., 2016; Basak and Narayanan, 2018). Therefore, we calculated the integrated calcium over the entire spine volume (area under the curve [AUC]; Gorman and Sejnowski, 1988; Heinrich et al., 2002; Atay and Skotheim, 2017) over 300 ms. Calcium AUC at 300 ms (Fig. 2 c) shows that all spines have slightly different accumulated calcium. Upon closer examination of the membrane fluxes of the sphere and ellipsoid spines (Fig. 2, e–m), we see that there is a complex relationship between the calcium ion concentration and nonlinear flux terms. In particular, the higher calcium influx in the ellipsoids due to the larger surface area leads to a higher efflux of calcium through the PMCA, NCX, and SERCA pumps and increased binding to fixed buffers on the PM. Therefore, the nonlinear effects of the fluxes associated with the pumps reduces the difference in calcium between spines of different shapes.

Another way to modulate the volume to surface area ratio is to consider spines with and without the SpApp. In our model, the SpApp serves both as a calcium sink through SERCA pumps and as an excluded volume. We observe that the presence of the SpApp in both spheres and ellipsoids (Fig. 3, a and d, top row) results in a steeper gradient from the PSD to the neck when compared with the spines without a SpApp (Fig. 3, a and d, bottom row). Additionally, at 10 ms, spines without a SpApp have a higher calcium concentration than spines with a SpApp. This is because the dynamics of calcium are altered in the presence of the SpApp by the SERCA fluxes, Eq. 3. As a result, regardless of the shape, spines with a SpApp have a faster decay of cytosolic calcium (Fig. 3, b and e). We found that in both geometries, the AUC of spines with a SpApp was lower at different time points when compared with spines without a SpApp (Fig. 3, c and f). At 300 ms, the spine without a SpApp has 47.75% more total calcium for the spherical geometry and 49.09% more for the ellipsoidal geometry. Thus the flux due to SERCA plays a significant role in altering the decay dynamics of calcium in spines with SpApp and as a result alters the total calcium.

In addition to spine shape, spine size (volume) is also known to change during development and plasticity related events (Knott et al., 2006; Kasai et al., 2010; Rangamani et al., 2016). How does changing the volume of the spine, while maintaining the same shape, affect calcium dynamics? To answer this question, we conducted simulations in spherical and ellipsoidal spines of different volumes. Recall that the control spine has a cytosolic volume, V_cyto, of 0.06 µm³ (Table S8). For each geometry (sphere and ellipsoid), we maintained the same size and shape of the SpApp as before and only changed the spine cytoplasm volume in the range of 0.5 V_cyto to 1.5 V_cyto (Tables S9 and S10). We found that the relationship between spine volume and calcium concentration is inversely proportional. For spine volumes smaller than the control, we observed an increase in calcium concentration for both geometries, whereas for larger volumes, calcium concentration decreases (Fig. 4, a and b). As expected, we found that for both geometries, an increase in spine volume resulted in an increase in cumulative calcium (Fig. 4 c). Furthermore, we found that the change in cumulative calcium has a direct but nonlinear relationship with the change in spine volume. For the range of volumes investigated, the peak calcium concentration and AUC show an exponential relationship with respect to volume. We see at all sizes the ellipsoid has higher peak concentrations but lower AUC compared with the sphere.

The complex architecture of the SpApp was recently elucidated in a focused ion beam scanning electron microscopy study by Wu et al. (2017). Since we cannot yet manipulate the shape of the SpApp in vivo, we varied the geometric features of the SpApp in silico to see how they affect calcium dynamics. Previously, we showed that SpApp shape and spine volume separately can alter the AUC and peak calcium. Here, for a given spine shape, we varied the volume of the SpApp to modulate the cytosolic volume to SpApp volume ratio. In this case, by varying the SpApp volume, we altered the spine volume to be 50% and 75% of V_cyto, the control spine volume. Here, we change the spine volume by the changing SpApp size (Fig. 5). We found that a larger SpApp leads to a decrease in calcium concentrations and AUC (Fig. 5, a and c) but an increase in peak concentration (Fig. 5 b). We note that as SpApp volume increases, AUC drops in a nonlinear manner in both geometries. The spherical spine has a 55.4% and 77.8% reduction in AUC from control for the 75% and 50% spines, while the ellipsoidal spine has a 57.2% and 79.4% reduction from control. For all sizes, the ellipsoid shows higher peak concentrations but lower AUC compared with the sphere. This effect is in part due to the shortened distance between the PSD and SpApp in the ellipsoid. The distance between the PSD and SpApp is an important lengthscale in the spine, since it controls the distance between the sources and sinks of calcium. Therefore, changing this lengthscale, as happens when changing the SpApp volume, influences calcium dynamics in nonlinear ways. From these observations, we conclude that spine volume coupled with SpApp volume and surface area is an important regulator of calcium dynamics in dendritic spines.

An intriguing feature of the SpApp in particular (Fig. 5 d; Wu et al., 2017) and the ER in general (Nixon-Abell et al., 2016) is the large surface-to-volume ratio occupied by this organelle. Therefore, we next considered the effect of the volume-to-surface area ratio (n; given in units of length) of the SpApp (Fig. 5, d–f). We modeled the boundary flux on the SpApp membrane such that this flux is proportional to n_SpApp (volume-to-surface area ratio for the SpApp). As a result, when we increase the “volume” of the SpApp by increasing n_SpApp, calcium flux into the SpApp will increase. We noticed that at lower n_SpApp values, the peak calcium concentration and to a less obvious extent the AUC (Fig. 5, e and f) plateau but decrease substantially at larger n_SpApp values.

From these observations, we conclude that the SpApp acts as a physical and spatial buffer for calcium dynamics by regulating the timescale through surface to volume regulation in the interior of the spine. The SpApp acts as a calcium sink in the timescale of interest (Fifková et al., 1983; Fifková, 1985; Jedlicka et al., 2008), and in the absence of the SpApp, the only way to remove calcium from this system is through CBPs and pumps. Furthermore, since the SpApp has been known to grow and retract from the dendritic spine in response to stimuli (Deller et al., 2006), regulation of SpApp surface area can also allow for rescaling calcium dynamics in the spine (Wilson et al., 1983).

The density of VSCCs and number of NMDARs that open in response to stimuli in dendritic spines varies (Sabatini and Svoboda, 2000). One of the primary determinants of calcium dynamics in the spine is the various membrane fluxes, because these fluxes serve as sources and sinks at the PM and sinks at the SpApp membrane. We investigated the effect of the spatial distribution of membrane fluxes on calcium dynamics in the dendritic spine by considering either only NMDARs or VSCCs as the calcium source (Fig. 6). We observed that if only NMDAR activity was present, then calcium concentration was high in the PSD region due to the localization of NMDARs to the PSD but the overall calcium concentration was small regardless of the spine head shape (Fig. 6 a). However, if only VSCCs were active, then the spatial gradient of calcium is mainly between the spine head and spine neck (Fig. 6 a). The temporal dynamics are also affected by the receptor and channel distributions (Fig. 6 b). With only VSCC present, there is a larger calcium peak, but with faster decay. With only NMDARs, we observe a lower calcium peak concentration but a prolonged transient. Therefore, membrane flux distribution can impact the spatial and temporal dynamics of calcium in a nonlinear manner. This agrees with experimental results stating that the various calcium sources behave supralinearly (Yuste and Denk, 1995) and a balance between various calcium fluxes is required for tightly regulating calcium concentrations in these small volumes.

Due to the vast number of buffers that are known to affect calcium dynamics, we investigated the effect of four different buffer conditions: (i) both fixed membrane-bound buffers and mobile cytosolic buffers (control), (ii) fixed buffers localized to the membrane, (iii) mobile buffers in the cytoplasm, and (iv) a uniform exponential decay applied across the whole cytoplasm (Figs. 7 and S9). We observe similar spatial dynamics for all buffer types (Fig. 7 a), but temporal dynamics differ greatly (Fig. 7 b) and, as a result, alter accumulated calcium (Fig. 7 c). We see that the decay behavior of calcium, and therefore total calcium, is highly dependent on buffer type. Peak calcium follows the same trend as AUC, with mobile buffers having the highest calcium concentrations and total calcium, then exponential decay, fixed buffers, and finally the control. We also consider how reaction dynamics versus diffusion rate govern calcium dynamics because buffer dynamics and diffusion rates of calcium are coupled (Yuste et al., 2000). We varied the diffusion rate of calcium and the concentration of mobile buffers to quantify how calcium dynamics are reaction or diffusion controlled (Fig. 8). We see that diffusion controls the spatial gradient seen within the spine (Fig. 8 a), while mobile buffer concentration controls the lifetime of the calcium transient (Fig. 8 b). Combining these effects, we see that the buffer concentration variation has the greatest effect at lower diffusion rates (bottom row of Fig. 8, a and b, and Fig. 8 c). The peak concentration of calcium is almost entirely dependent on the diffusion rate (Fig. 8 d). Therefore, based on the high diffusion rate of Ca²⁺ reported in the literature, we expect the system to be in a diffusion-dominated regimen.

Calcium is a fundamental player in both neuronal (cellular) and neural (systems) functionality (Siesjö, 1990; Keener and Sneyd, 1998; Yuste, 2010). Compartmentalized by dendritic spines, calcium has a vital role in triggering signaling pathways for long-term potentiation, long-term depression, synaptic plasticity, and other processes associated with learning and memory formation (Rangamani et al., 2016). However, while dendritic spines are known to form functional subcompartments, it is less understood how this specialized compartmentalization impacts calcium dynamics (Yuste et al., 2000). In this study, we explored the intricate relationship between calcium dynamics and the shape and size of dendritic spine structures. We found that while the relationship between spine geometry and calcium dynamics is quite complicated (Hering and Sheng, 2001; Rochefort and Konnerth, 2012; Berry and Nedivi, 2017), some general conclusions can be drawn from our study.

First, the volume-to-surface ratio, rather than the shape and size itself, seems to have a dramatic effect on spine calcium (Figs. 2, 3, 4, and 5). Of course, the volume-to-surface ratio itself can be dramatically altered by size, shape, and internal organization as a many-to-one function. Then, we can think of the ultrastructural organization of the dendritic spine (Spacek and Harris, 1997) as perhaps “optimized” to not only increase contacts with neighboring axons and neural circuit connection (Yuste and Denk, 1995; Yuste et al., 2000; Yuste, 2010) but also tune this volume-to-surface ratio dynamically (Murakoshi and Yasuda, 2012). This volume-to-surface ratio coupling further highlights the complex relationship between spatial sources and sinks of calcium, which becomes apparent when the distance between the spine PM and internal organelle becomes quite small. We note that in our model, we assume constant pump density, which highlights the volume-to-surface area ratios between various shapes (Matsuzaki et al., 2001; Noguchi et al., 2005). Experimental results have already shown different behavior in large versus small dendritic spines (Paulin et al., 2016), and additional studies on dendritic spine geometry have shown that stable, mature spines are usually larger spines that tend toward mushroom shapes as they grow around adjoining axons and are more likely to have a SpApp (Spacek and Harris, 1997). In comparison, younger, less stable spines tend to be smaller and more spherical (Spacek and Harris, 1997; Berry and Nedivi, 2017). Therefore, we predict that spine size and SpApp presence are coupled to control calcium dynamics. This result should be investigated further, in particular to make predictions on why stable spines tend to be larger and mushroom shaped. The inverse relationship between the volume-to-surface ratio that we found and a possible exponential relationship suggest a potential limiting mechanism for maintaining homeostasis of synaptic potentiation (Lee et al., 2010, 2012; Béïque et al., 2011; Turrigiano, 2011). By altering the spine size dynamically (Murakoshi and Yasuda, 2012) and the presence and absence of the SpApp dynamically (Deller et al., 2006), spines could maintain their optimal range of synaptic function (Lee et al., 2010, 2012; Béïque et al., 2011; Turrigiano, 2011).

Second, localization of membrane fluxes alters calcium transients (Figs. 5 and 6). These fluxes, which serve as boundary conditions, can be altered by changing the density and distribution of calcium sources and sinks. This idea is consistent with how calcium signal localization is a result of tuning the distance between sources and sinks (Augustine et al., 2003). We show here that in addition to distance, the strength of the fluxes is important. Thus, in various disease states that impact the distribution or strength of membrane components, we predict atypical spatial calcium gradients are possible that could impact downstream signaling pathways. For example, NMDAR dysfunction, whether leading to increased or decreased functionality, can potentially lead to central nervous system diseases such as stroke, Huntington’s disease, or schizophrenia (Zhou and Sheng, 2013).

Finally, the role buffers play in modulating calcium transients is not only by changing the decay time as previously thought but also by tuning the membrane fluxes, especially in the case of fixed buffers (Figs. 7 and 8; Higley and Sabatini, 2012; Matthews and Dietrich, 2015). Again, the timescale that we see is a combination of rate alterations at the membrane and rate alterations in the volume resulting in broader control of calcium dynamics. The crowded environment within the spine head also has consequences for calcium diffusion, and while it is possible for calcium to diffuse through a crowded space, particularly in the PSD, the exact mechanisms of such transport remain unclear (Santamaria et al., 2006; Hotulainen and Hoogenraad, 2010; Byrne et al., 2011). Thus, our study highlights the need for connecting biophysical features of the spine and molecule localization to the dynamics of calcium (Fig. 1).

We also note that our model has certain limitations. In particular, within the crowded environment of the dendritic spine cytoplasm is an abundance of actin, which has previously been shown to have the potential to create cytosolic flow through contraction following spine activation, leading to faster calcium diffusion (Holcman et al., 2004). We do not address this spine contraction in this model, but actin contributions are a focus of ongoing research in our group. While we touched upon the role of diffusion and CBPs, we also acknowledge that much work remains to be done on the true impact of the dense actin network and crowded environment within dendritic spines in regulating these processes (Ouyang et al., 2005). In addition, we modeled isolated spines, but the width of the spine neck has been showed as an important determinant of calcium dynamics when comparing larger and smaller spines (Noguchi et al., 2005; Arellano et al., 2007; Araya et al., 2014) and could play into communication to the dendrite. It is also possible that stochastic modeling will give better quantitative insights without altering the underlying physics (Kotaleski and Blackwell, 2010; Bartol et al., 2015b; Voorsluijs et al., 2019). The development of a combined stochastic–deterministic model can help combine these two regimes to address the fundamental physics that occurs in these small systems with complex membrane dynamics and few molecule situations.

Despite these shortcomings, we have identified several key features of the relationship between dendritic spine geometry and calcium dynamics. Current models of synaptic weight updates use calcium as the determinant for the synaptic weight vector and the learning rate (Malenka et al., 1988; Cummings et al., 1996; Shouval et al., 2002; Yeung et al., 2004). Here, we show that calcium in a spine, even in short timescales, is a function of geometry, ultrastructure, and buffers. Based on these insights, we speculate on what the biophysical features of the spine mean for neural systems-level functionality. It has long been considered that a neural circuit level model of synaptic weight updates can be informed by the calcium transients in the synapse. As a result, weighting functions have been proposed that consider calcium dynamics (Cummings et al., 1996; Shouval et al., 2002; Yeung et al., 2004), and these functions have been used to connect biophysical features of NMDARs to synaptic weight updates in models of spike time-dependent plasticity (STDP; Sejnowski, 1977; Song et al., 2000; Izhikevich, 2007; Bush and Jin, 2012; Standage et al., 2014). We now propose that the calcium transient is explicitly a function of the spine volume-to-surface area, ultrastructure, and buffers, and that the calcium functions that inform the synaptic weight vectors and synaptic learning rate (Fig. 9) must be updated to consider such geometric information. We anticipate that such new models can give us better insight into the neural circuitry of the brain and also better inform bioinspired engineering of neuromorphic circuits (Cruz-Albrecht et al., 2012). We also acknowledge that much work remains to be done in connecting the spatial signaling aspects in postsynaptic spines with neural circuit behavior but hope that this work will inspire more multiscale modeling efforts in this field. The spatial aspects of calcium dynamics are also fundamental toward understanding the downstream dynamics of critical molecules such as calcium/calmodulin-dependent protein kinase II (CaMKII), the small RhoGTPases (Cdc42, Rho, and Rac), and subsequently actin dynamics in dendritic spine remodeling (Oertner and Matus, 2005; Murakoshi and Yasuda, 2012; Miermans et al., 2017; Ohadi and Rangamani, 2019 ,Preprint; Ohadi et al., 2019 ,Preprint; Yasuda et al., 2003; Rangamani et al., 2016; Yasuda, 2017). Going beyond single spine dynamics, the propagation of the downstream mechanochemical activity to neighboring spines is a key step toward integrating single spine behavior to multiple spines, across the dendrite, and ultimately the whole cell (Majewska et al., 2000; Bloodgood and Sabatini, 2005; Herz et al., 2006; Schmidt et al., 2007; Yasuda, 2017). Thus, we posit that accounting for the spatial and physical aspects of calcium dynamics is the first step toward deciphering the complex shape–function relationships of dendritic spines.

We thank various members of the Rangamani Lab for valuable discussion, Mariam Ordyan for manuscript review, and Professor Brenda Bloodgood and Evan Campbell for their expertise. We thank Professor Pietro De Camilli and his laboratory for supplying images used to construct the mesh seen in Fig. 5 d.

This work was supported by the University of California, San Diego Interfaces Graduate Training Program and a San Diego Fellowship, Air Force Office of Scientific Research (AFOSR) Multidisciplinary University Research Initiative (MURI) grant FA9550-18-1-0051 (M. Bell, P. Rangamani, and T. Sejnowski), and the National Institutes of Health (grant T32EB009380).

The authors declare no competing financial interests.

Author contributions: M. Bell and P. Rangamani conceived the project. M. Bell conducted all the simulations. All authors conducted data analysis and contributed to writing the manuscript. M. Bell made all the figures.

José D. Faraldo-Gómez served as editor.