A Glutamatergic Spine Model to Enable Multi-Scale Modeling of Nonlinear Calcium Dynamics

Hu, Eric; Mergenthal, Adam; Bingham, Clayton S.; Song, Dong; Bouteiller, Jean-Marie; Berger, Theodore W.

doi:10.3389/fncom.2018.00058

ORIGINAL RESEARCH article

Front. Comput. Neurosci., 27 July 2018
Volume 12 - 2018 | https://doi.org/10.3389/fncom.2018.00058

A Glutamatergic Spine Model to Enable Multi-Scale Modeling of Nonlinear Calcium Dynamics

Eric Hu^*

Adam Mergenthal

Clayton S. Bingham

Dong Song

Jean-Marie Bouteiller^*

Theodore W. Berger

Department of Biomedical Engineering, University of Southern California, Los Angeles, CA, United States

In synapses, calcium is required for modulating synaptic transmission, plasticity, synaptogenesis, and synaptic pruning. The regulation of calcium dynamics within neurons involves cellular mechanisms such as synaptically activated channels and pumps, calcium buffers, and calcium sequestrating organelles. Many experimental studies tend to focus on only one or a small number of these mechanisms, as technical limitations make it difficult to observe all features at once. Computational modeling enables incorporation of many of these properties together, allowing for more complete and integrated studies. However, the scale of existing detailed models is often limited to synaptic and dendritic compartments as the computational burden rapidly increases when these models are integrated in cellular or network level simulations. In this article we present a computational model of calcium dynamics at the postsynaptic spine of a CA1 pyramidal neuron, as well as a methodology that enables its implementation in multi-scale, large-scale simulations. We first present a mechanistic model that includes individually validated models of various components involved in the regulation of calcium at the spine. We validated our mechanistic model by comparing simulated calcium levels to experimental data found in the literature. We performed additional simulations with the mechanistic model to determine how the simulated calcium activity varies with respect to presynaptic-postsynaptic stimulation intervals and spine distance from the soma. We then developed an input-output (IO) model that complements the mechanistic calcium model and provide a computationally efficient representation for use in larger scale modeling studies; we show the performance of the IO model compared to the mechanistic model in terms of accuracy and speed. The models presented here help achieve two objectives. First, the mechanistic model provides a comprehensive platform to describe spine calcium dynamics based on individual contributing factors. Second, the IO model is trained on the main dynamical features of the mechanistic model and enables nonlinear spine calcium modeling on the cell and network level simulation scales. Utilizing both model representations provide a multi-level perspective on calcium dynamics, originating from the molecular interactions at spines and propagating the effects to higher levels of activity involved in network behavior.

Introduction

The calcium ion is a key biochemical signaling molecule for cellular function, and a number of studies have demonstrated its importance in numerous cell types. Calcium is known to be involved in regulation of gene transcription factors (Bading et al., 1993; Dolmetsch et al., 1998), muscle contraction (Ebashi and Endo, 1968; Weber and Murray, 1973), bone formation (Zhu and Prince, 2012), cell metabolism (Contreras et al., 2010) and apoptosis (Mattson and Chan, 2003). Within neurons, calcium has an especially critical role in modulating communication and network activity (Zucker, 1999; Emptage et al., 2001). Calcium has been extensively investigated in postsynaptic spines due to its involvement in various signaling cascades that lead to synapse formation and plasticity, cellular mechanisms which underlie learning and memory.

Early experimental studies on spine calcium focused on measuring calcium concentration changes at postsynaptic spines of the CA1 pyramidal cell in response to presynaptic stimulations and action potential transients through fluorescence procedures (Sabatini et al., 2002). Other studies focused on identifying sources of calcium influx, which include voltage dependent calcium channels (Bloodgood and Sabatini, 2007), intracellular calcium stores (Holbro et al., 2009), and NMDA receptor channels (Bloodgood and Sabatini, 2009). Downstream calcium signaling pathways have also been investigated, where calcium microdomains and localized calcium signaling (Higley and Sabatini, 2012) can invoke signaling of the ubiquitous Calmodulin/CAMKII protein, which leads to AMPA receptor upregulation (Naoki et al., 2005; Zhabotinsky et al., 2006). These standalone studies have helped further our understanding of numerous processes that regulate calcium dynamics in spines, but more research is needed to study how such processes interact with each other and together influence synaptic transmission.

Computational models have also been successfully adapted to the study of calcium dynamics in spines and neurons (Shouval et al., 2002; Naoki et al., 2005; Bartol et al., 2015). The advantage of using computer models over experimental protocols is their inherent ability to provide a controlled environment and overcome limitations in size constraints, an issue common when studying subcellular spaces such as synaptic compartments. Calcium dynamics models vary in biophysiological detail and accuracy, ranging from simple phenomonological models that directly relate calcium concentration to synaptic plasticity (Shouval et al., 2002), to detailed and complete reconstruction of the calcium dynamics at a dendritic subsection of a hippocampal CA1 neuron using stochastic Monte Carlo simulations (Bartol et al., 2015). Typically, the degree of realism used for a model is often dependent on the scale of the study, where calcium models at the cellular or network scales have less physiological detail than molecular scale models.

Scientific computation has been trending toward multi-scale modeling, where models are developed to explore biological phenomena across multiple length or hierarchical scales (Yu and Bagheri, 2016; Seo and Jun, 2017). Neural computation spans molecular (Naoki et al., 2005; Bartol et al., 2015), cellular (Jarsky et al., 2005; Migliore and Migliore, 2012), network (Hendrickson et al., 2015), and cortical systems (Markram, 2006) hierarchical scales. Thus, multi-scale modeling is especially valuable for studying calcium dynamics in neurons, because calcium induced molecular signaling cascades can have dramatic effects on patterns of activity at the cellular/network level. Calcium activity in a CA1 pyramidal neuron is not distributed equally as evidenced in experimental studies (Higley and Sabatini, 2008). Rather, calcium processes can be categorized based on location and degree of influence (Figure 1). For example, the high-voltage activated (HVA) calcium channels present on postsynaptic spines can create localized calcium microdomains—brief, high concentrations of calcium in a small area—which can then lead to protein kinase activation and induce secondary messenger cascades, eventually resulting in synaptic plasticity and larger scale changes in cell properties and network activity (Higley and Sabatini, 2012). Meanwhile, intracellular calcium stores, such as the smooth endoplasmic reticulum (ER) and mitochondria, in dendritic compartments can be involved in more metabolic processes such as gene transcription and ATP production (Li et al., 2004). If possible, use of detailed, biologically accurate computational models in large-scale simulations would add the benefit of monitoring such significant molecular-level influences in large network interaction; however, simulating such a large number of complex mechanisms and interactions in a realistic model requires a prohibitively high computational cost. Unfortunately, this computational burden limits the capability of current models to elucidate critical calcium-dependent mechanisms associated with plasticity, learning, and memory, which emerge from network level activity.

FIGURE 1

Figure 1. Calcium is regulated within a CA1 pyramidal cell neuron on multiple scales—from spines, to dendrites, to the entire cell itself. Presented here is a schematic representation of the regulation of calcium within a CA1 pyramidal cell neuron, as well as some of its functional properties. At the molecular level, calcium dynamics at the postsynaptic spine are controlled by various ion channels and pumps that respond to synaptic activity. Opening of high-voltage activated (HVA) calcium channels can influence other key processes such as localization of calcium microdomains, which can lead to functional changes brought about by calcium-dependent protein kinases such as activation of secondary messenger pathways and synaptic plasticity. Calcium dynamics are also important in dendrites, where many intracellular calcium stores, such as the smooth endoplasmic reticulum (ER) and mitochondria, are located. The presence of calcium stores allows calcium to influence local regulation of factors such as gene transcription and ATP production and regulation. The various sub-cellular calcium dynamics integrated together make up the changes in cell and network activity of the neuron. Our current work in this manuscript is highlighted in the dashed red rectangle in the presented hierarchy: we focus on the development of a model of calcium at the postsynaptic spine, and how our model may be applied in future work to higher hierarchies.

To adequately model calcium-influenced cellular and network level behaviors, it is critical to construct a model that can efficiently and accurately replicate nonlinear calcium dynamics based on the numerous calcium processes on multiple levels, starting at the spine.

The focus of this article is to present a model describing the calcium dynamics at the postsynaptic spine of a CA1 pyramidal cell, as well as a methodology to adapt the previously defined model for multi-scale simulations. The calcium model presented aims to (1) consider the variety of subcellular processes that influence calcium at the spine, and (2) enable multi-scale simulations that include the influence of said subcellular processes on calcium dynamics on a larger scale. In (1), we develop a mechanistic model that consists of various kinetic state models of receptor channels and pumps pertaining to calcium regulation at the spine; the mechanistic model is validated with experimental data from the literature and is used to study the subcellular processes involved in spine calcium dynamics. In (2), we implement a “input-output” model using the Volterra functional series trained on the nonlinear calcium dynamics from the mechanistic model; the Volterra functional series is a nonlinear systems filter that has been adapted previously to successfully model dynamics of nonlinear systems with reduced computational cost (Marmarelis and Marmarelis, 1978; Bharathy et al., 2008; Song et al., 2009a,b; Berger et al., 2010; Tu et al., 2012; Hu et al., 2015). We demonstrate that the IO model requires less time to simulate than the mechanistic model, while still reproducing the complex nonlinear dynamics that arise from calcium interactions at the spine. Thus, we propose that the IO model is better suited for multi-scale modeling of calcium. In future studies, this will allow us to investigate the effects of how the various subcellular mechanisms in which affect spine calcium dynamics can influence cell to cell interactions at the cellular and network levels.

Materials and Methods

Here we first describe the models and parameters used for the mechanistic model; an overview diagram is presented in Figure 2, and a summary of the parameters is provided in Table 1. A precursor of the mechanistic model had been described previously in a conference publication (Hu et al., 2016), which was at the time an incomplete; this current model now describes the model in full with additional mechanisms, optimized parameters, and validation that is described in the results section. The development of the Input-Output model is described afterwards.

FIGURE 2

Figure 2. Postsynaptic calcium can be influenced by a number of factors. Our focus is on calcium dynamics at the postsynaptic spine, as was highlighted in the red dashed rectangle in Figure 1. Models for receptor channels and pumps that contribute to calcium dynamics have been developed in several studies. In our model, these sources are integrated into a single platform to construct a mechanistic model of postsynaptic calcium dynamics.

TABLE 1

Table 1. List of the parameters, specifications, and models used in the mechanistic model and their sources.

Mechanistic Model

The premise of the mechanistic model is to build a physiological representation of the postsynaptic spine taking into account the components that are known to significantly influence spine calcium dynamics. These components are identified based on a number of experiments, reviews, and models found in the literature (Sabatini et al., 2002; Bloodgood and Sabatini, 2007; Higley and Sabatini, 2012; Bartol et al., 2015). In the mechanistic model, influx and efflux components such as the calcium channels (NMDAr, VDCC) and calcium pumps (PMCA, NCX) are represented as calcium current sources which can add or remove calcium in the spine compartment. As calcium ions flow in and out of the spine, the concentration is determined by calculating the change in calcium divided by spine volume—the standard definition of concentration. Buffers and intracellular calcium stores interact directly with the calcium within the spine using reaction rate equations. A schematic diagram of the interactions between the components in the model is presented in Figure 3A.

FIGURE 3

Figure 3. A visual depiction of (A) mechanistic calcium model and (B) the Input-Output spine calcium model. (A) portrays an interaction diagram of the mechanistic calcium model presented in the manuscript, with references the equations and kinetic schemas which govern the dynamics of each of the individual components. (B) In the IO model, the inputs and output correspond to parameters in the mechanistic model. Inputs to the IO model are NMDAr conductance and postsynaptic potential. Output is the estimate on calcium concentration in the spine. Coefficients are estimated based on the spine calcium mechanistic model, and presented in Supplementary Table 6.

