Dendritic spine morphology regulates calcium-dependent synaptic weight change

Dendritic spines are small protrusions along the dendrites of neurons that compartmentalize postsynaptic biochemical, electrical, and mechanical responses. These subcompartments house the majority of excitatory synapses and are key for neuronal communication and function (Nishiyama and Yasuda, 2015; Murakoshi and Yasuda, 2012). Because of their unique biochemical compartmentation capabilities, spines are thought of as computational units that can modify their synaptic strength through a process called synaptic plasticity (Yuste and Denk, 1995; Nishiyama and Yasuda, 2015). There have been many different approaches, both experimental and computational, to understand how these small subcompartments of excitatory neurons can regulate learning and memory formation (Lee et al., 2012; Segal, 2010). These studies have helped identify a few key scientific threads: (1) the biochemical signal transduction cascades in spines span multiple time scales, but Ca²⁺ is the critical initiator of these events; (2) spines have distinct morphological features that can be categorized depending on physiological or pathological conditions (Table 1; Kasai et al., 2003); and (3) synaptic weight update is a measure of the strength of a synapse and represents the strength of the connection between neurons. Synaptic weight represents changes to synaptic connection strength that occurs during synaptic plasticity, such as during long-term potentiation (LTP) and long-term depression (LTD; Lee et al., 2012). In this work, we focus on bridging these different ideas by asking the following question: how does spine morphology affect synaptic weight update? To answer this question, we develop a computational model that focuses on the stochastic dynamics of Ca²⁺ in spines of different geometries and map the synaptic weight update to geometric parameters.

Calcium plays a key role as a second messenger in biochemical and physical modifications during synaptic plasticity, triggering downstream signaling cascades within dendritic spines and the entire neuron (Yuste and Denk, 1995; Yuste et al., 2000; Augustine et al., 2003). Theoretical efforts have linked calcium levels to synaptic plasticity change through a parameter called synaptic weight that essentially captures the strength of the synapse (Jahr and Stevens, 1993; Shouval et al., 2002; Graupner and Brunel, 2012; O’Donnell et al., 2011). An increase in synaptic weight is associated with synapse strengthening, while a decrease in synaptic weight is associated with synapse weakening (Lisman, 2017; Earnshaw and Bressloff, 2006). While changes in synaptic strength require a host of downstream signaling and mechanical interactions (Xia and Storm, 2005; Mäki-Marttunen et al., 2020), the level of calcium is a well-accepted indicator of synaptic plasticity and weight (Shouval et al., 2002; Jędrzejewska-Szmek et al., 2017). This led to the hypothesis that synaptic plasticity outcome could be determined from the calcium dynamics alone; this theory has been readily used for numerous models in computational neuroscience (Shouval et al., 2002; Mahajan and Nadkarni, 2019). Because of their probabilistic nature and discrete number, calcium ion channels and receptors appear to behave stochastically (Anwar et al., 2013; Dudman and Nolan, 2009; Faisal et al., 2005). This indicates that calcium dynamics in the spine leans toward stochasticity, and it has been suggested that synaptic plasticity itself relies on stochasticity for robustness (Cannon et al, 2010; O’Donnell et al, 2011; Anwar et al, 2013; Koumura et al, 2014; O’Donnell and Nolan, 2014; Fujii et al, 2017; Tottori et al, 2019). In this work, we seek to understand how spine morphology can modulate synaptic weight update predicted through stochastic calcium dynamics.

Dendritic spines have characteristic sizes and shapes that dynamically change over time in response to stimulus and are associated with their function and synaptic plasticity (Bourne and Harris, 2008; Holthoff et al., 2002). Just as whole-cell shape is known to influence signaling dynamics (Calizo et al., 2020; Neves et al., 2008; Rangamani et al., 2013; Héja et al., 2021; Bell and Rangamani, 2021; Scott et al., 2021), studies have specifically probed the interplay between calcium dynamics and dendritic spine morphology (Bell et al., 2019; Cugno et al., 2019; Yuste et al., 2000; Bartol et al., 2015). Because of the historical significance of dendritic spines as electrical subcompartments, the morphology of the spine neck has been implicated in regulating calcium signaling, and longer spine necks were found to decouple spine-dendrite calcium signaling (Volfovsky et al., 1999). Additional modeling work coupled actin–myosin contractions to cytoplasmic flow to identify two time scales of calcium motion, driven by flow and diffusion, respectively, that depend on spine geometry (Holcman et al., 2004). A combined analytical and numerical study showed how geometry and curvature gives rise to pseudo-harmonic functions that can predict the locations of maximum and minimum calcium concentration (Cugno et al., 2019). More recently, we used a deterministic reaction-diffusion model to investigate dendritic spine morphology and ultrastructure and found that dendritic spine volume-to-surface-area ratios and the presence of spine apparatus (SpApp) modulate calcium levels (Bell et al., 2019). As we have shown before, the natural length scale that emerges for reaction-diffusion systems with boundary conditions that have influx and efflux rates is the volume-to-surface-area ratio (Calizo et al., 2020; Cugno et al., 2019). What remains unclear is whether the trends from dimensional analysis of deterministic models continue to hold despite the stochastic nature of calcium influx and efflux across the wide range of spine shapes.

In this work, using idealized and realistic spine geometries, we investigate the impact of shape and stochasticity on calcium dynamics and synaptic weight change. We seek to answer the following question: How do specific geometric parameters—namely, shape and size of dendritic spines—influence calcium dynamics and therefore synaptic weight change? To address this question, we built a spatial, stochastic model of calcium dynamics in different dendritic spine geometries. We used idealized geometries informed by the literature to control for the different geometric parameters and then extended our calculations to realistic geometries. We probed the influence of spine shape, volume, and volume-to-surface-area ratio on calcium influx, variance of calcium dynamics, and the robustness of synaptic weight. We show that although calcium dynamics in individual spines are stochastic, synaptic weight changes proportionally with the volume-to-surface-area ratio of the spines, suggesting that there may exist deterministic relationships between spine morphology and strengthening of synapses.