Besides modeling calcium in itself, each synapse also provides postsynaptic current to the model neuron. The postsynaptic current influences the postsynaptic potential, which later influences spine calcium via voltage-dependent calcium mechanisms. The postsynaptic current is governed by both AMPA receptor channel and NMDA receptor channel kinetic rate models. The NMDA receptor channel kinetic rate model is the same model described in this Supplementary Figure 1. The AMPA receptor channel rate model is the model described in Robert and Howe (2003). These components are part of the synaptic platform which has been developed in our lab and further details can be seen in Bouteiller et al. (2011), Allam et al. (2015), and Hu et al. (2015).

Spine Volume and Diffusion

The concentration of any constituent depends on the volume of its compartment. Spines come in a variety of shapes and sizes, with different forms such as thin spines, stubby spines, and mushroom spines (Stewart et al., 2005). The postsynaptic calcium model presented here considers the composition of an average mushroom spine of a CA1 neuron, with a spine volume of 0.1 μm³ and the postsynaptic density volume set to 0.0032 μm³ (Stewart et al., 2005). Here, we model the postsynaptic density as separate compartment due to the rapid increase in local calcium concentration during ion channel activation; the calcium concentration then diffuses into the rest of the spine. Diffusion from the postsynaptic density (PSD) to the spine compartment was approximated using the spine diffusion model described by Naoki et al. (2005). The diffusion rate is modeled as:

\begin{array}{l} {(\frac{d [C a^{2 +}]}{d t})}_{d i f f u s i o n} = D \frac{A}{d^{*} V_{s p i n e}} ({[C a^{2 +}]}_{P S D} - [C a^{2 +}]) & (1) \end{array}

\begin{array}{l} \frac{d {[C a^{2 +}]}_{P S D}}{d t} = - \frac{I_{N M D A}}{2 F V_{P S D}} - D \frac{A}{d^{*} V_{P S D}} ({[C a^{2 +}]}_{P S D} - [C a^{2 +}]) & (2) \end{array}

where A is the surface area between the PSD and the rest of the spine volume, d is the distance between the midpoints of the PSD and the spine, V_spine is the total volume of the spine, and ${[C a^{2 +}]}_{P S D}$ and ${[C a^{2 +}]}_{S p i n e}$ is the calcium concentration at the postsynaptic density and the spine, respectively. The postsynaptic density is dependent on I_NMDA, which is the calcium influx current from NMDA receptors. The equation for I_NMDA is described in Equation (3). The influx depends on Faraday's constant F and the volume of the postsynaptic density V_PSD.

In our model, we consider the spine to be an isolated compartment from the rest of the cell, where the calcium does not flow out of the spine into the dendrite or vice versa. We make this assumption based on experimental evidence that spines are isolated electrically (Grunditz et al., 2008), and that less than 0.01% of the total calcium flux into the spine comes from the adjacent dendrite (Sabatini et al., 2002). As such, the model described here assumes that calcium exchange with the dendrite is insignificant and calcium dynamics from our spine model does not influence dendritic potential in our neuron model. We are aware that while dendritic calcium can be influenced by spine calcium dynamics, but in the current study, the focus of our model considers only the spine calcium and not dendrite calcium dynamics, as highlighted specifically as the red area in Figure 1. Future studies are planned that will expand our model beyond the spine level and integrate calcium across more hierarchies that had been presented in Figure 1. In the meanwhile, dendritic potential can influence spine calcium influx dynamics, as described in the following section.

Calcium Influx

Here, we first describe the models for NMDAr and VDCC, then explain protocols to simulate calcium influx mechanisms and activate the aformentioned models.

The calcium current contribution of the NMDA channel is calculated as:

\begin{array}{l} I_{N M D A} = g_{t o t a l N M D A} \times V_{C a} & (3) \end{array}

\begin{array}{l} g_{t o t a l N M D A} = g_{N M D A} \times n_{N M D A} & (4) \end{array}

I_NMDA indicates total calcium current that flows into the spine, g_NMDA is the NMDA conductance to be described in the following paragraph and n_NMDA represents the number of NMDA receptors with the given PSD volume and is set to 20 (Racca et al., 2000). V_Ca represents the calcium potential difference between intracellular and extracellular calcium, which again will be described in more detail later in the text.

The NMDA receptor channel is represented as a kinetic states model as described in Erreger et al. (2005). The model consists of 8 states that represent the resting, activation, opening, and desensitization states of the NMDA receptor channel. The magnesium ion blockade of the NMDAr channel also must be considered. To do so, the magnesium ion binding properties in the channel pore are described by:

\begin{array}{l} g_{0} = g_{1} + \frac{g_{2} - g_{1}}{1 + e^{α V}} & (5) \end{array}

\begin{array}{l} g_{max} = \frac{g_{0}}{1 + (\frac{{[M g^{2 +}]}_{o}}{K_{0}}) e^{- δ z F V / R T}} & (6) \end{array}

\begin{array}{l} g_{N M D A} = g_{max} \times O (t) & (7) \end{array}

Where g₀ represents the total conductance in the absence of any magnesium, g₁ and g₂ represent the open state conductances with one glutamate bound and two glutamate molecules bound, respectively. g₁is set at 40 pS while g₂ is set at 247 pS. The value α = 0.01 represents the steepness of the transition between g₁and g₂. ${[M g^{2 +}]}_{o}$ represents the external magnesium concentration and is set at 1 mM. K₀ is the equilibrium constant for magnesium set at 3.57, F is Faraday's Constant (9.64867.104 C mol⁻¹), R is the molecular gas constant (8.31434 J mol⁻¹ K⁻¹), z is the valency of the calcium ion (2), and T is the temperature at 299.5 K. V represents the membrane potential. δ is the affinity between NMDAr and magnesium, which is dependent on the postsynaptic potential of the synapse; the value is set to 0.8. The variable O(t) is the open state probability governed by the kinetic rate equations for the NMDA model. The kinetic schema of the NMDAr model is presented in Supplementary Figure 1 where O(t) is represented as “Open.” The rate constants which govern NMDA kinetics are presented in Supplementary Table 1. In the kinetic model, presynaptic action potentials cause vesicle glutamate release. Thus, in our model, a presynaptic event correlates to a glutamate spike in the NMDAr kinetic states model, moving the NMDAr channel kinetics away from resting state and causing them to open. For a more extensive description of the kinetic NMDAr model, please refer to Erreger et al. (2005).

NMDAr is differentially permeable to different ions when activated. Burnashev et al. (1995) reported that on average, calcium constitutes about 11% of the total ion current when NMDAr channels are opened. However, considering NMDA current alone lead to an incorrect representation of NMDAr mediated calcium influx, since the reversal potential of calcium (+ 50 mV) is considerably different from the reversal potential of NMDAr conductance (+0 mV). Therefore, in this model we calculate the influx using the Nernst equation, which instead depends on the difference between intracellular and extracellular calcium concentration. The Nernst equation used is as follows:

\begin{array}{l} V_{C a} = - \frac{R T}{2 F} \log (\frac{{[C a^{2 +}]}_{o}}{[C a^{2 +}]}) & (8) \end{array}

V_Ca represents the potential difference for calcium. R, T, and F represent the molecular gas constant, temperature, and Faraday's constant with the same values mentioned previously. ${[C a^{2 +}]}_{o}$ is the extracellular calcium set as a constant concentration of 2 mM, while [Ca²⁺] represents intracellular calcium concentration at the spine. Here, we assume the calcium at the spine dominates the driving force represented by the Nernst potential. It is also possible to use the calcium concentration at the postsynaptic density to calculate the potential, but due to its small volume, minor fluctuations in current can result in drastic changes in the concentration value, which would lead to more erratic changes in the potential. Therefore, we consider the use of the overall spine concentration to be a more adequate representation of the Nernst potential for calcium.

The Nernst potential is typically used to determine the potential when the spine is at during the resting state. However, for our purposes we use the Nernst potential as an estimate of the driving force for calcium influx. Our justification is thus: In our spine model, we assume that all ions besides calcium remain constant and the electrochemical force for all other ions is zero. As mentioned earlier, considering only a percentage of synaptic current influx as calcium is inaccurate, since the reversal potential between calcium is much higher than the reversal potential of synaptic current. This case is true even for the Goldman–Hodgkin–Katz flux equation, since it depends on the membrane potential—if the membrane potential moves from negative to positive, the flux is also reversed. On the other hand, using a constant value for calcium potential, i.e., when calcium is at rest at 50 nM, does not account for changes in flux induced by increased calcium levels in the spine. Therefore, the use of the Nernst potential is an estimation of the electrochemical gradient when calcium concentration is changed. The concentration of intracellular calcium is orders of magnitude lower than the concentration of extracellular calcium, so the change is not large, but still significant enough such that we believe it should be accounted for, i.e., the value of V_Ca when intracellular calcium is 50 nM is approximately −60 mV, but can reach −30 mV when intracellular calcium reaches at 10 μM.

Voltage dependent calcium channels (VDCCs) let calcium into the spine when there is a significant difference in membrane potential, such as action potential backpropagation from the postsynaptic neuron (Higley and Sabatini, 2012). Various types of VDCCs exist—each having different channel properties, mechanics, and functional roles—and are found on different types of cells (Catterall, 2011). For CA1 pyramidal cells, experimental evidence suggests that T-type VDCCs contribute the most to overall calcium concentration within dendritic spines (Bloodgood and Sabatini, 2007). It should be noted that L type and R type VDCCs are also present. However, the calcium influx contribution of L-types and R-types to the overall calcium concentration within spines was found to be insignificant. Instead, these channels tend to be concentrated into microdomains and activate secondary messenger pathways (Higley and Sabatini, 2012). Consequently, we consider the VDCC influx through T-type channels only.

The more specific details of the VDCC model are described in Supplementary Table 2. In general, the calcium contribution of the VDCC channel is calculated as:

\begin{array}{l} I_{V D C C} = g_{V D C C} \cdot m^{2} \cdot h \cdot f_{d r i v e} & (9) \end{array}

Where I_VDCC is the calcium current from the voltage dependent calcium channel. m and h are part of the Hodgkin-Huxley equation, with parameters as defined from Jaffe et al. (1994). f_drive is the driving force of the internal and external calcium dynamics, considered through modifications to the Hodgkin Huxley equation as reported by Poirazi et al. (2003) (equivalent to dvf in their model). Once again, the equations and parameters used in the model related to VDCC are presented in Supplementary Table 2. g_VDCC is the he single channel conductance for VDCC is set to be 7.5 pS (Perez-Reyes et al., 1998) and the average number of VDCCs opened for each AP-evoked transient is 5 (Sabatini and Svoboda, 2000).

Calcium Efflux

Experimental evidence indicates that calcium is removed from the intracellular space through pumps and exchangers, but mechanistic details concerning calcium efflux at the postsynaptic spine are not yet fully understood. Generally, Plasma Membrane Calcium Pumps (PMCA) and Sodium-Calcium Exchangers (NCX) are the two prominent elements that participate in calcium extrusion in spines and small dendrites of CA1 pyramidal cell neurons (Scheuss et al., 2006). Overall it is thought that the constant active pumping by PMCA helps maintain the standard basal levels of calcium at ~50 nM (Carafoli, 1991), while NCX helps to quickly extrude calcium in a short amount of time, such as during an action potential (Carafoli et al., 2001). One source in the literature suggests the PMCA isoform is type PMCA2w, although details on the dynamics and extrusion rates of the isoform are lacking (Burette et al., 2010). Details on NCX at spines are even less studied, with only one source indicating that NCX is present in larger numbers in dendritic shafts than in spines, though exact numbers are unknown (Lörincz et al., 2007).

Previously published models of spine calcium use calcium extrusion models as a calibration factor to fit simulations to experimental results (Naoki et al., 2005; Bartol et al., 2015). For our model, we are interested in using more physiologically accurate representations of extrusion, but such models specifically relating to the extrusion pumps and channels in spines are currently absent from the literature. Therefore, we have decided to use extrusion models from models of other physiological systems (Gall and Susa, 1999; Fridlyand et al., 2003), and adjust parameters according to the geometry of the synapse based on the surface area density of the models presented within these papers. The adjusted parameters are detailed in the following paragraph.

Both PMCA and NCX are represented through the Hill equation, which is typically used to describe binding properties of a ligand to the receptor:

\begin{array}{l} I = I_{M A X} \frac{{[C a^{2 +}]}^{h_{c}}}{K^{h_{c}} + {[C a^{2 +}]}^{h_{c}}} & (10) \end{array}

Here, I_MAX is the maximum calcium current extruded from the model, K is the equilibrium constant and h_c is the hill coefficient. For PMCA model parameters, the Hill coefficient is set to 2 and an equilibrium constant of 1 μM (Fridlyand et al., 2003). The maximum calcium current was optimized through gradient descent based on experimental protocols highlighted in the results section; the optimized maximal PMCA current was set to be 0.25 pA. The NCX model has a Hill coefficient of 5, equilibrium constant value of 1.5 μM and a conductance of 0.0117 pS (Gall and Susa, 1999); maximal current is then calculated from the sodium/calcium gradient and the given conductance. The higher Hill coefficient represented in the NCX model indicates a higher affinity for calcium ions when more calcium is bound. A table describing the parameter values is presented in Table 2.

TABLE 2

Table 2. List of calcium efflux parameters.

Additionally, the NCX model accounts for the sodium/calcium gradient in which three sodium ions are exchanged for one calcium ion (Fridlyand et al., 2003):

\begin{array}{l} V_{N a, C a} = \frac{R T}{F} (3 ln (\frac{{[N a^{+}]}_{o}}{{[N a^{+}]}_{i}}) - ln (\frac{{[C a^{2 +}]}_{o}}{{[C a^{2 +}]}_{i}})) & (11) \end{array}

R, T, and F represent the molecular gas constant, temperature, and Faraday's constant with the same values mentioned previously. I_MAX for NCX is determined based on the gradient difference. In the current model, intracellular and extracellular sodium concentrations are considered constant at 10 and 140 mM respectively. Similar to V_Ca from the Nernst equation, V_{Na, Ca} is an estimate of the driving force to consider how calcium efflux will change depending on calcium concentration, providing a slightly better estimate compared to keeping V_{Na, Ca} constant.

Buffers and Intracellular Calcium Stores

Buffering is an integral part of calcium dynamics at the postsynaptic spine, as up to 95% of the total intracellular calcium is bound to buffers (Sabatini et al., 2002). There are numerous types of buffers that can bind calcium; in our model, we specify two types of buffers, calmodulin and calbindin, while the other possible buffers are represented as a collection of generic calcium binding proteins (CBP). In our platform, the buffer models directly influence the free calcium concentration in the spine using reaction rate equations. The kinetic schemas and descriptions of the parameters and equations are presented in Supplementary Figure 2, Table 1, and Supplementary Table 3. Calmodulin is a ubiquitous calcium buffer which plays a role in AMPA receptor upregulation and synaptic potentiation when found in the postsynaptic spine. The calmodulin buffering parameters are defined in accordance to Zhabotinsky et al. (2006). Calbindin is a binding protein with four calcium binding sites; here it is defined as a 9 states kinetic model, with parameters defined in the calcium model by Bartol et al. (2015). The CBP were calibrated after the previous two buffers were implemented, where the total buffered calcium at steady state reaches approximately 95% in the presence of all three buffers.

Intracellular calcium stores are known to play a large role in calcium dynamics, but current evidence on its impact particularly on dendritic spines in CA1 neurons remains controversial. So far it is found that approximately 19% of the total spine count contain endoplasmic reticulum (ER), with a majority of the ER-containing spines having a larger volume than others (approx. 0.06 μm³) (Holbro et al., 2009). The ER apparatus within the spine was shown to have no IP3 receptors present while retaining a number of ryanodine receptors (Paula-Lima et al., 2014). As such, we have included in our model state representations for SERCA pumps and ryanodine receptors, but omit IP3 receptors. Just like the buffer models, the models pertaining to the intracellular calcium stores directly interact with the free calcium concentration in the spine. We describe the kinetic schema, parameters, and equations of SERCA and ryanodine receptors in Supplementary Figure 3, Supplementary Tables 4, 5. The SERCA pump is a 2 states model with equations and parameters derived from Higgins et al. (2006). We also included the ryanodine receptors model proposed by Williams et al. (2011).

Inputs Into the Mechanistic Spine Calcium Model

In Figure 3A, we provide a diagram of the interactions in the mechanistic model and describe the components which can influence calcium activity. Input stimulation predominantly occurs in two ways: (1) through synaptically activated transients, where presynaptic release of neurotransmitter activates the NMDA receptor channels on the postsynaptic density; and (2) AP-evoked transients, where stimulation of the postsynaptic neuron triggers an action potential, which is then backpropagated to the spines. Simulation of (1) is represented through presynaptic event-based activation of the NMDA receptor model, where a single event triggers the opening kinetics of the NMDAr model. The protocol for (2) is slightly more complex: in order to simulate AP-evoked transients in our model, we stimulate a number of synapses on the neuron model to invoke an action potential. Calcium concentration can then be measured on a stimulated or non-stimulated synapse, where the back-propagating action potential opens VDCCs and prompts calcium entry into the spine.

Input-Output Model

The development of the Input-Output (IO) model for postsynaptic calcium dynamics follows a protocol similar to the IO models that had been covered in Berger et al. (2010) and Hu et al. (2015). To describe briefly, the IO model uses the Volterra functional power series, with Laguerre functions as the basis equations of the Volterra series. The single input, single output (SISO) model of the Volterra functional series and the Laguerre equations have been previously described in Hu et al. (2015). In brief, we propose to use a simplified functional representation of the system under consideration.

In the current work, a multi-input, single output (MISO) model was developed with the inputs being the two major sources of calcium influx for the calcium model: (1) postsynaptic potential (V), and (2) NMDAr channel conductance (g_totalNMDA) based on glutamate-based calcium influx. One notable difference is that we herein consider the inputs to be continuous. Although we use spike trains as input to our synapse platform, the conductance and potential which are calculated from the spike trains are continuous. The use of continuous inputs contrasts with IO models that have been published previously assumed the inputs to be of binary nature (Hu et al., 2015). A pictorial representation of the mechanistic model and IO model is shown in Figure 3B, to highlight which parameters the two models have in common.

The inputs of the IO model were chosen because they represent outside activity that drives changes in calcium dynamics. The first input, the postsynaptic potential, is a factor known to drive calcium influx due to its effects on voltage dependent mechanisms (NMDAr associated channel and calcium channels). In both mechanistic and IO models, the postsynaptic potential depends on the neuron cell model. In our simulations, the Izhikevich model (Izhikevich, 2003) is used for most simulations except simulations that involve distance measurements, where the Migliore neuron model is used instead (Migliore and Migliore, 2012). Since the Izhikevich model has no geometry component, measurements with respect to distance cannot be performed using the Izhikevich model. The Migliore model, on the other hand, is a reconstruction of a hippocampal CA1 neuron, which consequently has dendritic geometry that can be used for our simulation study. In the neuron model, we simulate synaptic activity which results in firing of the cell; the action potential is then backpropagated to the spine model, which we use as input to either the mechanistic or the IO model. The second input, i.e., NMDA-R channel conductance, is a critical measure of pre- and postsynaptic activity; we determine the total NMDAr conductance based on both pre- and post-synaptic activity of the spine.

In our platform, opening of the NMDA receptor channel model: (1) Allows calcium influx into the spine, and (2) produces postsynaptic currents, which are passed to the neuron model and influences the postsynaptic potential. We chose the NMDAr channel conductance as an input parameter to the input-output model to account for the calcium influx while still allowing NMDAr channels to influence postsynaptic activity. The output response of the IO model is calcium concentration, calibrated using the calcium concentration obtained with the mechanistic model. Thus, the influence on calcium dynamics of all other components besides the NMDAr model are captured in the IO model. The use of the IO model then allows us to model complex nonlinear calcium dynamics without requiring the large number of components that would otherwise be necessary when using a mechanistic model.

To describe the structure of the IO model, we first begin with a description of the SISO Volterra series:

\begin{array}{l} u_{S I S O} (t) = c_{0} + \sum_{j = 1}^{L} c_{1} (j) v_{j} (t) + \sum_{j_{1} = 1}^{L} \sum_{j_{2}}^{j_{1}} c_{2 s} (j_{1}, j_{2}) v_{j_{1}} (t) v_{j 2} (t) \\ + \sum_{j_{1} = 1}^{L} \sum_{j_{2} = 1}^{j_{1}} \sum_{j_{3} = 1}^{j_{2}} c_{3 s} (j_{1}, j_{2}, j_{3}) v_{j_{1}} (t) v_{j 2} (t) v_{j 3} (t) & (12) \end{array}

\begin{array}{l} v_{j} (t) = \sum_{τ = 0}^{M} b_{j} (τ) x (t - τ) & (13) \end{array}

Where u_SISO is the single input Volterra series up to 3rd order. L refers to the total Laguerre functions in the SISO model, and c₀, c₁, c_2s, and c_3s refer to the coefficients associated with the 0th, 1st, 2nd, and 3rd order response of the series, respectively. v_j(t) is the convolution of the input to the IO model with the basis function; x(t−τ) is the input to the IO model; and b_j(τ) is the j-th basis function. M is the memory window of the IO model, set to 5 s. The coefficients are determined during the training process to best fit the nonlinear response of calcium dynamics. For our case, we use the Laguerre basis functions for our model. The Laguerre equations are used for their orthogonality and convergence properties, which are characteristic of many biophysiological systems (Berger and Song, 2010; Ghaderi et al., 2011).

The structure of the Volterra functional series in the MISO model is similar to the SISO model. In the case of two inputs, the series consists of the summation of two SISO model components that account for two different inputs, then adding a cross-kernel component accounting for possible nonlinear interactions that may occur due to the presence of two inputs. The equations become thus:

\begin{array}{l} u_{M I S O} (t) = c_{0} + u_{1} (t) + u_{2} (t) \\ + \sum_{k_{1} = 1}^{L_{1}} \sum_{k_{2} = 1}^{L_{2}} c_{2 r} (k_{1}, k_{2}) v_{k 1}^{u_{1}} (t) v_{k 2}^{u_{2}} (t) \\ + \sum_{k_{1} = 1}^{L_{1}} \sum_{k_{2} = 1}^{L_{2}} \sum_{k_{3} = 1}^{L_{1}} c_{3 r_{1}} (k_{1}, k_{2}, k_{3}) v_{k 1}^{u_{1}} (t) v_{k 2}^{u_{2}} (t) v_{k 3}^{u_{1}} (t) \\ + \sum_{k_{1} = 1}^{L_{1}} \sum_{k_{2} = 1}^{L_{2}} \sum_{k_{4} = 1}^{L_{2}} c_{3 r_{2}} (k_{1}, k_{2}, k_{4}) v_{k 1}^{u_{1}} (t) v_{k 2}^{u_{2}} (t) v_{k 4}^{u_{2}} (t) & (14) \end{array}

\begin{array}{l} v_{j}^{u_{i}} (t) = \sum_{τ = 0}^{M} b_{j} (τ) x^{u_{i}} (t - τ) & (15) \end{array}

Here, u_MISO is the equation for the multi-input Volterra series, and u₁ and u₂ are the single input Volterra series based on u_SISO with the inputs the postsynaptic potential (V), and the NMDA receptor conductance (g_totalNMDA), respectively. The cross kernel components involve the number of Laguerre functions L₁ and L₂ from the first and second SISO components, respectively, and at all orders. $v_{k}^{u} (t)$ is the associated basis functions convolved with inputs: postsynaptic potential (V), and the NMDA receptor conductance (g_totalNMDA). c_r is then the associated coefficients for these terms. As a result, the nonlinearities influenced by having two different inputs are considered.