Ca²⁺ dynamics in dendritic spines have been previously studied using computational models (Bartol et al., 2015; Friedhoff et al., 2021; Cugno et al., 2019; Bell et al., 2019; Holcman et al., 2005). In this work, we focused on modeling effort on the early, rapid influx on Ca²⁺ for spines of different sizes and shapes with the goal of identifying relationships between spine geometry and early synaptic weight change. Our model is based on previous works (Bartol et al., 2015; Bell et al., 2019; Mahajan and Nadkarni, 2019) with some modifications and simplifications to enable us to identify the relationship between spine morphology and synaptic weight change. Inspired by Bartol et al. (2015), we converted a previous deterministic model of calcium influx (Bell et al., 2019) to a spatial, particle-based stochastic model constructed in Monte Carlo Cell (MCell; Stiles and Bartol, 2001; Stiles et al., 1996; Kerr et al., 2008), to capture the stochastic nature of Ca²⁺ dynamics in the small spine volumes. We specifically focus on dendritic spine geometries and calcium dynamics representative of hippocampal pyramidal neurons (Bell et al., 2019). We detail the steps below.

Here we list the main assumptions in the model and describe the components of the model shown in Fig. 1.

 ∙ Geometries: We investigate how spine geometry (spine volume, shape, and neck geometry) and ultrastructure (spine apparatus) can influence synaptic weight change, with the goal of drawing relationships between these different morphological features and synaptic weight (Fig. 1 e).

 ∙ Idealized geometries: Idealized geometries of thin, mushroom-shaped, and filopodia-shaped spines were selected from Alimohamadi et al. (2021), and the different geometric parameters are given in Tables S1 and S2 and Fig. 1 b.

– Postsynaptic density (PSD): For each control geometry, the PSD area was set as a fixed proportion of the spine volume.

– Size variations: For each spine geometry, we varied the volume of the control geometry to consider the impact of different morphological features.

– Spine apparatus: A SpApp is included in the thin and mushroom idealized spines by scaling the spine geometry to a smaller size and including it within the plasma membrane (PM) geometry. These variations were included to modify the volume of the spine in the presence of these organelles.

 ∙ Realistic geometries: We also investigated how spine morphology affected synaptic weight change in realistic morphologies. Realistic spine morphologies were reconstructed from 3-D EM images (Wu et al., 2017) at sufficient mesh quality to import into MCell (Lee and Laughlin, 2020). Realistic spines were selected to have a variety of morphologies to reflect filopodia-shaped, thin, and mushroom spines. PSDs were denoted based on the segmentation of the original 3-D electron micrographs.

 ∙ Time scales: Our goal is to consider the initial changes in synaptic weight due to a single calcium pulse, consistent with prior studies (Bartol et al., 2015; Shouval et al., 2002). Our focus is on early time scale events associated with synaptic weight, rather than the induction of LTP/LTD specifically. The time scale of calcium transients is rapid, on the millisecond time scale (Bartol et al., 2015; Holcman et al., 2004; Holcman et al., 2005), owing to the single activation pulse, various buffering components, and the SpApp acting as a calcium sink, rather than a source (Basnayake et al., 2019; Basnayake et al., 2021). Each spine geometry is initiated with a basal concentration of calcium as shown in Table 2. For each geometry, these concentrations were converted to numbers of Ca²⁺ ions to initialize the particle-based simulations.

 ∙ Calcium model stimulus: The stimulus used in the model is an excitatory postsynaptic potential (EPSP) and back-propagating action potential (BPAP) offset by 10 ms and a glutamate release that activates N-methyl-D-aspartate receptor (NMDAR; Bartol et al., 2015) as shown in Fig. 1 a, inset. We include the presynapse as a surface from which glutamate is released from a central location. α-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptor (AMPARs), which competes with NMDARs, to bind glutamate, are also included in the model but do not contribute to the calcium influx.

 ∙ Channel and receptor dynamics: We assume that the surface density of the receptors and channels on the membrane of the spine is constant and uniformly distributed (Sabatini and Svoboda, 2000; Holthoff et al., 2002). This assumption is based on experimental observations (Sabatini and Svoboda, 2000) and has been used in other computational models of calcium dynamics (Bell et al., 2019; Keller et al., 2008; Franks et al., 2002). An important consequence of this assumption is that when the surface area of the spine changes, the total number of receptors will also change. How calcium influx scales with spine volume is an important consideration with implications on how calcium concentration scales with spine size (O’Donnell et al., 2011). The constant receptor density assumption means that calcium influx undercompensates for increases in spine volume (O’Donnell et al., 2011).

 ∙ Boundary conditions: Calcium ion influx occurs through NMDAR localized to the PSD region and voltage-sensitive calcium channels (VSCCs) on the PM, based on Bartol et al. (2015). Calcium ions leave the spine volume through the pumps on the PM, plasma membrane Ca²⁺-ATPase (PMCA), and sodium–calcium exchanger (NCX), and into the SpApp, when present, through sarco/endoplasmic reticulum Ca²⁺-ATPase (SERCA) (Bell et al., 2019). We consider the spine as an isolated geometric compartment and do not consider the effect of calcium influx from the dendrite at this time scale. The base of the spine neck has a Dirichlet boundary condition of calcium clamped to zero and acts as a calcium sink, which represents Ca²⁺ leaving the spine into the dendrite owing to the sudden increase in calcium in the spine (Holcman et al., 2005).

 ∙ Buffers: We do not model the different buffer species but rather use a lumped parameter approach as was done before (Mahajan and Nadkarni, 2019; Bell et al., 2019). Because free Ca²⁺ is rapidly buffered in cells, we consider both mobile buffers in the cytoplasm and immobile buffers on the PM (Bell et al., 2019; Schmidt and Eilers, 2009; Schmidt, 2012). We also include an exponential decay of calcium throughout the cytoplasm to capture the complex cytosolic buffering dynamics without including explicit buffers. Addition of more species introduces many more free parameters and can make the model computationally intractable; therefore, we focus on a lumped parameter approach.

 ∙ Stochastic trials: Each simulation condition was run with 50 random seeds, and these individual runs were averaged to obtain mean and SD (Friedhoff et al., 2021; Bartol et al., 2015).

 ∙ Model readouts: We report Ca²⁺ dynamics in terms of the number of ions rather than concentration. This is because the total number of ions in the spine reflects total signal coming into the spine and is the natural output from these particle-based simulations. The total number of calcium ions is used as input to calculate the synaptic weight change.

 ∙ Synaptic weight: The synaptic plasticity model developed by Shouval et al. (2002) was adapted to have dependence on total calcium ions rather than calcium concentration. The rate of synaptic weight update depends on a learning rate, τ_w, and a thresholding function, Ω_w, that are both dependent on calcium ion levels (Fig. 1, c and d). The learning rate determines the rate of synaptic weight change, while Ω_w determines if the weight increases or decreases. Thresholds for LTP and LTD, θ_P and θ_D, are set so that an intermediate level of calcium leads to a weakening of a synapse and LTD, while an elevated level of calcium leads to the strengthening of a synapse and LTP (Shouval et al., 2002; Cho et al., 2001; Cormier et al., 2001).