The IO model must be trained to tune its parameters and minimize the error with respect to the original system—in this case, the mechanistic model. The training process consisted in using the response of the original mechanistic model to random Poisson train events given at both low average frequency (2 Hz) and high average frequency (10 Hz), with a total of 1,000 events to ensure a wide range of activity. The low average frequency of 2 Hz was chosen to reflect the typical firing rate of hippocampal CA1 and CA3 neurons during resting state (Berger et al., 1988). Meanwhile, the high average frequency of 10 Hz is used to account for situations with higher levels of neuron activity in physiological conditions (Ranck, 1973). The coefficients associated with the basis functions were determined through the pseudoinverse matrix multiplication, and the optimal decay values associated with the Volterra series were estimated using gradient descent to find the parameters that resulted in the lowest root mean square error value in the comparison between the estimated and the actual result from the original model.

Platforms and Computational Tools

The mechanistic calcium model is described in the Systems Biology Markup Language (SBML) and was constructed using CellDesigner, a visual diagram editor for SBML (Funahashi et al., 2008). Our models were then run in MEMORY, a python-based platform designed for synapse-based simulations using the libSBML library and the simulation engine libroadrunner (Somogyi et al., 2015). In our studies which involve measuring calcium and dendritic potential based on distance from the soma, we required a neuron model with a detailed dendritic structure where synapses can be placed along the dendritic arbor. To achieve this, we use the CA1 pyramidal cell model designed by Migliore (Migliore and Migliore, 2012) simulated within NEURON cell simulation platform (Carnevale and Hines, 2006). The IO calcium model first required the coefficients to be determined, which were calculated within MATLAB; afterwards, the model was implemented into the MEMORY platform to run for simulation. Unless otherwise indicated, the input stimulus provided to synapses during simulation were Poisson 2 Hz randomized input trains, where each synapse was provided a unique randomized input. All simulations were conducted on a computer using the Fedora OS, with Intel quad-core 2.67 GHz processor and 8 Gb RAM.

Results

Calcium Dynamics Calibration and Validation With Published Experimental Data

In the resting phase, the average cytosolic calcium concentration in spines is typically kept at around ~50 nM, maintained by various pumps and buffers (Higley and Sabatini, 2012). Calcium influx occurs during activation of the various channels present on the spine. There are two major sources of calcium influx: Glutamate-dependent calcium influx, where calcium flows into the spine via NMDA receptor channels, and Voltage-dependent calcium influx, primarily through voltage-dependent calcium channels. For calibration of the mechanistic calcium model at the spine, we consider two scenarios: (1) calcium influx due to presynaptic activation, where the presynaptic terminal releases glutamate in response to a presynaptic action potential, and (2) calcium levels when the postsynaptic neuron is fired, leading to a backpropagating action potential (bAP). The response in (1) has been measured both in the presence and absence of postsynaptic depolarization (Sabatini et al., 2002; Higley and Sabatini, 2012). Our calcium model has a response of approximately 9 and 0.9 μM in the presence and absence of postsynaptic depolarization, respectively, which is in line with the experimental data presented in the literature (Figure 4A). Similarly, calcium levels in response to bAP were simulated and compared to calcium levels measured in Sabatini et al. (2002) (Figure 4B). The backpropagation factor assumed that the spine was a distance of ~150 μm from the soma, similar to the synapses recorded in the literature. The simulated results showed an amplitude of approximately 700 nM, comparable to measurements from Sabatini et al. (2002).

FIGURE 4

Figure 4. Validation of single event responses of the mechanistic model of calcium concentration at the postsynaptic spine. (A) The response of the postsynaptic spine model to a single presynaptic event with an elicited action potential (in light blue, dashed) that reaches ~9 μM, and without any resulting action potential (in blue, solid) that reaches ~1 μM, in line with expected amplitudes that have been reported in experimental studies and reviews (Sabatini et al., 2002; Higley and Sabatini, 2012). Neurotransmitter release following a presynaptic event opens NMDA receptor channels, causing calcium influx; the degree of influx is also dependent on membrane potential due to the voltage-dependent magnesium block. (B) Response to a single postsynaptic backpropagation event by the postsynaptic spine model (in red) in comparison to reported calcium response seen in the study by Sabatini et al. (2002) (in black). Voltage-dependent calcium channels (VDCC) open in response to a change in calcium levels—the result shown here is the response due to low-voltage activated VDCC-T type channels in the postsynaptic spine model.

Calcium Fluctuations Vary According to the Inter-spike Intervals Between Presynaptic and Postsynaptic Activity of the Neuron

Due to the number of mechanisms in place, the integration of glutamate-dependent calcium influx and voltage-dependent calcium influx leads to further complex dynamics in calcium concentration at the spine. When the membrane potential increases as a result of postsynaptic events, the magnesium block is removed in the NMDA receptor channel pore, increasing the influx of ions when the receptor is activated (Jahr and Stevens, 1990; Ambert et al., 2010). In our model, we measured the maximum amplitude of calcium concentration at varying pre-post intervals to determine its influence on calcium dynamics in the spine (Figure 5). Through the study, we found that the highest calcium concentration peak reached was 9.88 μM, when the postsynaptic event occurred 3 ms after the presynaptic event. It was also noted that for pre-post intervals, where the presynaptic event precedes the postsynaptic event, there is a notable increase in the magnitude in comparison to the calcium response when there is only presynaptic activation with no postsynaptic activation. The increase in maximum calcium amplitude is present when the pre-post interval ranges from 0 to 80 ms; beyond 80 ms, there appears to be no significant change in maximum amplitude compared to presynaptic activation only. In the case of post-pre intervals, where the presynaptic event follows the postsynaptic event (corresponding to negative delays in Figure 5), calcium amplitude begins rising sharply starting at 20 ms all the way to 0 ms. It is interesting to note that while our model does not account for synaptic plasticity, the pre-post interval time-scale dependency of calcium amplitude on pre-post intervals resembles the well-known spike timing dependent plasticity (STDP) curve presented by Bi and Poo (2001), particularly in regards to synaptic strengthening when post-synaptic activation follows pre-synaptic activation. The correlation in timescale dependence between the STDP curve and the presented results from the mechanistic calcium model suggests the mechanistic model can be considered a plausible model basis for a plasticity model in the future.

FIGURE 5

Figure 5. The magnitude of the calcium response changes depending on the interval between presynaptic and postsynaptic events. Timing interval between presynaptic and postsynaptic events can influence the amplitude of calcium dynamics in the mechanistic model of calcium at the postsynaptic spine. Pre-Post interval is defined as the amount of time in milliseconds the postsynaptic event (presynaptic neuron activation, NT release is triggered) occurs after the presynaptic event (postsynaptic neuron activation, the action potential is backpropagated to the spine). Black dots indicate measured maximum calcium response for simulations with the designated pre-post interval. In the inset, A., B., and C. show the simulated calcium profile in response to different pre-post intervals. The change in maximum calcium amplitude reflects the kinetics of the NMDAr channel dynamics, where the voltage dependent magnesium block results in the nonlinear behavior of the calcium response.

Postsynaptic Calcium Activity Depends on the Distance Between the Synapse and the Soma

The experimental studies conducted by Sabatini et al. (2002) observed synapses located approximately 150 μm from the soma. Studies have shown that backpropagation signal properties of the synapse and dendrite can change depending on distance (Golding et al., 2001). To observe the influence of distance on the calcium response, we used the CA1 pyramidal neuron model by Migliore and Migliore (2012) to simulate 50 randomly placed synapses on the stratum radiatum, with the closest synapse having a distance of 72.62 μm and the furthest synapse has a distance of 407.03 μm; these values are close to the minimum and maximum range specified for the stratum radiatum (Megías et al., 2001). When synapses were randomly placed on the pyramidal cell model, we measured the dendritic diameter of the synapse locations. It was found that 45 (90%) of the synapse locations had a diameter of 0.5 μm. This corresponds to the measured diameters of the thin dendrites located in the stratum radiatum in experimental findings (Megías et al., 2001). Other diameters were 2, 2, 1.2, 1.2, and 0.18 μm. Our simulations from these synapses showed no influence of the diameter on our measured results. In our first set of simulations, the synapses were stimulated with a single presynaptic event, then a fixed single postsynaptic event following 10 ms afterwards. Figure 6A shows a schematic of the simulation setup. In Figure 6B, we summarize the results where for each synapse we consider (Figure 6Bi) the resting potential, Figure 6Bii maximum amplitude of the backpropagating action potential during stimulation with a pre-post interval of 10 ms, and Figure 6Biii maximum calcium response during stimulation with a pre-post interval of 10 ms. We note that both the resting potential from Figure 6Bi and the max bAP from Figure 6Bii are properties inherent to the neuron model described by Migliore and Migliore (2012). For Figure 6Biii, we simulate calcium dynamics at the spine using our mechanistic calcium model described in this manuscript. In Figure 6C, we repeated the simulations with different pre-post intervals and measured the maximum calcium response at the spine of each synapse location.

FIGURE 6

Figure 6. Simulated results of dendritic potential and maximum calcium amplitudes as a function of distance from the soma. Fifty random locations were chosen within the stratum radiatum sections of a pyramidal CA1 cell model by Migliore and Migliore (2012). (A) Diagram of the simulation protocol for stimulating and measuring calcium levels for synapses located at different distances from the dendrite. Presynaptic stimulation indicates simulated presynaptic release on the respective synapse, while postsynaptic stimulation indicates simulation of the soma, inducing a backpropagating action potential. (B) Synapses were stimulated with a presynaptic release event at time 0, followed by a simulated current injection into the soma of the neuron 10 ms following the presynaptic event, triggering an action potential and invoking backpropagation. (i) The dendritic resting potential values from the 50 locations with respect to distance. On average, the potential is slightly higher in distal than in proximal locations by approximately 2 mV. (ii) In contrast to the baseline potential, the maximum potential reached with the backpropagation of a single postsynaptic spike becomes more attenuated when further away from the soma. (iii) The influence of the attenuation is seen in the calcium measurements, where maximum calcium is reduced when moving further away from the soma, with 300 μm and beyond the backpropagation becomes negligible. (C) The pre-post interval was varied between simulations, and the maximum calcium response was measured for each simulation at each synapse location. The darker lines indicate calcium responses from synapse locations closer to the soma, while lighter lines are responses from locations further away from the soma. The blue dotted line is where the pre-post interval is 10 milliseconds, reflecting results that are seen in (Bi) through (Biii).

Our simulation results indicated that baseline potential is slightly increased in distal dendrites than in proximal dendrites, with a range lying between −73.4 mV and −71 mV. Such a difference appears rather small, but still constitutes a notable trend with respect to distance. Differences up to 2 mV between somatic and dendritic resting potential have been observed experimentally, and our model falls in line within these constraints (Golding et al., 2001). Conversely, the maximum potential reached after pre-post event stimulations decreases with respect to distance; proximal dendrites are more likely to reach a higher maximum potential than distal synapses.

Observations of the calcium concentration levels at the postsynaptic spine indicated that at a pre-post interval of 10 ms, the amplitude of the calcium concentration peak in spines decreased as distance with the soma increased. Again, beyond 300 μm we found that calcium amplitudes do not extend much beyond 1 μM, likely due to the reduced influence of the postsynaptic activity. Subsequent simulations where the pre-post intervals are changed demonstrated that the influence of the timing between presynaptic and postsynaptic events is more prominent in proximal synapses than in distal synapses. In particular, there no longer appears to be any dependency of the max calcium response on the pre-post interval timing for the synapses farthest away from the soma. These simulations are in line with studies on the influence of distance on STDP, where it was found that backpropagation induces LTP more commonly in proximal synapses, while at distal synapses LTD occurs more frequently in response to the same backpropagating potential (Sjöström and Häusser, 2006).

Results presented here suggest that the bAP is significantly attenuated in distal spines to the degree that a single pre-post events does not trigger much calcium influx. At the time of writing, experimental data of spine calcium levels based on distance have not been documented in the literature. However, the results presented on backpropagation are in line with what has been observed in experimental studies, such as Golding et al. (2001), where it was observed the bAP amplitude is reduced when further away from the soma, and especially beyond 300 μm. The reason for attenuation is likely due to two factors: (1) changes in active conductance with respect to distance from the soma, where at distal dendrites there is a higher density of potassium channels and low density of calcium and sodium channels (Bikbaev et al., 2016); and (2) the branching of the dendritic arbors, which has also been seen to contribute to the attenuation of the bAP (Golding et al., 2001). In distal synapses, there may be other mechanisms at play that may more prominently influence signaling and plasticity to compensate for attenuated bAP, such as modulation by glial cells and neurotransmitters (acetylcholine, brain derived neurotrophic factor, dopamine noradrenaline) (Edelmann et al., 2017).

A Third Order, Multi-Input Input-Output Calcium Model Closely Replicates the Response of the Mechanistic Calcium Model at Lower Frequency

Each individual component in the mechanistic calcium spine model has its own degree of computational complexity, and the integration of all the components also compounds the overall computational burden. As a result, simulation of the calcium dynamics for a larger number of spines becomes increasingly difficult. We demonstrate here the use of an input-output model based on the Volterra functional series that reduces computational cost of simulating calcium dynamics. The output of the model is calcium concentration. The inputs to the proposed model are membrane potential and NMDA receptor conductance. The NMDA receptor model is the only component in the calcium dynamics model that we consider outside of the IO calcium model, since an IO NMDA model had been developed before and can be utilized in its place (Hu et al., 2015). In brief, the IO NMDA model is a single-input-single output model which uses the Volterra series with Laguerre basis functions to predict the open state probability of the NMDA receptor channel (“Open” from Supplementary Figure 1). The open state probability is then used to calculate the conductance based on the magnesium blockade equation (Equations 3–5). However, for the purposes of consistency, in our simulations we use the kinetic NMDA receptor model to properly compare results only between the IO calcium model and mechanistic calcium model. Training the IO model requires keeping track and replicating calcium concentration profiles from the mechanistic model in three separate conditions. First, when only presynaptic stimulation is applied; then, when only postsynaptic stimulation is applied (back-propagated action potential); and finally, when both presynaptic and postsynaptic stimulations are applied to the spine. For each type of stimulation, we used Poisson random interval trains for 1,000 events at 2 Hz and 1,000 events at 10 Hz for a total of 2,000 events, as each input. This gives us a total of 8,000 events for the model to be trained on: 4,000 events in total for presynaptic stimulation, and 4,000 events total for postsynaptic stimulation. The large number of events at lower and higher frequencies gives the IO model an adequate range of nonlinear dynamics to be trained on. At the end of the training phase, the root mean square (RMS) difference is calculated and normalized to the maximum value and minimum of mechanistic model response. The difference between the mechanistic model and the trained IO model was 6.79% with the given training data.

We then validated the trained IO model with a naive train of presynaptic and postsynaptic stimulations for both 2 and 10 Hz and compared the results with the ones obtained with mechanistic model. The trained IO model was found to be more accurate at lower frequencies: the validation error at 2 Hz was 8.15%, while the error for 10 Hz reached 16.9% (Figure 7). We have also provided a comparison between the mechanistic model response and the response from the linear calcium model by Shouval et al. (2002), which is presented in Figure 8. The linear calcium model is presented as:

\begin{array}{l} \frac{d [C a^{2 +}]}{d t} = I_{N M D A} (t) - (\frac{1}{τ_{C a}} [C a^{2 +}]) & (16) \end{array}

FIGURE 7

Figure 7. Comparisons between the calcium responses from the mechanistic model and the IO calcium model. (A,B) Shows the responses of the mechanistic model (in blue) and the trained IO calcium model (in red) over the course of 20 s given Poisson randomized presynaptic and postsynaptic events with an average of 2 and 10 Hz, respectively. The difference between the mechanistic and the IO model are plotted beneath each response. The calculated RMS difference between the two models is 8.15% for the 2 Hz response and 16.9% for the 10 Hz response. (C,D) Shows a direct comparison between the calcium response values from the mechanistic model (x axis) and the IO calcium model (y axis) from the 2 and 10 Hz responses, respectively.

FIGURE 8

Figure 8. Comparison between the linear calcium model (Shouval et al., 2002) and the mechanistic calcium model. (A,B) Shows the responses of the mechanistic model (in blue) and the linear calcium model (in orange) over the course of 20 s given Poisson randomized presynaptic and postsynaptic events with an average of 2 and 10 Hz, respectively. The root mean square difference between the two models is 80.52% for the 2 Hz average response, and 89.75% for the 10 Hz average response. (C,D) Shows a direct comparison between the calcium response values from the mechanistic model (x axis) and the linear calcium model (y axis) from the 2 and 10 Hz responses, respectably.

[Ca²⁺] represents calcium concentration. I_NMDA(t) is the contribution of calcium current provided by the NMDA receptors; in our simulations, we used the NMDAr kinetic rate model described by Erreger et al. (2005) to determine calcium current from NMDA receptors. This is the same NMDAr model used for our mechanistic calcium model and IO calcium model. τ_Ca is the time constant for calcium decay in the linear model. As a result, all nonlinearities associated with NMDAr kinetics are also being accounted for in the simulation with the linear calcium model. We calibrated the parameter to 20 ms, which best approximates the decay of the first order response to the mechanistic model. After calibration, we simulated Poisson random input events with an average frequency of 2 and 10 Hz to the linear model to compare with the mechanistic calcium model response. The difference is shown in Figure 8. The root mean square difference between the two models is very large: 80.52% for the 2 Hz average response, and 89.75% for the 10 Hz average response. This demonstrates that, even with the nonlinear dynamics of the NMDAr model accounted for, considerable nonlinearities in the mechanistic calcium model exist that cannot be replicated in a linear calcium model. In particular, the results from the linear model significantly undershoots the calcium dynamics seen in the mechanistic model as demonstrated in Figures 8C,D. Because the nonlinear dynamics of the NMDA receptor channel have been accounted for, these differences in nonlinearity are more likely a result of the buffers, VDCC, and calcium influx dynamics simulated within the mechanistic model but not in the linear model. Meanwhile, the nonlinear dynamics are reproduced in the IO model, where the RMS error was shown to be much lower.

Following validation, the computational time to run the IO model and the mechanistic model was determined based on number of spines with 2 Hz Poisson random interval train inputs. We find that the IO model finished the simulations faster than the mechanistic model, where the runtime of the IO model required around half the time to finish a simulation compared to the mechanistic model (Figure 9).

FIGURE 9

Figure 9. Comparison between the Mechanistic and IO model runtime vs. number of spines. Each spine represents the kinetics for synaptic transmission as well as calcium dynamics, which is either represented by the mechanistic calcium model or the IO calcium model. Benchmarking was conducted at 5, 10, 15, 20, 30, 40, and 50 spines, and simulations were repeated 10 times to derive the standard deviation of the runtimes. For each simulation, a 2 Hz Poisson random interval train input was given to all synapses. Overall, the IO spine model required on average a little less than half the runtime needed to finish a simulation compared to the mechanistic model.

Another advantage of the IO calcium model is that the framework of the input-output model is easily implemented and adaptable to other simulation platforms. To test the performance of the IO calcium model simulated as an embedded mechanism within the NEURON engine, we adapted the IO model into a module file for the NEURON platform and compared cell level simulations based on number of spine instances. Two types of models were simulated with different spine configurations: (1) spines using both an NMDA 8 state model and the IO calcium model, and (2) spines with only the NMDA 8 state model and no IO calcium model. Simulation protocols with 10, 100, 500, 1,000, 5,000, and 10,000 spine instances were conducted; simulations were run in fixed time step of 0.1 ms, a randomized poisson input train of 2 Hz frequency, and the overall simulated time period is for 20 s. Simulation times were then benchmarked to determine how much of a computational burden is added when including IO calcium model within spines. Simulations were repeated 10 times each to derive the standard deviation in the simulation time. Results are shown in Figure 10. Our simulations concluded that at 10,000 spines, the computational burden increases from 9.3 to 13.3 h. For details on the variation in the benchmarks, the standard deviation of the benchmarking data from Figures 9, 10 are plotted in Supplementary Figure 4. While the increased required simulation time is not insignificant, the IO framework still gives a viable option for simulating complex postsynaptic calcium dynamics on a larger scale—with numerous spines on a neuron or a neuron network.

FIGURE 10

Figure 10. Benchmarking the IO model within the NEURON simulation platform. The IO model was re-implemented as a module file usable within the NEURON platform and simulated on a compartmental neuron model for 20 simulated seconds at 0.1 ms timesteps. The number of spines was varied for each simulation ranging from 10 to 10,000, where a spine was defined as a kinetic NMDA 8 state model, either with or without the IO calcium model and the simulation runtime is plotted here. Each simulation was repeated 10 times each to derive the standard deviation, shown as error bars in the figure.

Discussion

This article describes the development and simulation of a model of the postsynaptic calcium concentration in the spine. The model presented is an integration of various mechanisms which shape the dynamics of calcium concentration at the postsynaptic spine, comprising elements that contribute to calcium influx, calcium extrusion, and buffering. Experimental studies on spine signaling have focused on calcium more than any other signaling molecule within active spines (Higley and Sabatini, 2012), as calcium dynamics and its effectors (NMDA, VDCC, etc.) have been repeatedly shown to strongly influence plasticity and learning. Calcium has also been implicated as a role player in neurodegenerative diseases such as Alzheimer's Disease (Khachaturian, 1994; Alberdi et al., 2010). Our goal in modeling calcium is to: (1) explore mechanisms and details underlying calcium dynamics that would otherwise be difficult to achieve with experimental studies alone (i.e., influence of pre-post timing and distance on spine calcium, where the researcher must take multiple time measurements at multiple spine locations, would be difficult to measure in experimental setups), and (2) reduce computational complexity of the calcium model to enable multi-scale simulations. We have presented a viable model which is supported by experimental data. We configured the model to replicate a particular type of synapse—a glutamatergic CA3-CA1 synapse of a pyramidal cell neuron; it incorporates the nonlinear dynamics that result from interactions between the components that contribute to spine calcium concentration.

Beyond experimental validation, we presented simulated results with the mechanistic model which show changes in spine calcium activity as a function of presynaptic and postsynaptic intervals—a standard protocol for inducing spike-timing dependent plasticity. In STDP, the weight of a synapse changes after repeated identical pre-post stimulations at given pre-post intervals. In hippocampal CA1 glutamatergic spines, intervals where presynaptic stimulation precedes postsynaptic stimulation induce synaptic potentiation, with the strength of the potentiation inversely proportional to the interval distance between the pre- and post-stimulation (Bi and Poo, 2001). Our model demonstrates that similarly, calcium influx is significantly amplified when presynaptic stimulation precedes postsynaptic stimulation, and that the amplitude is also inversely proportional to the interval size. Many plasticity associated signaling cascades are activated by calcium—for example, AMPA receptor upregulation into the spine is a known indicator for synaptic strengthening (Zhabotinsky et al., 2006). This process is initiated by spine calcium binding with CaMKII and triggering secondary messenger pathways. Likewise, recruitment of actin will lead larger spines—this, too has been associated with calcium interaction (Araya, 2014).