We summarize the main reactions for Ca²⁺ in the volume. These reaction models were obtained from Bartol et al. (2015) and Bell et al. (2019) and are discussed in detail below. Model parameters are given in Table 2. We found that our calcium dynamics are comparable to previously published models (Bell et al., 2019; Bartol et al., 2015; Rubin et al., 2005; Hu et al., 2018) and experimental observations (Sabatini et al., 2002; Hoogland and Saggau, 2004; Segal and Korkotian, 2014; Fig. S8).

Reaction rates for mobile and fixed buffers are provided in Table 2.

The primary influx of calcium through the PM occurs through NMDARs and VSCCs, and calcium is pumped out of the volume through PMCA and NCX. In this model, NMDARs depend on both voltage and glutamate and are localized to the PSD region. VSCCs are voltage dependent and localized throughout the PM. PMCA and NCX are calcium-dependent pumps and are also located throughout the PM surface.

NMDARs are localized to the PSD area with a surface density of 150 molecule µm⁻² (Bartol et al., 2015). The activation of NMDAR is modeled with an asymmetric trapping block kinetic scheme as proposed in Vargas-Caballero and Robinson (2004). The activation of NMDAR depends on the diffusion of glutamate through the synaptic cleft and its binding to inactive receptors. A surface identical to the top of the spine head is translated 2 µm above the head to represent the synaptic cleft. At time t = 0 in each simulation, 500 molecules of glutamate are released at the center of this synaptic cleft. The released glutamate molecules diffuse through the cleft at a rate of 2.2 × 10⁻⁶ cm² s⁻¹ and bind to membrane-bound proteins. On the postsynaptic membrane, NMDARs compete with the glutamate receptor AMPAR for glutamate; thus, AMPARs are also included in the simulation to model this competition, but do not play a role in calcium influx. AMPAR is also localized to the PSD area at a density of 1,200 molecule µm⁻² (Bartol et al., 2015). The binding of glutamate to AMPAR is modeled according to the kinetic scheme proposed in Jonas et al. (1993).

The influx of Ca²⁺ through VSCCs is also dependent on the activation kinetics of VSCCs. The initial conditions for all the VSCCs is the closed state, and the activation of the channels is modeled here with a five-state kinetic scheme as used in Bartol et al. (2015). The parameters for Ca²⁺ influx through VSCCs are the same as in Bartol et al. (2015). VSCCs are located on the PM with a density of 2 molecules µm⁻².

Because the transmembrane potential varies with time (Fig. 1 a, inset) and the rate constants for NMDAR and VSCC are voltage dependent, the values of these rate constants at each simulation step were precomputed and passed into MCell. The voltage stimulus representing a single EPSP starting at time t = 0, followed by a single BPAP occurring at an offset of 10 ms, was obtained from Bartol et al. (2015). Note that this time offset is within the typical window for spike timing–dependent plasticity (STDP) to induce LTP (Bartol et al., 2015; Griffith et al., 2016).

PMCA and NCX are located on the PM with areal density 998 and 142 molecule µm⁻², respectively (Bartol et al., 2015), forcing calcium to flow out of the cell. These pumps are modeled using the set of elementary reactions and reaction rates from Bartol et al. (2015).

Extracellular calcium was not explicitly modeled for ease of computational tractability. We assumed a constant extracellular calcium concentration (2 mM) that is negligibly impacted by calcium influx to and efflux from the spine cytoplasm. The dynamics of Ca²⁺ are explicitly modeled once they enter the cell through channels located on the PM and cease to be explicitly represented once they are pumped out of the cell.

Synaptic weight update was calculated using the classic model from Shouval et al. (2002). The governing equations were modified to take total number of Ca²⁺ ions rather than a concentration. Here, we use total number of ions to highlight the details available from a stochastic simulation and to consider the consequences of using this global readout on synaptic weight.