It is thought that the major influence on the calcium brought about by presynaptic/postsynaptic interactions is the NMDA receptor channel kinetics, but we have demonstrated in our simulations that the NMDAr channel alone is not the sole contributor of the nonlinear dynamics of calcium in the spine. The role of NMDAr in synaptic activity is considerably important: it has been shown experimentally that NMDA contributes to synaptic plasticity and LTP (Sakimura et al., 1995; Grover et al., 2009; Larson and Munkácsy, 2015). Many LTP models are based around this hypothesis, where NMDAr models are used to represent calcium influx, and repeated stimulation leads to calcium induced plasticity (Shouval et al., 2002; Standage et al., 2014). However, the NMDA representation in such models is a simple, linear representation where the NMDA-based calcium influx is represented as a ratio proportional to the bAP. This poorly reflects on nonlinear calcium dynamics in two ways: (1) the simplified version of NMDAr dynamics ignores important dynamical features that are known to be associated in NMDAr channels, such as desensitization (Mayer et al., 1989), which is included as a state in the NMDAr kinetic model used in our mechanistic calcium platform, and the magnesium blockade (Calabresi et al., 1992) which is instead roughly approximated using a BPAP curve; (2) there is no influence or contribution from other elements or properties from the spine, which can drastically alter the calcium response. Meanwhile, our model integrates validated channel kinetics within the confines of the spine compartment. Thus, our mechanistic model can consider the nonlinear aspects of calcium influx which are influenced by NMDAr channels, along with other channels, pumps, and buffers that regulate spine calcium concentration—all of which influence observed calcium levels at the spine. Especially to note is the spine volume and the buffers within the spine. Changes to volume can result in undercompensation (in larger volumes) or overcompensation (in smaller volumes) of calcium concentration, unless the mechanisms which govern calcium dynamics are scaled appropriately (O'Donnell et al., 2011). Furthermore, experimental evidence has shown that changes in volume also can result in change in AMPAr expression and upregulation (Noguchi et al., 2011), which is also important for synaptic plasticity. Such mechanisms revolved around volume are ignored in linear calcium concentration models, but can be studied in future simulations using the mechanistic model. Buffering provides another layer of complexity that will influence the amplitude and decay of the calcium response, yet there is no consideration for buffers in the linear calcium models. Our results demonstrate that, in the study of spine calcium, these components are important contributors to the nonlinear response of calcium.

The investigation into calcium dynamics at the postsynaptic spine is our latest study in the use of computational synaptic modeling platform to better understand synaptic activity and signaling: (1) We modeled the impact of astrocytic glutamate uptake (Allam et al., 2012) and ionotropic receptor distribution (Allam et al., 2015) on synaptic transmission in glutamatergic CA1 synapses; (2) Our synapse platform has been adapted to cellular and network levels in simulation (Bouteiller et al., 2011) to observe effects of nonlinear activity of synapses in a network simulation; (3) Large scale simulation models containing millions of neurons have also been developed within our research lab (Hendrickson et al., 2015); (4) and efforts had been made to adapt the complex nonlinear postsynaptic response of mechanistic synapses to large scale simulations using input-output modeling (Hu et al., 2015). Our modeling platform is consistently under expansion, with current projects considering the effects of modulators such as acetylcholine and how intracellular calcium stores influence metabolism and pump activity. The mechanistic and IO calcium model we describe in this manuscript expands our modeling platform to simulate not only synaptic transmission, but complex calcium dynamics as well. From here on, we plan to investigate and implement mechanisms that are based on spine calcium (the CaMKII signaling pathway, for example) and move to the next hierarchical level of calcium dynamics, the calcium response at the dendrite.

Expansion into large-scale, multi-scale modeling with complex biologically accurate synapse dynamics requires reduction of the computational burden while minimizing loss in accuracy. Spine calcium plays a key role in synaptic plasticity and influences communication between neurons, and understanding how calcium dynamics change network properties on a large scale can give us a better sense of the mechanisms that give rise to plasticity. Key downstream processes are influenced by the slightest changes in calcium dynamics (timing, magnitude, frequency, decay) (Evans and Blackwell, 2015). For example, an increase in spine calcium levels can activate signaling cascades that lead to either LTP or LTD induction (Lisman, 1989; Artola et al., 1990; Malenka and Bear, 2004); however, it is also observed in experiments that there is not a simple threshold that distinguishes when LTP or LTD occurs during calcium influx (Neveu and Zucker, 1996)—emphasizing even more the need of an integrated model of varying calcium dynamics, not just a linear model based on thresholds. Hence, we believe accurate representations of nonlinear calcium are required not only at the subcellular scale models, but on larger network-level models as well. Our model can provide an accurate reflection on the magnitude, duration, and location of spine calcium response—nonlinear dynamics that have been described to be more and more important for calcium based synaptic plasticity (Evans and Blackwell, 2015). Furthermore, neurodegenerative disease are often accompanied by an imbalance in calcium levels (Arundine and Tymianski, 2003). Nonlinear calcium models can be modified to represent pathological conditions, and multi-scale modeling can help identify the network level changes that occur with neurodegeneration and disease.

However, the computational cost of using many instances of the mechanistic calcium model in full exceeds the computational capacity of even the most recent high-performance computers; a method to improve computational efficiency is needed. Our previous work considers the use of the Volterra functional series to develop an IO model for the postsynaptic response to a synapse (Hu et al., 2015). Using the same input-output framework, we adapted this method to reduce the computational burden of modeling calcium dynamics. We have shown that using the IO model reduces the required simulation time by two to three-fold compared to the mechanistic model within our MEMORY platform. We also demonstrate that the IO framework can be easily adapted into other platforms such as NEURON, where, in our setup, 10,000 instances of the IO calcium model with the NMDAr kinetic model can be simulated on a complex, morphological, compartmental cell model, resulting in 1 h additional simulation time compared to the same protocol using the NMDAr kinetic model but without the IO calcium model. However, the current IO model as described has a limitation that must be addressed: it is limited to a 3rd order model, with higher order models requiring exponentially increasing memory requisites. Spine calcium dynamics become progressively more nonlinear when given higher frequency input, such as high frequency stimulation protocols often used in LTP induction. As such, other IO model frameworks are being investigated, such as the Laguerre-Volterra network structure (Geng and Marmarelis, 2016), as possible enhancements leading to even more efficient computational modeling of complex dynamic systems.

The concept of using computer simulations to study postsynaptic calcium dynamics is not new. There have been several computational models of spine calcium that have been developed previously, and their work has provided useful insights on the dynamics and importance of calcium (Shouval et al., 2002; Standage et al., 2014; Bartol et al., 2015). Generally, current calcium models are divided into two categories: (1) Phenomenological models which describe very few mechanistic aspects of calcium dynamics, but help understand its influence on synaptic plasticity and LTP (Shouval et al., 2002; Naoki et al., 2005; Zhabotinsky et al., 2006; Standage et al., 2014); and (2) detailed, stochastic models which describe calcium all the way down to each individual ion (Bartol et al., 2015). In (1), the models do not extensively consider calcium dynamics at length and may even consider calcium as a linear system. Instead, models from (1) evaluate the downstream effects of calcium on important synaptic processes such as plasticity. In contrast, the detailed calcium model in (2) considers a complete reconstruction of a small area on the CA1 pyramidal cell neuron, including specific channel densities and exact volume and shape reconstruction of a 6 × 6 × 5 μm³ cube of neuropil in a Monte Carlo based stochastic simulator. However, such a model is difficult to adapt outside of the scope of the reconstructed area and is computationally intensive, making it unsuitable for larger scale simulations. The models in this article help bridge this discrepancy: (1) the mechanistic model is capable of replicating complex non-linear interactions between the elements that shape spine calcium dynamics, and (2) the input-output model provides a method to simulate these complex calcium dynamics on a larger scale.

The spine is a constantly changing organelle as a result of development, plasticity, and/or pathological conditions. The current spine calcium model as described here represents only a snapshot of a particular spine, constrained by static parameters based on the averaged responses from experimental data. Future renditions of our model will not be limited to a single type of synapse. Our spine calcium model has potential to be adapted to varying physiological states (i.e., different morphologies and channel distributions) and pathological conditions (such as Alzheimer's disease). Furthermore, the model has potential applications in drug discovery, for in silico testing of compounds that modulate calcium (either directly or via channel/pump interactions). The calcium model is an expansion of the synapse model framework that is constantly being built upon to provide extensive and detailed multi-level models that can help explore the pathways and processes of that take place in the spine and influence synaptic plasticity, neuron communication, and pathological processes.

Availability

Model scripts and code will be made available in the future on the EONS synaptic platform modeling site, www.synapticmodeling.com.

Author Contributions

EH and J-MB: concept and design of study; EH and AM: data acquisition; EH, DS, and J-MB: analysis and/or interpretation of data; EH and J-MB: drafting of the manuscript; CB, DS, J-MB, and TB: critical revision; EH, AM, CB, DS, J-MB, and TB: approval of the manuscript to be published.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Work supported by NIBIB grant U01 GM104604 to TB.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2018.00058/full#supplementary-material

References

Alberdi, E., Sánchez-Gómez, M. V., Cavaliere, F., Pérez-Samartín, A., Zugaza, J. L., Trullas, R., et al. (2010). Amyloid beta oligomers induce Ca2+ dysregulation and neuronal death through activation of ionotropic glutamate receptors. Cell Calcium 47, 264–272. doi: 10.1016/j.ceca.2009.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Allam, S. L., Bouteiller, J. M. C., Hu, E. Y., Ambert, N., Greget, R., Bischoff, S., et al. (2015). Synaptic efficacy as a function of ionotropic receptor distribution: a computational study. PLoS ONE 10:e0140333. doi: 10.1371/journal.pone.0140333

PubMed Abstract | CrossRef Full Text | Google Scholar

Allam, S. L., Ghaderi, V. S., Bouteiller, J. M. C., Legendre, A., Ambert, N., Greget, R., et al. (2012). A computational model to investigate astrocytic glutamate uptake influence on synaptic transmission and neuronal spiking. Front. Comput. Neurosci. 6:70. doi: 10.3389/fncom.2012.00070

PubMed Abstract | CrossRef Full Text | Google Scholar

Ambert, N., Greget, R., Haeberlé, O., Bischoff, S., Berger, T. W., Bouteiller, J. M., et al. (2010). Computational studies of NMDA receptors: differential effects of neuronal activity on efficacy of competitive and non-competitive antagonists. Open Access Bioinformatics 2, 113–125. doi: 10.2147/OAB.S7246

PubMed Abstract | CrossRef Full Text | Google Scholar

Araya, R. (2014). Input transformation by dendritic spines of pyramidal neurons. Front. Neuroanat. 8:141. doi: 10.3389/fnana.2014.00141

PubMed Abstract | CrossRef Full Text | Google Scholar

Artola, A., Bröcher, S., and Singer, W. (1990). Different voltage-dependent thresholds for inducing long-term depression and long-term potentiation in slices of rat visual cortex. Nature 347, 69–72. doi: 10.1038/347069a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Arundine, M., and Tymianski, M. (2003). Molecular mechanisms of calcium-dependent neurodegeneration in excitotoxicity. Cell Calcium. 34, 325–337 doi: 10.1016/S0143-4160(03)00141-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bading, H., Ginty, D., and Greenberg, M. (1993). Regulation of gene expression in hippocampal neurons by distinct calcium signaling pathways. Science 260, 181–186. doi: 10.1126/science.8097060

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartol, T. M., Keller, D. X., Kinney, J. P., Bajaj, C. L., Harris, K. M., Sejnowski, T. J., et al. (2015). Computational reconstitution of spine calcium transients from individual proteins. Front. Synaptic Neurosci. 7:17. doi: 10.3389/fnsyn.2015.00017

PubMed Abstract | CrossRef Full Text | Google Scholar

Berger, T., and Song, D. (2010). “The neurobiological basis of cognition: identification by multi-input, multioutput nonlinear dynamic modeling,” in Proceedings of the IEEE Institute of Electrical and Electronics Engineers, Vol. 98. Available online at: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=5424199

Berger, T. W., Eriksson, J. L., Ciarolla, D. A., and Sclabassi, R. J. (1988). Nonlinear systems analysis of the hippocampal perforant path-dentate projection. ii. effects of random impulse train stimulation. J. Neurophysiol. 60, 1077–1094. doi: 10.1152/jn.1988.60.3.1077

PubMed Abstract | CrossRef Full Text

Berger, T. W., Song, D., Chan, R. H. M., and Marmarelis, V. Z. (2010). The neurobiological basis of cognition: identification by multi-input, multioutput nonlinear dynamic modeling: a method is proposed for measuring and modeling human long-term memory formation by mathematical analysis and computer simulation of nerve-cell. Proc. IEEE Inst. Elect. Electr. Eng. 98, 356–374. doi: 10.1109/JPROC.2009.2038804

PubMed Abstract | CrossRef Full Text | Google Scholar

Bharathy, C., Sachdeva, P., Parthasarthy, H., and Tayal, A. (2008). “Advanced intelligent computing theories and applications. with aspects of contemporary intelligent computing techniques,” in An Introduction to Volterra Series and Its Application on Mechanical Systems, Communications in Computer and Information Science, Vol. 15, eds D.-S. Huang, D. C. Wunsch, D. S. Levine, and K.-H. Jo (Berlin; Heidelberg: Springer), 478.

Google Scholar

Bi, G., and Poo, M. (2001). Synaptic modification by correlated activity: hebb's postulate revisited. Annu. Rev. Neurosci. 24, 139–66. doi: 10.1146/annurev.neuro.24.1.139

PubMed Abstract | CrossRef Full Text | Google Scholar

Bikbaev, A., Duménieu, M., Lopez-Rojas, J., and Heine, M. (2016). “Localising receptors and channels across the dendritic arbour,” in Dendrites, eds K. Emoto, R. Wong, E. Huang, and C. Hoogenraad (Tokyo: Springer), 387–424.

Google Scholar

Bloodgood, B. L., and Sabatini, B. L. (2007). Nonlinear Regulation of unitary synaptic signals by CaV(2.3) voltage-sensitive calcium channels located in dendritic spines. Neuron 53, 249–260. doi: 10.1016/j.neuron.2006.12.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloodgood, B. L., and Sabatini, B. L. (2009). NMDA Receptor-Mediated Calcium Transients in Dendritic Spines. Boca Raton, FL: CRC Press.

Google Scholar

Bouteiller, J. M. C., Allam, S. L., Hu, E. Y., Greget, R., Ambert, N., Keller, A. F., et al. (2011). Integrated multiscale modeling of the nervous system: predicting changes in hippocampal network activity by a positive AMPA receptor modulator. IEEE Trans. BioMed. Eng. 58, 3008–3011. doi: 10.1109/TBME.2011.2158605

PubMed Abstract | CrossRef Full Text | Google Scholar

Burette, A. C., Strehler, E. E., and Weinberg, R. J. (2010). A plasma membrane Ca2+ ATPase isoform at the postsynaptic density. Neuroscience 169, 987–993. doi: 10.1016/j.neuroscience.2010.05.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Burnashev, N., Zhou, Z., Neher, E., and Sakmann, B. (1995). Fractional calcium currents through recombinant GluR channels of the NMDA, AMPA and kainate receptor subtypes. J. Physiol. 485(Pt 2), 403–418. doi: 10.1113/jphysiol.1995.sp020738

PubMed Abstract | CrossRef Full Text | Google Scholar

Calabresi, P., Pisani, A., Mercuri, N. B., and Bernardi, G. (1992). Long-term potentiation in the striatum is unmasked by removing the voltage-dependent magnesium block of NMDA receptor channels. Eur. J. Neurosci. 4, 929–935. doi: 10.1111/j.1460-9568.1992.tb00119.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Carafoli, E. (1991). Calcium pump of the plasma membrane. Physiol. Rev. 71, 129–153. doi: 10.1152/physrev.1991.71.1.129

PubMed Abstract | CrossRef Full Text | Google Scholar

Carafoli, E., Santella, L., Branca, D., and Brini, M. (2001). Generation, control, and processing of cellular calcium signals. Crit. Rev. Biochem. Mol. Biol. 36, 107–260. doi: 10.1080/20014091074183

PubMed Abstract | CrossRef Full Text | Google Scholar

Carnevale, N. T., and Hines, M. L. (2006). The NEURON Book. Cambridge University Press. Available online at: http://books.google.com/books?hl=ja&lr=&id=YzcOyjKBPHgC&pgis=1

Google Scholar

Catterall, W. A. (2011). Voltage-Gated calcium channels. Cold Spring Harb. Perspect. Biol. 3:a003947. doi: 10.1101/cshperspect.a003947

PubMed Abstract | CrossRef Full Text | Google Scholar

Contreras, L., Drago, I., Zampese, E., and Pozzan, T. (2010). Mitochondria: the calcium connection. Biochim. Biophys. Acta 1797, 607–618. doi: 10.1016/j.bbabio.2010.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolmetsch, R. E., Xu, K., and Lewis, R. S. (1998). Calcium oscillations increase the efficiency and specificity of gene expression. Nature 392, 933–936. doi: 10.1038/31960

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebashi, S., and Endo, M. (1968). Calcium ion and muscle contraction. Prog. Biophys. Mol. Biol. 18, 123–183. doi: 10.1016/0079-6107(68)90023-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Edelmann, E., Cepeda-Prado, E., and Leßmann, V. (2017). Coexistence of multiple types of synaptic plasticity in individual hippocampal CA1 pyramidal neurons. Front. Synaptic Neurosci. 9:7. doi: 10.3389/fnsyn.2017.00007

PubMed Abstract | CrossRef Full Text | Google Scholar

Emptage, N. J., Reid, C. A., and Fine, A. (2001). Calcium stores in hippocampal synaptic boutons mediate short-term plasticity, store-operated Ca2+ entry, and spontaneous transmitter release. Neuron 29, 197–208. doi: 10.1016/S0896-6273(01)00190-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Erreger, K., Geballe, M. T., Dravid, S. M., Snyder, J. P., Wyllie, D. J. A., and Traynelis, S. F. (2005). Mechanism of partial agonism at NMDA receptors for a conformationally restricted glutamate analog. J. Neurosci. 25, 7858–7866. doi: 10.1523/JNEUROSCI.1613-05.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, R. C., and Blackwell, K. T. (2015). Calcium: amplitude, duration, or location? Biol. Bull. 228, 75–83. doi: 10.1086/BBLv228n1p75

CrossRef Full Text | Google Scholar

Fridlyand, L. E., Tamarina, N., and Philipson, L. H. (2003). Modeling of Ca2+ flux in pancreatic beta-cells: role of the plasma membrane and intracellular stores. Am. J. Physiol. Endocrinol. Metab. 285, E138–E154. doi: 10.1152/ajpendo.00194.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Funahashi, A., Matsuoka, Y., Jouraku, A., Morohashi, M., Kikuchi, N., and Kitano, H. (2008). CellDesigner 3.5: a versatile modeling tool for biochemical networks. Proc. IEEE 96, 1254–1265. doi: 10.1109/JPROC.2008.925458

CrossRef Full Text | Google Scholar

Gall, D., and Susa, I. (1999). Effect of Na/Ca exchange on plateau fraction and [Ca]i in models for bursting in pancreatic beta-cells. Biophys. J. 77, 45–53. doi: 10.1016/S0006-3495(99)76871-2

CrossRef Full Text | Google Scholar

Geng, K., and Marmarelis, V. Z. (2016). Methodology of recurrent laguerre-volterra network for modeling nonlinear dynamic systems. IEEE Trans. Neural Netw. Learn. Syst. 28, 1–13. doi: 10.1109/TNNLS.2016.2581141

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghaderi, V. S., Allam, S. L., Ambert, N., Bouteiller, J. M. C., Choma, J., and Berger, T. W. (2011). “Modeling neuron-glia interactions: from parametric model to neuromorphic hardware,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Annual Conference (Boston, MA), 3581–3584.

PubMed Abstract | Google Scholar

Golding, N. L., Kath, W. L., and Spruston, N. (2001). Dichotomy of action-potential backpropagation in CA1 pyramidal neuron dendrites. J. Neurophysiol. 86, 2998–3010. doi: 10.1152/jn.2001.86.6.2998

PubMed Abstract | CrossRef Full Text | Google Scholar

Grover, L. M., Kim, E., Cooke, J. D., and Holmes, W. R. (2009). LTP in hippocampal area CA1 is induced by burst stimulation over a broad frequency range centered around delta. Learn. Mem. 16, 69–81. doi: 10.1101/lm.1179109

PubMed Abstract | CrossRef Full Text | Google Scholar

Grunditz, A., Holbro, N., Tian, L., Zuo, Y., and Oertner, T. G. (2008). Spine neck plasticity controls postsynaptic calcium signals through electrical compartmentalization. J. Neurosci. 28, 13457–13466. doi: 10.1523/JNEUROSCI.2702-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Hendrickson, P. J., Yu, G. J., Song, D., and Berger, T. W. (2015). “A million-plus neuron model of the hippocampal dentate gyrus: dependency of spatio-temporal network dynamics on topography,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Annual Conference 2015 (Milan: NIH Public Access), 4713–4716.

PubMed Abstract | Google Scholar

Higgins, E. R., Cannell, M. B., and Sneyd, J. (2006). A buffering SERCA pump in models of calcium dynamics. Biophys. J. 91, 151–163. doi: 10.1529/biophysj.105.075747

PubMed Abstract | CrossRef Full Text | Google Scholar

Higley, M. J., and Sabatini, B. L. (2008). Calcium signaling in dendrites and spines: practical and functional considerations. Neuron 59, 902–913. doi: 10.1016/j.neuron.2008.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Higley, M. J., and Sabatini, B. L. (2012). Calcium signaling in dendritic spines. Cold Spring Harb. Perspect. Biol. 4:a005686. doi: 10.1101/cshperspect.a005686

PubMed Abstract | CrossRef Full Text | Google Scholar

Holbro, N., Grunditz, A., and Oertner, T. G. (2009). Differential distribution of endoplasmic reticulum controls metabotropic signaling and plasticity at hippocampal synapses. Proc. Natl. Acad. Sci. 106, 15055–15060. doi: 10.1073/pnas.0905110106

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, E. Y., Bouteiller, J. M. C., Song, D., Baudry, M., and Berger, T. W. (2015). Volterra representation enables modeling of complex synaptic nonlinear dynamics in large-scale simulations. Front. Comput. Neurosci. 9:112. doi: 10.3389/fncom.2015.00112

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, E. Y., Bouteiller, J. M. C., Song, D., and Berger, T. W. (2016). “Development of a detailed model of calcium dynamics at the postsynaptic spine of an excitatory synapse,” in 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (Orlando, FL: IEEE), 6102–6105.

Google Scholar

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

PubMed Abstract | CrossRef Full Text | Google Scholar

Jaffe, D. B., Ross, W. N., Lisman, J. E., Lasser-Ross, N., Miyakawa, H., and Johnston, D. (1994). A model for dendritic Ca2+ accumulation in hippocampal pyramidal neurons based on fluorescence imaging measurements. J. Neurophysiol. 71, 1065–1077. doi: 10.1152/jn.1994.71.3.1065

PubMed Abstract | CrossRef Full Text | Google Scholar

Jahr, C. E., and Stevens, C. F. (1990). Voltage dependence of NMDA-activated macroscopic conductances predicted by single-channel kinetics. J. Neurosci. 10, 3178–3182. doi: 10.1523/JNEUROSCI.10-09-03178.1990

PubMed Abstract | CrossRef Full Text | Google Scholar

Jarsky, T., Roxin, A., Kath, W. L., and Spruston, N. (2005). Conditional dendritic spike propagation following distal synaptic activation of hippocampal CA1 pyramidal neurons. Nat. Neurosci. 8, 1667–1676. doi: 10.1038/nn1599

PubMed Abstract | CrossRef Full Text | Google Scholar

Khachaturian, Z. S. (1994). Calcium hypothesis of Alzheimer's disease and brain aging. Ann. N. Y. Acad. Sci. 747, 1–11.

PubMed Abstract | Google Scholar

Larson, J., and Munkácsy, E. (2015). Theta-burst LTP. Brain Res. 1621, 38–50. doi: 10.1016/j.brainres.2014.10.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Okamoto, K. I., Hayashi, Y., and Sheng, M. (2004). The importance of dendritic mitochondria in the morphogenesis and plasticity of spines and synapses. Cell 119, 873–887. doi: 10.1016/j.cell.2004.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Lisman, J. (1989). A mechanism for the hebb and the anti-hebb processes underlying learning and memory. Proc. Natl. Acad. Sci. U.S.A. 86, 9574–9578. doi: 10.1073/pnas.86.23.9574

PubMed Abstract | CrossRef Full Text | Google Scholar

Lörincz, A., Rózsa, B., Katona, G., Vizi, E. S., and Tamás, G. (2007). Differential distribution of NCX1 contributes to spine-dendrite compartmentalization in CA1 pyramidal cells. Proc. Natl. Acad. Sci. U.S.A. 104, 1033–1038. doi: 10.1073/pnas.0605412104

PubMed Abstract | CrossRef Full Text | Google Scholar

Malenka, R. C., and Bear, M. F. (2004). LTP and LTD: an embarrassment of riches. Neuron 44, 5–21. doi: 10.1016/j.neuron.2004.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Markram, H. (2006). The blue brain project. Nat. Rev. Neurosci. 7, 153–160. doi: 10.1038/nrn1848

PubMed Abstract | CrossRef Full Text | Google Scholar

Marmarelis, P., and Marmarelis, V. (1978). Analysis of Physiological Systems: The White-Noise Approach. New York, NY: Plenum. doi: 10.1007/978-1-4613-3970-0

CrossRef Full Text | Google Scholar

Mattson, M. P., and Chan, S. L. (2003). Calcium orchestrates apoptosis. Nat. Cell Biol. 5, 1041–1043. doi: 10.1038/ncb1203-1041

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayer, M. L., Vyklicky, L., and Clements, J. (1989). Regulation of NMDA receptor desensitization in mouse hippocampal neurons by glycine. Nature 338, 425–427. doi: 10.1038/338425a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Megías, M., Emri, Z., Freund, T., and Gulyás, A. (2001). Total number and distribution of inhibitory and excitatory synapses on hippocampal CA1 pyramidal cells. Neuroscience 102, 527–540. doi: 10.1016/S0306-4522(00)00496-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Migliore, M., and Migliore, R. (2012). Know your current I(h): interaction with a shunting current explains the puzzling effects of its pharmacological or pathological modulations. PLoS ONE 7:e36867. doi: 10.1371/journal.pone.0036867

CrossRef Full Text | Google Scholar

Naoki, H., Sakumura, Y., and Ishii, S. (2005). Local signaling with molecular diffusion as a decoder of Ca2+ signals in synaptic plasticity. Mol. Syst. Biol. 1:2005.0027. doi: 10.1038/msb4100035

PubMed Abstract | CrossRef Full Text | Google Scholar

Neveu, D., and Zucker, R. S. (1996). Postsynaptic levels of [Ca2+]i needed to trigger LTD and LTP. Neuron 16, 619–629. doi: 10.1016/S0896-6273(00)80081-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Noguchi, J., Nagaoka, A., Watanabe, S., Ellis-Davies, G. C. R., Kitamura, K., Kano, M., et al. (2011). In vivo two-photon uncaging of glutamate revealing the structure-function relationships of dendritic spines in the neocortex of adult mice. J. Physiol. 589(Pt 10):2447–2457. doi: 10.1113/jphysiol.2011.207100

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Donnell, C., Nolan, M. F., and Van Rossum, M. C. W. (2011). Dendritic spine dynamics regulate the long-term stability of synaptic plasticity. J. Neurosci. 31, 16142–16156. doi: 10.1523/JNEUROSCI.2520-11.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Paula-Lima, A. C., Adasme, T., and Hidalgo, C. (2014). Contribution of Ca ²⁺ release channels to hippocampal synaptic plasticity and spatial memory: potential redox modulation. Antioxid. Redox Signal. 21, 892–914. doi: 10.1089/ars.2013.5796

CrossRef Full Text | Google Scholar

Perez-Reyes, E., Cribbs, L. L., Daud, A., Lacerda, A. E., Barclay, J., Williamson, M. P., et al. (1998). Molecular characterization of a neuronal low-voltage-activated T-type calcium channel. Nature 391, 896–900. doi: 10.1038/36110

PubMed Abstract | CrossRef Full Text | Google Scholar

Poirazi, P., Brannon, T., and Mel, B. W. (2003). Arithmetic of subthreshold synaptic summation in a model CA1 pyramidal cell. Neuron 37, 977–987. doi: 10.1016/S0896-6273(03)00148-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Racca, C., Stephenson, F. A., Streit, P., Roberts, J. D. B., and Somogyi, P. (2000). NMDA receptor content of synapses in stratum radiatum of the hippocampal ca1 area. J. Neurosci. 20, 2512–2522. doi: 10.1523/JNEUROSCI.20-07-02512.2000

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranck, J. B. (1973). Studies on single neurons in dorsal hippocampal formation and septum in unrestrained rats. Exp. Neurol. 41, 462–531. doi: 10.1016/0014-4886(73)90290-2

CrossRef Full Text | Google Scholar

Robert, A., and Howe, J. R. (2003). How AMPA receptor desensitization depends on receptor occupancy. J. Neurosci. 23, 847–858. doi: 10.1523/JNEUROSCI.23-03-00847.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sabatini, B. L., Oertner, T. G., and Svoboda, K. (2002). The life cycle of Ca2+ ions in dendritic spines. Neuron 33, 439–452. doi: 10.1016/S0896-6273(02)00573-1

CrossRef Full Text | Google Scholar

Sabatini, B. L., and Svoboda, K. (2000). Analysis of calcium channels in single spines using optical fluctuation analysis. Nature 408, 589–593. doi: 10.1038/35046076

PubMed Abstract | CrossRef Full Text | Google Scholar

Sakimura, K., Kutsuwada, T., Ito, I., Manabe, T., Takayama, C., Kushiya, E., et al. (1995). Reduced hippocampal LTP and spatial learning in mice lacking NMDA receptor epsilon 1 subunit. Nature 373, 151–155. doi: 10.1038/373151a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheuss, V., Yasuda, R., Sobczyk, A., and Svoboda, K. (2006). Nonlinear [Ca2+] signaling in dendrites and spines caused by activity-dependent depression of Ca2+ extrusion. J. Neurosci. 26, 8183–8194. doi: 10.1523/JNEUROSCI.1962-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Seo, H., and Jun, S. C. (2017). Multi-scale computational models for electrical brain stimulation. Front. Hum. Neurosci. 11:515. doi: 10.3389/fnhum.2017.00515

PubMed Abstract | CrossRef Full Text | Google Scholar

Shouval, H. Z., Bear, M. F., and Cooper, L. N. (2002). A unified model of NMDA receptor-dependent bidirectional synaptic plasticity. Proc. Natl. Acad. Sci. U.S.A. 99, 10831–10836. doi: 10.1073/pnas.152343099

PubMed Abstract | CrossRef Full Text | Google Scholar

Sjöström, P. J., and Häusser, M. (2006). A cooperative switch determines the sign of synaptic plasticity in distal dendrites of neocortical pyramidal neurons. Neuron 51, 227–238. doi: 10.1016/j.neuron.2006.06.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Somogyi, E. T., Bouteiller, J. M., Glazier, J. A., König, M., Medley, J. K., Swat, M. H., et al. (2015). LibRoadRunner: a high performance SBML simulation and analysis library: table 1. Bioinformatics 31, 3315–3321. doi: 10.1093/bioinformatics/btv363

CrossRef Full Text | Google Scholar

Song, D., Marmarelis, V. Z., and Berger, T. W. (2009a). Parametric and non-parametric modeling of short-term synaptic plasticity. Part I: Computational Study. J. Comput. Neurosci. 26, 1–19. doi: 10.1007/s10827-008-0097-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, D., Wang, Z., Marmarelis, V. Z., and Berger, T. W. (2009b). Parametric and non-parametric modeling of short-term synaptic plasticity. Part II: experimental study. J. Comput. Neurosci. 26, 21–37. doi: 10.1007/s10827-008-0098-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Standage, D., Trappenberg, T., Blohm, G., Bear, M., and Huganir, R. (2014). Calcium-dependent calcium decay explains STDP in a dynamic model of hippocampal synapses. edited by miguel maravall. PLoS ONE 9:e86248. doi: 10.1371/journal.pone.0086248

CrossRef Full Text | Google Scholar

Stewart, M. G., Medvedev, N. I., Popov, V. I., Schoepfer, R., Davies, H. A., Murphy, K., et al. (2005). Chemically induced long-term potentiation increases the number of perforated and complex postsynaptic densities but does not alter dendritic spine volume in CA1 of adult mouse hippocampal slices. Eur. J. Neurosci. 21, 3368–3378 doi: 10.1111/j.1460-9568.2005.04174.x

CrossRef Full Text | Google Scholar

Tu, C., Song, D., Breidt, F., Berger, T., and Wang, H. (2012). “Functional model selection for sparse binary time series with multiple inputs,” in Economic Time Series: Modeling and Seasonality, eds R. William Bell and H. Scott Holan (Boca Raton, FL: Chapman and Hall/CRC), 477–498.

Google Scholar

Weber, A., and Murray, J. M. (1973). Molecular control mechanisms in muscle contraction. Physiol. Rev. 53, 612–673. doi: 10.1152/physrev.1973.53.3.612

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, G. S. B., Chikando, A. C., Tuan, H. T. M., Sobie, E. A., Lederer, W. J., and Jafri, M. S. (2011). Dynamics of calcium sparks and calcium leak in the heart. Biophys. J. 101, 1287–1296. doi: 10.1016/j.bpj.2011.07.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J. S., and Bagheri, N. (2016). Multi-class and multi-scale models of complex biological phenomena. Curr. Opin. Biotechnol. 39, 167–173. doi: 10.1016/j.copbio.2016.04.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhabotinsky, A. M., Camp, R. N., Epstein, I. R., and Lisman, J. E. (2006). Role of the neurogranin concentrated in spines in the induction of long-term potentiation. J. Neurosci. 26, 7337–7347. doi: 10.1523/JNEUROSCI.0729-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, K., and Prince, R. L. (2012). Calcium and bone. Clin. Biochem. 45, 936–942. doi: 10.1016/j.clinbiochem.2012.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Zucker, R. S. (1999). Calcium- and activity-dependent synaptic plasticity. Curr. Opin. Neurobiol. 9, 305–313. doi: 10.1016/S0959-4388(99)80045-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: multi-scale modeling, spine calcium, computational model, nonlinear dynamical systems, glutamatergic synapse

Citation: Hu E, Mergenthal A, Bingham CS, Song D, Bouteiller J-M and Berger TW (2018) A Glutamatergic Spine Model to Enable Multi-Scale Modeling of Nonlinear Calcium Dynamics. Front. Comput. Neurosci. 12:58. doi: 10.3389/fncom.2018.00058

Received: 28 December 2017; Accepted: 05 July 2018;
Published: 27 July 2018.

Edited by:

Florentin Wörgötter, Georg-August-Universität Göttingen, Germany

Reviewed by:

Christian Tetzlaff, Max-Planck-Institut für Dynamik und Selbstorganisation, Germany
Michael Fauth, University of Edinburgh, United Kingdom

Copyright © 2018 Hu, Mergenthal, Bingham, Song, Bouteiller and Berger. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Eric Hu, ehu@usc.edu
Jean-Marie Bouteiller, jbouteil@usc.edu

ORIGINAL RESEARCH article

A Glutamatergic Spine Model to Enable Multi-Scale Modeling of Nonlinear Calcium Dynamics

Introduction

Materials and Methods

Mechanistic Model

Spine Volume and Diffusion

Calcium Influx

Calcium Efflux

Buffers and Intracellular Calcium Stores

Inputs Into the Mechanistic Spine Calcium Model

Input-Output Model

Platforms and Computational Tools

Results

Calcium Dynamics Calibration and Validation With Published Experimental Data

Calcium Fluctuations Vary According to the Inter-spike Intervals Between Presynaptic and Postsynaptic Activity of the Neuron

Postsynaptic Calcium Activity Depends on the Distance Between the Synapse and the Soma

A Third Order, Multi-Input Input-Output Calcium Model Closely Replicates the Response of the Mechanistic Calcium Model at Lower Frequency

Discussion

Availability

Author Contributions

Conflict of Interest Statement

Acknowledgments

Supplementary Material

References

This article is part of the Research Topic

People also looked at