Calcium simulations were conducted for a total simulation time of 35 ms with a 500 ns time step. Each geometry is simulated in MCell with 50 distinct seeds to generate an appropriate sample size of results. All simulations use a write-out frequency of once per iteration for reproducibility of results. Longer write-out frequencies introduce nondeterminism to the trajectories arising from the MCell reaction scheduler. At the beginning of each simulation, membrane proteins are randomly distributed over specified regions of the spine geometry surface area according to an assigned count or concentration. System configuration and analysis scripts are all available on Github (https://github.com/RangamaniLabUCSD/StochasticSpineSimulations).

Fig. S1 presents the profiles of the idealized geometries used and various geometric parameters. Fig. S2 provides an artificial calcium input to demonstrate how the synaptic weight function depends on calcium transient dynamics. Figs. S3 and S4 provide additional plots of calcium transient dynamics against volume and as concentration in temporal plots, respectively. Figs. S5, S6, and S7 provide additional plots of the thin spine neck variations, mushroom spine neck variations, and thin spine with SpApp variations, respectively. Figs. S8 and S9 also compare our results to previous calcium transients and calcium dynamics, respectively. Figs. S10 and S11 compare the use of total ion number versus concentration to determine synaptic weight. Fig. S12 provides a matrix of two-tailed t tests between all simulations, and Fig. S13 provides synaptic weight predictions for a pulse train. Videos 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 show the spatiotemporal calcium dynamics. Table S1 and Table S2 list all geometric variations. Table S3 lists values for realistic geometries. Details on how to use realistic geometries in these simulation modalities can be found in the supplemental text at the end of the PDF.

We first demonstrate the coupling between stochastic dynamics of calcium in a single spine to the deterministic model of synaptic weight update (Fig. 2). We then investigate whether spine size has any effect on synaptic weight change of filopodia-shaped spines (Fig. 3), thin spines (Fig. 4), and mushroom-shaped spines (Fig. 5). Next, we consider the role of the SpApp (Fig. 6). Finally, we investigate the relationship between spine morphology and synaptic weight update in realistic geometries (Fig. 7). Our results predict that synaptic weight change through calcium dynamics is a deterministic function of geometric parameters of the spines (Fig. 8). We discuss these results in detail below.

We demonstrate how the stochastic calcium model informs the deterministic synaptic weight predictions. We consider a single-seed trial for a single geometry, in this case seed 1 from a realistic thin spine (Fig. 2 a). We observed that in response to the voltage trace and glutamate release (Fig. 2 b), the NMDARs and VSCCs stochastically opened and closed (Fig. 2, c and d), resulting in a noisy calcium transient (Fig. 2 e). We compared our calcium transients to both published experimental (Sabatini et al., 2002; Hoogland and Saggau, 2004; Segal and Korkotian, 2014) and computational (Bell et al., 2019; Bartol et al., 2015; Rubin et al., 2005; Hu et al., 2018) results, and found a reasonable agreement (Fig. S8). We next calculated the learning rate (Fig. 2 f) and Ω_w (Fig. 2 g) terms in the synaptic weight model as a function of this calcium transient. We noticed that the noisy calcium dynamics were filtered into a smoother synaptic weight prediction (Fig. 2 h). For another example of how calcium pulse magnitude and width translates to synaptic weight update, see Fig. S2. With this understanding of how the two models integrate, we next investigated how spine geometry influences calcium transients and subsequent synaptic weight predictions.

We began our analysis with a simple question: does spine size alter synaptic weight change? To answer this question, we first examined filopodia-shaped spines. Dendritic filopodia are precursors of dendritic spines and serve to bridge the gap between the dendrite and an axon that is passing by during synapse formation (Ozcan, 2017). These are highly motile, elongated structures that resemble tubules (lengths of 2–20 µm and neck diameters <0.3 µm). The simplicity of this geometry allows us to focus on the role of size alone in a simple spine geometry. We used spine geometries of three different volumes (0.017, 0.058, and 0.138 µm³). Simulations revealed that the calcium dynamics in these tubule-shaped spines appeared to follow a “plug-flow” behavior in which, at 15 ms, all the calcium was localized to one region (Fig. 3 a). This behavior is because of the narrow geometry of the spine, preventing dispersion of the calcium (see also Video 1). Next, we examined the temporal dynamics of calcium and noted that the larger spines had larger numbers of calcium ions (Fig. 3 b) but also a larger variance of calcium ions (Fig. 3 c). We further characterized the dynamics by considering the peak calcium values and decay time constants of the calcium transients versus the spine volume-to-surface-area ratio. We chose the volume-to-surface-area ratio as a geometric metric of spine morphology because it encompasses both the cytosolic volume through which calcium diffuses and the surface area of the spine membrane through which calcium can enter and leave the system. Additional analyses with respect to spine volume are shown in Fig. S3 and example calcium transients are shown in Fig. S4.

We noted that, indeed, increasing spine size, and therefore the volume-to-surface-area ratio, causes a linearly proportional and significant increase in peak calcium ions (Fig. 3 d). We also found that the decay time of calcium from the peak decreased with increasing volume-to-surface-area ratios and satisfied an exponential dependence (Fig. 3 e). As spine size increases, the decay time constant decreases, showing that it takes longer for calcium to clear out of larger spines and spines with larger volume-to-surface-area ratios. Finally, we calculated the synaptic weight change (see Synaptic weight change) and compared this value at 35 ms across volume-to-surface-area ratios for the filopodia-shaped spines (Fig. 3 f). We observed that while the smallest spine had no observable weight change, presumably because of the net low calcium influx, the weight change increases with increase in spine volume-to-surface-area ratio (Fig. 3 f). Thus, we found that even for a shape as simple as a filopodia-shaped spine, changes in spine volume-to-surface-area ratio can dramatically alter calcium dynamics and synaptic weight change in stochastic conditions, suggesting a close coupling between spinogenesis and calcium handling.

We next asked whether the relationships of spine size and synaptic weight change observed for filopodia-shaped spines (Fig. 3) also held for thin and mushroom-shaped spines. Thin and mushroom-shaped spines emerge from filopodia-shaped spines as spinogenesis progresses (Ozcan, 2017; Freire, 2010). It has been proposed that spines exist in a continuum of shapes (Ofer et al., 2021), but historically it has been useful to categorize spines into specific categories of shapes (Yuste, 2010). Thin spines, with small heads and thin necks, have been classified as “write-enabled” or learning spines due to their high motility. Mushroom spines, on the other hand, with bulbous heads and relatively wider necks, are termed “write-protected” or memory spines due to their stability (Jasinska et al., 2019). Thin spines are characterized by a spherical head, and we repeated the calcium influx simulations in a small thin spine with three different spine neck lengths (0.07, 0.06, and 0.04 µm) and thin spines of three different volumes (0.035, 0.119, and 0.283 µm³) that were informed by the ranges found in the literature (Fig. 4). We observed that, in thin spines, the calcium ions were concentrated in the head at 15 ms but dispersed more uniformly by 30 ms (Fig. 4, a and b; and Video 2). We did not observe plug–flow-like behavior as we did for filopodia-shaped spines, likely because of the differences in both shape and volume of the thin spines. The thin spines with different neck lengths showed very similar calcium transients and variance (Fig. 4, c and d), except for the thin-necked spine, which showed much more variance during its decay dynamics. For a closer look at the thin spine neck variation dynamics, see Fig. S5.

Calcium dynamics in thin spines follows the expected temporal dynamics (Fig. 4 e), with larger spines having larger peak calcium and increased time to decay. Larger thin spines also have larger variance in calcium ion transients over time (Fig. 4 f). Next, we found that the maximum calcium ions per spine was significantly larger in larger spines, with statistically different values for the different-sized spines. The peak calcium increased linearly compared with spine volume-to-surface area but with a smaller slope compared with the filopodia-shaped spines (maximum peak values in filopodia-shaped spines increased three times faster than those in thin spines; Fig. 4 g). This suggests that the size dependence of calcium grows slower in thin spines than in filopodia-shaped spines. The decay time also showed an exponential decay in thin spines with increasing volume-to-surface-area ratio (Fig. 4 h). The exponent was smaller for thin spines compared with filopodia-shaped spines (47.9 versus 22.43), suggesting that the decay rate with respect to volume-to-surface-area ratio was slower in thin spines. Finally, the synaptic weight change showed an increase with volume-to-surface-area ratio in thin spines (Fig. 4 i), indicating that larger spines are capable of stronger learning outcomes.

Finally, we repeated our analysis for mushroom-shaped spines of different neck length (0.13, 0.10, and 0.08 µm) and increasing volume (0.080, 0.271, and 0.643 µm³; Fig. 5). The effect of the shape of the spines is evident in the spatial dynamics of calcium (Fig. 5, a and b; and Video 3). Even at 15 ms, we noted that while a vast majority of calcium ions were localized in the spine head, there was spillover of calcium into the neck; this is particularly evident in the spines of larger volume (Fig. 5, a and b). Spine neck length showed similar increase and decay dynamics to each other (Fig. 5, c and d), but the thick-neck mushroom spine in particular showed a reduced variance. For a closer look at the mushroom spine neck variation dynamics, see Fig. S6. The effect of increases in volume, and therefore increases in volume-to-surface-area ratio, on the temporal dynamics of calcium is an increase in peak calcium (Fig. 5, e and g) and variance (Fig. 5 f) and a decrease in the decay time constant (Fig. 5 h). Interestingly, changing the neck length in the mushroom spine creates the opposite trend with respect to volume-to-surface-area ratio (Fig. 5, g and h), suggesting that spine neck length is an additional geometric factor for mushroom spines to tune their calcium response. The synaptic weight change in mushroom spines increases with spine volume-to-surface-area ratio and is larger for these mushroom spines than the filopodia-shaped and thin spines (Fig. 5 i). We observed that the peak calcium showed a linear increase with volume-to-surface-area ratio, with a slope that lies between the thin spines and filopodia-shaped spines. Finally, the decay time constant decreased with spine volume-to-surface-area ratio but with a smaller exponential decay compared with thin spines and filopodia-shaped spines. These two results point to the following conclusions. First, an increase in spine volume results in an increase in critical readouts of synaptic plasticity, and second, the shape of the spine alters the quantitative relationships of synaptic plasticity by allowing access to different volume-to-surface-area ratios.

Approximately 14% of dendritic spines have specialized ER called SpApp, which are preferentially present in larger, mature spines (Chirillo et al., 2019; Spacek and Harris, 1997; Bell et al., 2019). Furthermore, recent studies have shown that the SpApp and the ER are dynamic structures in the dendrite and dendritic spines (Perez-Alvarez et al., 2020). Previously, we showed that the SpApp modulates calcium transients in deterministic models of calcium influx (Bell et al., 2019) by altering the net fluxes (Cugno et al., 2019). Here, we investigated how these relationships are altered in stochastic models in thin and mushroom spines (Fig. 6). When a SpApp is present in the spine head, it effectively reduces the volume of the spine cytosol and, in the time frame of our consideration, acts as a calcium sink (by the action of the SERCA pumps; Sabatini et al., 2002). One example trajectory in a mushroom spine with a SpApp is shown in Video 4. We varied SpApp size in the small-sized thin spine and medium-sized mushroom spine (Fig. 6, a and b; and Table S2). Calcium transients and variance showed much smoother dynamics for the mushroom spines than the thin spines (compare Fig. 6, e versus c). Peak calcium values were all statistically different for the different SpApp sizes in the mushroom spines but not the thin spines. Decay time constants were fitted with an exponential relationship (Fig. 6 h), but there were no statistical differences across different mushroom spines. All different sizes of the SpApp produce synaptic weight changes that are statistically different in the mushroom spines; increases in SpApp size result in smaller spine volume (and smaller volume-to-surface-area ratio) and therefore produce smaller weight changes (Fig. 6 i). The thin spines had a more complex trend and did not have statistically significant differences in predicted synaptic weight. For a closer look at the variations of the SpApp within thin spines, see Fig. S7. We conclude that the presence of SpApp alters the volume-to-surface-area ratio for spines and therefore tunes calcium levels and synaptic weight updates in the large mushroom spines with an inverse relationship to SpApp size.

Thus far, we focused on idealized geometries of spines to identify relationships between key synaptic variables and key geometric variables. We found that the peak calcium value, decay time constant, and synaptic weight depend on the volume-to-surface-area ratio within each shape classification. Do these relationships hold for realistic geometries as well? To answer this question, we selected realistic geometries from mesh models (Lee and Laughlin, 2020) informed by electron micrographs from Wu et al. (2017).

Realistic spines have more complex geometries that do not fall into the exact morphological categories that we used for idealized spines. To test the significance of these variations, we selected two spines of each shape (thin, mushroom, and filopodia) and conducted simulations with the exact same parameters as the idealized simulations (Fig. 7 a). We chose realistic geometries that were within the range of sizes of the idealized geometries. The PSDs in the realistic spines were annotated during the segmentation process, and no modifications were made to the PSD marked regions. To capture filopodia-shaped protrusions, we selected long, thin spines (with minimal differentiation between the head and neck) that had marked PSD, because we did not include dendritic filopodia in the dendrite section. Details on how to use realistic geometries in these simulation modalities can be found in the supplemental text at the end of the PDF. We showed the spatial distribution of calcium ions for a single seed for filopodia, thin, and mushroom spines (Fig. 7 b) and found that due to the complexity of realistic morphologies, the calcium distribution was more complicated than those observed in the idealized spines.

For filopodia-shaped spines, we found that peak calcium and variance varied with volume, but the variance was not appreciably different for the two spines that we used to conduct simulations (Fig. 7, c i and ii; and Videos 5 and 6). The realistic thin spines we chose had volumes similar to the filopodia-shaped spines, and they also exhibited calcium dynamics proportional to their volume (Fig. 7, c iii and iv; and Videos 7 and 8). Mushroom spines had larger volumes and larger PSD areas compared with thin or filopodia spines (Fig. 7, c v and vi; and Videos 9 and 10). Again, the calcium dynamics was proportional to the volume and showed that larger spines have higher peak calcium values. Thus, the relationships of spine geometry and calcium dynamics hold in realistic geometries as well.

Dendritic spines have been extensively studied as biochemical signaling compartments, and their role in calcium sequestration has been theorized (Bell et al., 2019; Cugno et al., 2019; Yuste et al., 2000; Kotaleski and Blackwell, 2010; Murakoshi and Yasuda, 2012; Yasuda, 2017; Friedhoff et al., 2021). Their unique morphological traits and the classification of spine sizes and shapes with respect to function suggest possible structure–function relationships at the level of individual spines. In this work, we used stochastic modeling of calcium transients in dendritic spines of different geometries to understand how spine size and shape affect the change in synaptic weight. Using a stochastic simulation is important to investigate variance among spine shape and size, as dendritic spines have small volumes and probabilistic channel dynamics. Using idealized and realistic geometries, we found that geometric properties, specifically volume-to-surface-area ratio, affected key properties of calcium transients including peak calcium, decay time constants, and synaptic weight change. We discuss these findings in the context of different aspects of synaptic plasticity.

Our models predict that despite the individual calcium transients being stochastic, there is a predictive deterministic trend that appears to carry through the different sizes and shapes of the spines used in our model (Fig. 8). One of the advantages of our modeling approach here is that we can directly compare across the entire range of idealized and realistic geometries. By considering all the data from our models, for a total of 18 geometries with 50 simulations in each, we found that the peak calcium number is more or less linear with the volume, surface area, and volume-to-surface-area ratio (Fig. 8, d, g, and j). The decay time constant for calcium transients shows an exponential decay for larger volume-to-surface-area ratios, volumes, and surface areas, with quite a bit of variability for smaller ratios (Fig. 8, e, h, and k). We note that both peak calcium and decay time constants show clearer trends within the same spine protrusion type (i.e., comparing within the same color). Finally, the synaptic weight change increases as volume-to-surface-area ratio, volume, and surface area increase (Fig. 8, f, i, and l). We emphasize that our goal is to demonstrate a trend in the data as opposed to building numeric functions. Although we fitted the various data, the r² is often weak, indicative of the complexities that underlie such efforts.

We want to highlight two takeaways from the synaptic weight trends with respect to volume-to-surface-area ratio. The first takeaway is that within a spine shape group (comparing within a specific color in Fig. 8), there are clear increasing trends with respect to volume-to-surface-area ratio. The second takeaway is that while there are general trends in the data highlighted by the fit lines in Fig. 8, there appear to be three regions with slightly different synaptic weight trends at small, intermediate, and large volume-to-surface-area ratios. We discuss the possible consequences of these trends in more detail below.

In the idealized geometries, the PSD area is a manually fixed proportion of the spine volume, but realistic geometries do not have this artificial constraint. Therefore, we redid our analysis using volume to PSD area ratio and PSD area-to-surface-area ratios (PSD to PM ratio). Interestingly, we did not see a clear trend within the plots against volume to PSD area ratio (Fig. 8, a–c). In comparison, the PSD area to PM area ratio shows the same relationships overall as volume-to-surface-area ratio (Fig. 8, m–o), but this time with clustering of data around some ratios. This indicates that the PSD area is an important additional degree of freedom for synaptic weight change that must be considered for interpretation of geometric features, and using realistic geometries with boundary markings allows us to investigate this. It is important to note that there is a lot more variability in the smaller volume-to-surface-area ratios, suggesting that the response of smaller spines may be more variable than that of larger spines. This feature can work as a double-edged sword: it may provide an advantage during the development of spines or be an disadvantage in the case of loss of spines (Yuste and Bonhoeffer, 2004; Stein and Zito, 2019).

We interpret our predictions in the context of spine shapes. Filopodia are prevalent during early synaptogenesis and can transition into dendritic spines based on synaptic activity (Ozcan, 2017). Additionally, various disease states produce modified dendritic spines that appear more like filopodia (Ruhl et al., 2019). The lack of significant weight changes for the smallest filopodia-shaped spine indicates that there is a volume threshold at which filopodia receive enough stimulus to trigger synaptic weight change and transition toward more stable, mature dendritic spines. Importantly, the early synaptic weight changes emphasize how the increase in spine volume changes the weight outcome from LTD to LTP. This increase in synaptic weight emphasizes how an increase in spine size can push a thin spine to transition into a stable, larger mushroom spine.

The difference in peak calcium level, decay dynamics, and synaptic weight changes as different spine shapes are scanned across different sizes can also provide insight into spine shape transitions during development and maturation. Filopodia-shaped spines have larger increases in peak calcium levels and synaptic weight updates and faster decreases in decay time constants as their volume-to-surface-area ratios and volumes increase, compared with both thin and mushroom spines (Figs. 3, 4, and 5). This suggests that filopodia spines can very quickly alter their calcium levels and therefore are well suited for initially identifying possible synaptic partners and subsequently directing resources to those filopodia that are good candidates to transition to dendritic spines (Lohmann and Bonhoeffer, 2008). Once filopodia are established, their linear calcium increase with volume might be unsustainable and might lead to reduced levels of increase for thin spines of comparable volume-to-surface area (and volume). This suggests that larger stimuli might be necessary to push thin spines toward more excitation, perhaps to prevent excessive numbers of thin spines from maturing and leading to resource depletion and excess neural connectivity (Sorra and Harris, 2000). Mushroom spines once again show more of an increase in synaptic weight as they increase in volume-to-surface-area ratio (and volume), but at volumes shifted from the filopodia-shaped spines, perhaps highlighting their role as key communication hubs (Sorra and Harris, 2000). The volume shift seen in mushroom spines versus filopodia-shaped spines might serve to limit the number of mature, highly excitable dendritic spines as both a key neuronal network and resource regulation feature. When the SpApp acts as a sink, its presence dampens synaptic weight changes in mushroom spines, potentially acting to stabilize the spine from future changes, as suggested by others (Jasinska et al., 2019; Mahajan and Nadkarni, 2019).

When considering why these trends hold across volume-to-surface-area ratios, it is important to note that Ca²⁺ influx is through receptors and channels with constant densities at the PM, or in the case of NMDARs, localized to the PSD. Therefore, as spines get larger, they have more surface area and more Ca²⁺ influx, which leads to higher numbers of total Ca²⁺ ions. This increase in total ions due to constant receptor and channel densities explains the increasing trend in peak Ca²⁺ number. When considering decay dynamics, Ca²⁺ efflux is due to pumps of constant density on the PM or the SpApp. Additionally, Ca²⁺ ions decay everywhere in the cytoplasm, bind to mobile buffers in the cytoplasm and fixed buffers on the PM, and bind to the spine neck base, which acts as a sink. Therefore, since many efflux or binding terms are either on the PM through pumps or at the base of the spine neck, larger volume-to-surface-area ratios mean that ions must diffuse further to reach the neck base or PM, explaining why decay time constants seem to decrease with increasing volume-to-surface-area ratio.

While changing geometric features can occur when spines increase and decrease in volume, they can also modify their volume, surface area, and volume-to-surface-area ratio by having a SpApp or through changes in spine neck geometry. We investigated these additional features (Figs. 4, 5, 6, S5, S6, and S7) and found that spine neck and SpApp size had volume-dependent effects. The smaller thin spine neck length and SpApp variations did not show much influence on peak calcium, decay rate, or synaptic weight, while the larger mushroom spine neck length and SpApp variations did have some impact on these readouts. Therefore, there are various means by which a spine can modify its synaptic weight response.

There has been substantial debate on deterministic versus stochastic studies for spine signaling (Kummer et al., 2005; Rüdiger, 2014; Skupin et al., 2010). Numerous studies have looked at the importance of stochastic calcium dynamics (Friedhoff et al., 2021; Friedhoff and Ramlow, 2021), and we agree that the consideration of stochasticity is important as noise often leads to efforts to average out its effects (Anwar et al., 2013). However, comparing our findings here to our previous deterministic results (Bell et al., 2019) shows that geometric factors play a critical role in determining Ca²⁺ dynamics; both approaches show that Ca²⁺ characteristics depend on the volume-to-surface-area ratio. Additionally, our hybrid approach of stochastic calcium dynamics and deterministic synaptic weight update is becoming increasingly common (Anwar et al., 2013; Rüdiger, 2014; Rodrigues et al., 2022,Preprint). However, care should still be taken in assuming model type as the dynamics of the species, not just particle number, plays an important role in the stochasticity of the system (Kummer et al., 2005).

We note that our study is only a small piece of the puzzle with respect to synaptic plasticity. There are many open questions remaining. Of particular interest and needing additional exploration is whether one should use total number of calcium ions or calcium concentration in evaluating synaptic weight change. For instance, we found that when calcium results are converted from total ions to average concentration along with the phenomenological synaptic weight equations, we got different trends in synaptic weight update results (Figs. S10 and S11). We note that this model of synaptic weight change has been used previously for concentration studies (Mahajan and Nadkarni, 2019; Shouval et al., 2002). We also observe that converting our previous results (Bell et al., 2019) into total ions shows the same trends for maximum Ca²⁺ peak and decay time constants as the current study (Fig. S9). Thus, a simple unit consideration can lead to conflicting results in spatial models and indicates that we need further discussion and investigation on the structure of phenomenological equations for synaptic weight to understand which factors of calcium dynamics matter and to what degree. Additional investigation is also needed in experimental data to relate fluorescence readouts to concentration or molecule numbers. However, we do compare our calcium transients to previously published experimental and computational results and find reasonable agreement (Fig. S8).

Related to these conflicting findings when considering ion total versus concentration, previous studies have considered the assumptions between calcium influx and spine geometry, more specifically the assumption of how calcium influx scales with spine volume (O’Donnell et al., 2011). Here the constant receptor and channel density assumption leads to an under-compensation scenario in which calcium influx does not scale with spine volume, leading to lower calcium concentrations for larger spines. See Fig. S4 for examples of our model results in terms of calcium concentration. Because the synaptic weight model depends on calcium influx, when using concentration to determine synaptic weight, larger volume spines have less synaptic weight increase. Therefore, the sublinear calcium influx assumption leads to this discrepancy in synaptic weight predictions based on total calcium ions versus ion concentration. See Figs. S10 and S11 for examples of this discrepancy. Further research is needed to determine how exactly calcium influx scales with dendritic spine volume in vivo, as it is currently unknown which assumption is correct (O’Donnell et al., 2011). Regardless of the relationship between dendritic spine volume and calcium influx, the spine uses various means to modify its calcium transients, including internal organelles such as the SpApp acting as either a calcium source or sink (Matsuzaki, 2007). More research is needed to explore the relationship of geometry-dependent calcium trends and their consequences on phenomenological synaptic weight predictions.

An additional limitation of this study is the use of traditional P values for statistical analysis of the data (see Fig. S12 for details on h and P values), since the statistics field has suggested moving away from null-hypothesis significance testing (Wasserstein et al., 2019). We also note that our current focus is on very early events, and these models must be extended to longer-time-scale events to explore the biochemical and geometric interplay for downstream signaling (Bhalla, 2004; Ohadi and Rangamani, 2019; Ohadi et al., 2019; Mäki-Marttunen et al., 2020). Associated with these longer-time-scale events, calcium often occurs in pulse trains owing to high-frequency stimulation of the dendritic spine (Chen et al., 1999; Zhu et al., 2015). We compared synaptic weight predictions for a single calcium transient to those due to a pulse train of activation at a single frequency (Fig. S13). However, further investigation should be done to more closely consider the role of stimulus magnitude and frequency on synaptic weight update. In addition, it is important to note that this calcium model and these dendritic spine geometries are representative of hippocampal pyramidal neurons. Calcium signaling, dendritic spine structure, and synaptic weight induction are neuron type specific, and other studies, including some MCell simulations, have investigated calcium signaling and synaptic plasticity in other neuron types (Friedhoff et al., 2021; Antunes and Simoes-de-Souza, 2018; Antunes and Simoes de Souza, 2020; Koumura et al., 2014). In some neuron types, including Purkinje cells, calcium release from the ER can play a vital role in calcium dynamics and subsequent synaptic plasticity (Koumura et al., 2014); thus, care must be taken when considering different neuron types.

In summary, our computational models using idealized and realistic geometries of dendritic spines have identified potential relationships between spine geometry and synaptic weight change that emerge despite the inherent stochasticity of calcium transients. We predict that dendritic spine morphology alters calcium dynamics to achieve their characteristic functions; in particular, so that filopodia can quickly change their synaptic weight, large mushroom spines can solidify their synaptic connections, and intermediate-sized spines require more activation to achieve larger synaptic weight changes. Additionally, we predict that within a certain spine shape, increasing volume (and increasing volume-to-surface-area ratio), while assuming receptors and channels are also recruited, allows for a larger future increase in synaptic weight, suggesting that the volume change associated with LTP and LTD serves to reinforce the biochemical changes during synaptic plasticity. Therefore, spine morphology tunes synaptic response. The advances in computational modeling and techniques have set the stage for a detailed exploration of biophysical processes in dendritic spines (Miermans et al., 2017; Basnayake et al., 2019; Ohadi and Rangamani, 2019). Such efforts are critical for identifying emergent properties of systems behavior and also eliminating hypotheses that are physically infeasible (Bell and Rangamani, 2021; Lee et al., 2021). Models such as this and others can set the stage for investigating longer-time-scale events in spines, including the downstream effectors of calcium (Jędrzejewska-Szmek et al., 2017; Mäki-Marttunen et al., 2020; Hayer and Bhalla, 2005; Ordyan et al., 2020), and actin remodeling for structural plasticity (Bonilla-Quintana et al., 2020; Rangamani et al., 2014).

Jeanne M. Nerbonne served as editor.

We thank Dr. Tom Bartol for helpful discussions on the use of MCell and Drs. Lingxia Qiao and Ali Khalilimeybodi for comments and proofreading.

This work was supported by a National Defense Science and Engineering Graduate Fellowship to M.K. Bell, a Hartwell Foundation Postdoctoral Fellowship and Kavli Institute of Brain and Mind Innovative Research Grant #2021-1755 to C.T. Lee, and Air Force Office of Scientific Research FA9550-18-1-0051 to P. Rangamani. MCell development is supported by the National Institute of General Medical Sciences-funded (P41- GM103712) National Center for Multiscale Modeling of Biological Systems (MMBioS).

The authors declare no competing financial interests.

Author contributions: Conceptualization: M.K. Bell, C.T. Lee, P. Rangamani; data curation: M.V. Holst, M.K. Bell, C.T. Lee; formal analysis: M.K. Bell, M.V. Holst; funding acquisition: P. Rangamani; investigation: M.K. Bell, M.V. Holst; methodology: M.K. Bell, C.T. Lee, P. Rangamani; project administration: M.K. Bell, C.T. Lee, P. Rangamani; resources: C.T. Lee, P. Rangamani; software: M.V. Holst; supervision: M.K. Bell, C.T. Lee, P. Rangamani; validation: M.K. Bell, M.V. Holst; visualization: M.K. Bell, M.V. Holst; writing - original draft: M.K. Bell, M.V. Holst, C.T. Lee, P. Rangamani; writing - review & editing: M.K. Bell, M.V. Holst, C.T. Lee, P. Rangamani.