Hence, by changing the buffer concentration, a synapse can go from a facilitating state to a depressing state. To investigate the role of the readily-releasable-pool organization in the release probability, a first step is to use a vesicular replacement rate at the AZ by considering the time τDock for a vesicle from the readily releasable pool of vesicles to replace a vesicle that has just fused: in the extreme case τDock = 0, which corresponds to an immediate refilling of vesicles, and τDock = ∞ which corresponds to no refilling of vesicles: **Figure 4C** shows the behavior of a vesicular release with no refilling of vesicles, and a slow unbinding rate. Simulations show an increase in release probability after the second spikes, due to the calcium accumulation at the binding sites. After the third spike, the release probability is decreased, due to the lack of vesicles docked at the AZ, which would be interpreted as short-term depression. To conclude, a low unbinding rate is responsible for calcium accumulation in the sensor binding site that increases the release probability and defines short-term facilitation.

#### 6. CONCLUSION AND PERSPECTIVE

The lesson from modeling and numerical simulations of diffusion in the first hundred nanometers between docked vesicles and the plasma membrane is that this boundary layer is crucial for computing the vesicular release probability due to the critical position of vesicular sensors. This space is difficult to access experimentally and its role has been underestimated in shortterm plasticity. However, experimental approaches using calcium chelator such as BAPTA or EGTA confirm the role of calcium ions traveling from far away compared to the ones entering directly through VGCC located underneath a vesicle to trigger release.

Another key feature relevant for short-term plasticity is the structural organization and the spatial correlation between the distribution of VGCCs and vesicles. Do vesicles contain the same amount of close VGCCs? What defines the exact location of VGCC underneath vesicles? Can this number fluctuate? What happens after vesicular fusion? How are VGCC redistributed? It is possible that VGCCs are constantly moving to find local optimal sites location (Schneider et al., 2015).

Short-term facilitation is classically thought as the accumulation of calcium in the synaptic bulk (calcium hypothesis

of Katz and Miledi, 1968), due to various possibilities such as slow and fast buffers. But computational evidences (Guerrier and Holcman, 2017) reviewed here suggest that calcium accumulation at the sensor binding sites and not in the bulk is actually the determinant effect to pre-activate vesicular release (by binding a certain fraction of the sensor sites). A similar conclusion was reached in Bornschein et al. (2013) using fluctuation analysis, calcium imaging and numerical simulation analysis indicating that the residual calcium bound to the release sensors (see also Schneggenburger and Neher, 2000), after the first AP could cause Paired Pulse Facilitation at Purkinje neuron synapses.

A byproduct of facilitation should be asynchronous release, because the random calcium accumulation at sensor increases the time window when a vesicular can be released, due to ions arriving at random time from the bulk, thus leading to a high variability in the calcium arrival and the vesicular release times. However, synaptic facilitation requires a low concentration of calcium buffer, suggesting that for facilitating synapses, calcium buffers should be maintained at a low level. In general, the relation between facilitation and asynchronous release triggered by residual calcium can be mediated by the specific molecular composition of the vesicle: indeed synaptotagmin 7 (Syt7) (Turecek and Regehr, 2018) present at cerebellar granule cell synapses onto stellate cells and Purkinje cells (mice) plays such a role: in Syt7 KO, facilitation, and asynchronous release are smaller and shorter lived than in WT, although the residual calcium was unchanged. This is in contrast with synchronous released mediated by synaptotagmin I. It would be certainly interesting to model the distribution of various synaptotagmin on vesicles to determine when a vesicle will be more likely to be released asynchronously vs synchronously. In particular, is the distribution of these two SytI and 7 different on vesicles? Do they have very different calcium buffer affinity? These properties could be explored using numerical simulations.

In terms of calcium residual regulation, a low concentration level can be achieved by preventing ER or mitochondria to come in too close proximity of the AZ. In contrast, synchronous release is associated with a high calcium buffer concentration, preventing calcium ions to travel from far away. More specifically, twodimensional numerical simulations (Delvendahl et al., 2015), modeling essentially the synaptic bulk, revealed that fixed endogenous buffers with low affinity, characterized by a low calcium-binding ratio, mixed with mobile buffers with high affinity, results in a fast AZ calcium clearance. This results in synchronous high-frequency transmission (at 200 Hz). But it remains unclear how calcium fluctuations is maintained low (Modchang et al., 2010; Nadkarni et al., 2010; Dittrich et al., 2013; Weinberg, 2015) especially near the vesicular calcium sensor.

Finally, how the rate of vesicular release can vary over 6 orders of magnitude for the same synapse (Kochubey et al., 2011) also remains enigmatic. The cusp geometry and rare binding events may hold the key to the solution of this spectacular modulation of the vesicular release rate (Guerrier and Holcman, 2015).

#### REFERENCES

Modeling vesicular trafficking and recycling at various synapses including ribbon synapses (Thoreson et al., 2016) should clarify the organization of the pre-synaptic terminal (Graydon et al., 2011, 2014) and the effect of VGCC trafficking (Schneider et al., 2015).

## AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

## ACKNOWLEDGMENTS

We thank S. Hallermann for reviewing this manuscript and for pointing us his articles, that we now discuss here. CG would like to thank the Fyssen foundation for support.

**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.

Copyright © 2018 Guerrier and Holcman. 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.

# 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\* and 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.

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

#### 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

#### \*Correspondence:

Eric Hu<EMAIL_ADDRESS>Jean-Marie Bouteiller<EMAIL_ADDRESS>
Received: 28 December 2017 Accepted: 05 July 2018 Published: 27 July 2018

#### 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

# 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 multiscale 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.

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 multiscale 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.

### 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**.

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

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 | List of the parameters, specifications, and models used in the mechanistic model and their sources.

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<sup>3</sup> and the postsynaptic density volume set to 0.0032 µm<sup>3</sup> (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

Frontiers in Computational Neuroscience | www.frontiersin.org

Supplementary Table 6.

model described by Naoki et al. (2005). The diffusion rate is modeled as:

$$\left(\frac{d\left[Ca^{2+}\right]}{dt}\right)\_{\text{diffusion}} = D \frac{A}{d^{\*}V\_{spin}} (\left[Ca^{2+}\right]\_{\text{PSD}} - \left[Ca^{2+}\right])$$

(1)

$$\frac{d\left[\text{Ca}^{2+}\right]\_{\text{PSD}}}{dt} = -\frac{I\_{\text{N/DDA}}}{2FV\_{\text{PSD}}} - D \frac{A}{d^\*V\_{\text{PSD}}} (\left[\text{Ca}^{2+}\right]\_{\text{PSD}} - \left[\text{Ca}^{2+}\right]) \tag{2}$$

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, Vspine is the total volume of the spine, and - Ca2<sup>+</sup> PSD and - Ca2<sup>+</sup> Spine is the calcium concentration at the postsynaptic density and the spine, respectively. The postsynaptic density is dependent on INMDA, which is the calcium influx current from NMDA receptors. The equation for INMDA is described in Equation (3). The influx depends on Faraday's constant F and the volume of the postsynaptic density VPSD.

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:

$$I\_{\rm NMDA} = \mathfrak{g}\_{total \rm NMDA} \times V\_{\rm Ca} \tag{3}$$

$$\mathfrak{g}\_{total\text{NMDA}} = \mathfrak{g}\_{\text{NMDA}} \times \mathfrak{n}\_{\text{NMDA}} \tag{4}$$

INMDA indicates total calcium current that flows into the spine, gNMDA is the NMDA conductance to be described in the following paragraph and nNMDA represents the number of NMDA receptors with the given PSD volume and is set to 20 (Racca et al., 2000). VCa 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:

$$\mathbf{g}\_0 = \mathbf{g}\_1 + \frac{\mathbf{g}\_2 - \mathbf{g}\_1}{1 + e^{\alpha V}} \tag{5}$$

$$g\_{\text{max}} = \frac{\text{go}}{1 + \left(\frac{[Mg^{2+}]\_{\text{o}}}{K\_{\text{0}}}\right)e^{-\delta zFV/RT}}\tag{6}$$

$$\mathbf{g}\_{\text{NMDA}} = \mathbf{g}\_{\text{max}} \times \mathbf{O}(t) \tag{7}$$

Where g<sup>0</sup> represents the total conductance in the absence of any magnesium, g<sup>1</sup> and g<sup>2</sup> represent the open state conductances with one glutamate bound and two glutamate molecules bound, respectively. g1is set at 40 pS while g<sup>2</sup> is set at 247 pS. The value <sup>α</sup> <sup>=</sup> 0.01 represents the steepness of the transition between <sup>g</sup>1and g2. - Mg2<sup>+</sup> o represents the external magnesium concentration and is set at 1 mM. K<sup>0</sup> is the equilibrium constant for magnesium set at 3.57, F is Faraday's Constant (9.64867.104 C mol−<sup>1</sup> ), R is the molecular gas constant (8.31434 J mol−<sup>1</sup> K −1 ), 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:

$$V\_{Ca} = -\frac{RT}{2F} \log\left(\frac{[Ca^{2+}]\_o}{[Ca^{2+}]}\right) \tag{8}$$

VCa 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. [Ca2+]<sup>o</sup> is the extracellular calcium set as a constant concentration of 2 mM, while [Ca2+] 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 VCa 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 Rtypes 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:

$$I\_{\rm VDCC} = \lg\_{\rm VDC} \cdot m^2 \cdot h \cdot f\_{\rm drive} \tag{9}$$

Where IVDCC 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). fdrive 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. gVDCC 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:

$$I = I\_{MAX}\frac{\left[Ca^{2+}\right]^{h\_c}}{K^{h\_c} + \left[Ca^{2+}\right]^{h\_c}}\tag{10}$$

Here, IMAX is the maximum calcium current extruded from the model, K is the equilibrium constant and h<sup>c</sup> 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**.

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):

$$V\_{Na,Ca} = \frac{RT}{F} \text{(3ln( } \frac{[Na^{+}]\_o}{[Na^{+}]\_i} \text{)} - \ln\left(\frac{[Ca^{2+}]\_o}{[Ca^{2+}]\_i}\right) \tag{11}$$

R, T, and F represent the molecular gas constant, temperature, and Faraday's constant with the same values mentioned previously. IMAX 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 VCa from the Nernst equation, VNa,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 VNa,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<sup>3</sup> ) (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 backpropagating 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 (gtotalNMDA) 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:

$$\mu\_{\rm SISO}(t) = \ c\_0 + \sum\_{j=1}^{L} c\_1 \begin{pmatrix} j \end{pmatrix} \nu\_j(t) + \sum\_{j\_1=1}^{L} \sum\_{j\_2}^{j\_1} c\_{2s} \begin{pmatrix} j\_1, j\_2 \end{pmatrix} \nu\_{j\_1}(t) \ \nu\_{j2}(t)$$

$$+\sum\_{j\_1=1}^{L} \sum\_{j\_2=1}^{j\_1} \sum\_{j\_3=1}^{j\_2} c\_{3s} \begin{pmatrix} j\_1, j\_2, \ j\_3 \end{pmatrix} \boldsymbol{\nu}\_{j\_1}(t) \ \boldsymbol{\nu}\_{j2}(t) \boldsymbol{\nu}\_{j3}(t) \tag{12}$$

$$\nu\_{\dot{\jmath}}(t) = \sum\_{\tau=0}^{M} b\_{\dot{\jmath}}(\tau) \, \varkappa(t-\tau) \tag{13}$$

Where uSISO is the single input Volterra series up to 3rd order. L refers to the total Laguerre functions in the SISO model, and c0, c1, c2<sup>s</sup> , and c3<sup>s</sup> refer to the coefficients associated with the 0th, 1st, 2nd, and 3rd order response of the series, respectively. v<sup>j</sup> (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<sup>j</sup> (τ )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{aligned} u\_{\text{MISO}}(t) &= \left. c\_0 + u\_1(t) + u\_2(t) \right. \\ &\quad + \sum\_{k\_1=1}^{L\_1} \sum\_{k\_2=1}^{L\_2} c\_{2r} \left( k\_1, k\_2 \right) \nu\_{k1}^{\mu\_1}(t) \, \nu\_{k2}^{\mu\_2}(t) \\ &\quad + \sum\_{k\_1=1}^{L\_1} \sum\_{k\_2=1}^{L\_2} \sum\_{k\_3=1}^{L\_1} c\_{3r\_1} \left( k\_1, k\_2, k\_3 \right) \nu\_{k1}^{\mu\_1}(t) \, \nu\_{k2}^{\mu\_2}(t) \, \nu\_{k3}^{\mu\_1}(t) \\ &\quad + \sum\_{k\_1=1}^{L\_1} \sum\_{k\_2=1}^{L\_2} \sum\_{k\_4=1}^{L\_2} c\_{3r\_2} \left( k\_1, k\_2, k\_4 \right) \nu\_{k1}^{\mu\_1}(t) \, \nu\_{k2}^{\mu\_2}(t) \nu\_{k4}^{\mu\_2}(t) \end{aligned} \tag{14}$$

$$\nu\_j^{\mu\_i}(t) = \sum\_{\mathfrak{r}=0}^M b\_j(\mathfrak{r}) \, \mathfrak{x}^{\mu\_i}(t-\mathfrak{r}) \tag{15}$$

Here, uMISO is the equation for the multi-input Volterra series, and u<sup>1</sup> and u<sup>2</sup> are the single input Volterra series based on uSISO with the inputs the postsynaptic potential (V), and the NMDA receptor conductance (gtotalNMDA), respectively. The cross kernel components involve the number of Laguerre functions L<sup>1</sup> and L<sup>2</sup> from the first and second SISO components, respectively, and at all orders. v u k (t) is the associated basis functions convolved with inputs: postsynaptic potential (V), and the NMDA receptor conductance (gtotalNMDA). c<sup>r</sup> 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 pythonbased 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 <sup>∼</sup>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).

# 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 prepost 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.

# 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.

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 <sup>∼</sup>9µM, and without any resulting action potential (in blue, solid) that reaches <sup>∼</sup>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.

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.

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

FIGURE 6 | 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).

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.

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:

$$\frac{d\left[Ca^{2+}\right]}{dt} = I\_{\text{NMDA}}\left(t\right) - \left(\frac{1}{\tau\_{Ca}}\left[Ca^{2+}\right]\right) \tag{16}$$

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**).

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.

# 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 multiscale 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

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.

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.

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.

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 spiketiming 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 NMDAbased 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 threefold 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 <sup>×</sup> <sup>6</sup> <sup>×</sup> <sup>5</sup> <sup>µ</sup>m<sup>3</sup> 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

#### REFERENCES

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.

#### 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

calcium transients from individual proteins. Front. Synaptic Neurosci. 7:17. doi: 10.3389/fnsyn.2015.00017

located in dendritic spines. Neuron 53, 249–260. doi: 10.1016/j.neuron.2006. 12.017

Engineering in Medicine and Biology Society. Annual Conference (Boston, MA), 3581–3584.

**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.

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.

# Dopamine Receptor Activation Is Required for GABAergic Spike Timing-Dependent Plasticity in Response to Complex Spike Pairing in the Ventral Tegmental Area

#### Ludovic D. Langlois, Matthieu Dacher † and Fereshteh S. Nugent\*

Department of Pharmacology, F. Edward Hebert School of Medicine, Uniformed Services University of the Health Sciences, Bethesda, MD, United States

#### Edited by:

Vito Di Maio, Istituto di Scienze Applicate e Sistemi Intelligenti "Eduardo Caianiello" (ISASI), Italy

#### Reviewed by:

Fatuel Tecuapetla, Instituto de Fisiología Celular (IFC), Mexico Jose Bargas, Universidad Nacional Autónoma de México, Mexico

\*Correspondence:

Fereshteh S. Nugent<EMAIL_ADDRESS>
#### †Present address:

Matthieu Dacher, Université Pierre et Marie Curie, Physiologie de l'Insecte, Signalisation et Communication—UMR 1272, Versailles, France

> Received: 20 July 2018 Accepted: 30 August 2018 Published: 21 September 2018

#### Citation:

Langlois LD, Dacher M and Nugent FS (2018) Dopamine Receptor Activation Is Required for GABAergic Spike Timing-Dependent Plasticity in Response to Complex Spike Pairing in the Ventral Tegmental Area. Front. Synaptic Neurosci. 10:32. doi: 10.3389/fnsyn.2018.00032 One of the most influential synaptic learning rules explored in the past decades is activity dependent spike-timing-dependent plasticity (STDP). In STDP, synapses are either potentiated or depressed based on the order of pre- and postsynaptic neuronal activation within narrow, milliseconds-long, time intervals. STDP is subject to neuromodulation by dopamine (DA), a potent neurotransmitter that significantly impacts synaptic plasticity and reward-related behavioral learning. Previously, we demonstrated that GABAergic synapses onto ventral tegmental area (VTA) DA neurons are able to express STDP (Kodangattil et al., 2013), however it is still unclear whether DA modulates inhibitory STDP in the VTA. Here, we used whole-cell recordings in rat midbrain slices to investigate whether DA D1-like and/or D2-like receptor (D1R/D2R) activation is required for induction of STDP in response to a complex pattern of spiking. We found that VTA but not Substantia nigra pars compact (SNc) DA neurons exhibit long-term depression (LTDGABA) in response to a combination of positive (pre-post) and negative (post-pre) timing of spiking (a complex STDP protocol). Blockade of either D1Rs or D2Rs prevented the induction of LTDGABA while activation of D1Rs did not affect the plasticity in response to this complex STDP protocol in VTA DA neurons. Our data suggest that this DA-dependent GABAergic STDP is selectively expressed at GABAergic synapses onto VTA DA neurons which could be targeted by drugs of abuse to mediate drug-induced modulation of DA signaling within the VTA, as well as in VTA-projection areas, thereby affecting reward-related learning and drug-associated memories.

Keywords: ventral tegmental area, VTA, spike-timing dependent plasticity, STDP, synaptic plasticity, long-term depression, LTD, GABAergic synapses

**Abbreviations:** AKAP, A kinase anchoring protein; CaN, Calcineurin; DA, dopamine; D1R, dopamine D1-like receptor; D2R, dopamine D2-like receptor; IPSC, Inhibitory postsynaptic current; LTD, long-term depression; LTP, Long-term potentiation; LFS, Low frequency stimulation; NAc, Nucleus accumbens; bAPs, propagating action potentials; PKA, protein kinase A; RMTg, rostromedial tegmental area; STDP, Spike-timing dependent plasticity; SNc, Substantia pars compacta; VTA, Ventral tegmental area.

# INTRODUCTION

In the past few decades, synaptic plasticity has emerged as an important candidate mechanism for drug-induced changes in reward-related neural circuits where addictive drugs usurp synaptic mechanisms underpinning reward/motivational or aversive behavioral learning to alter dopamine (DA) signaling from the ventral tegmental area (VTA), a critical brain region involved in reward and motivation (Pignatelli and Bonci, 2015; Langlois and Nugent, 2017). Drug-induced synaptic plasticity has attracted considerable interest in studies of drug addiction because strong and durable memories associated with drug experience are demonstrated to promote compulsive drug taking, craving and relapse.

Experimental synaptic plasticity can be induced using traditional induction paradigms as well as spike-timingdependent plasticity (STDP) protocols. STDP is considered a physiological form of plasticity that relies on the relative timing of neuronal activity (Feldman, 2012). Specifically, STDP at glutamatergic and GABAergic synapses within DA-related neural circuits likely presents a critical synaptic learning rule for encoding natural reward-related learning and memory (Pawlak et al., 2010; Langlois and Nugent, 2017; Foncelle et al., 2018). The Hebbian learning rules of STDP dictates that long-term potentiation (LTP) is induced when presynaptic activity precedes postsynaptic spiking (pre-post spiking, positive timing), whereas reversing the order induces long-term depression (LTD, post-pre spiking, negative timing; Dan and Poo, 2006; Caporale and Dan, 2008). Recently, it has become clear that GABAergic synaptic plasticity in addition to the widely studied glutamatergic plasticity is targeted by drugs of abuse and drug-induced modulation of this plasticity could critically influence DA neuronal activity and DA release in VTA projection areas, as well as in local VTA microcircuits (Langlois and Nugent, 2017). We previously demonstrated that GABAergic synapses onto VTA DA neurons are able to exhibit a Hebbian heterosynaptic STDP where near-coincident and correlated activities of presynaptic glutamatergic neurons with postsynaptic DA neurons result in expression of LTP and LTD (that we call LTPGABA and LTDGABA). More importantly, we found that GABAergic synapses are predisposed to undergo LTDGABA in response to spike trains (a complex STDP protocol that includes bursts of a combination of pre-post and post-pre spiking; Kodangattil et al., 2013). Somatodendritically released DA within the VTA acts on DA D2-like receptors (D2Rs) and presents an important mechanism for controlling the excitability of DA neurons (Beckstead et al., 2004). There is also an emerging neuromodulatory role for DA in STDP (Pawlak et al., 2010). In fact, our previous study (Dacher and Nugent, 2011) demonstrated that the induction of LTD in response to a traditional LTD pairing protocol (a combination of low frequency stimulation (LFS) with modest depolarization) is dependent on D2R activation. Here we sought to explore the potential effects of pharmacological manipulation of DA transmission within the VTA on the induction of GABAergic STDP in response to complex STDP protocols. We found that DA action through either D1Rs or D2Rs is necessary and sufficient for the induction of LTDGABA in VTA DA neurons and that this DA-dependent plasticity is limited to GABAergic synapses in the VTA. Given that drugs of abuse increase DA neurotransmission within the VTA (Bradberry and Roth, 1989; Klitenick et al., 1992; Campbell et al., 1996; Rahman et al., 2003), our present findings provide an inhibitory synaptic mechanism by which drug-induced alteration of local VTA DA signaling could affect DA cell excitability and subsequently DA release in VTA DA circuits.

## MATERIALS AND METHODS

Brain slice preparation and electrophysiological recordings were conducted as described previously from 14 days to 21 days old Sprague-Dawley rats (Dacher et al., 2013; Kodangattil et al., 2013). Briefly, animals were anesthetized using isoflurane and quickly decapitated. The brain was rapidly dissected and placed into ice-cold artificial cerebrospinal fluid (ACSF) containing (in mM): 126 NaCl, 21.4 NaHCO3, 2.5 KCl, 1.2 NaH2PO4, 2.4 CaCl2, 1.00 MgSO4, 11.1 glucose, 0.4 ascorbic acid, saturated with 95% O2/5% CO2. Horizontal midbrain slices containing the substantia nigra pars compact (SNc) and VTA were cut (250 µm) and incubated in ACSF during at least 1 h at 34◦C. Slices were then transferred into a recording chamber in ascorbic acid-free ACSF at 28◦C. All experiments were carried out in accordance with the National Institutes of Health (NIH) guidelines for the care and use of laboratory animals and were approved by the Uniformed Services University Institutional Animal Care and Use Committee. All efforts were made to minimize animal suffering, and to reduce the number of animals used.

GABA<sup>A</sup> inhibitory post-synaptic currents (IPSCs) were recorded using a patch amplifier (Multiclamp 700B) under infrared-differential interference contrast microscopy. Data acquisition and analysis were performed using DigiData 1440A and pClamp 10 (Molecular Devices, Union City, CA, USA). In all experiments, 6,7-dinitroquinoxaline-2,3-dione (DNQX, 10 µM) and strychnine (1 µM) obtained from Sigma were added to block AMPA- and glycine-mediated synaptic currents, respectively to pharmacologically isolate GABA<sup>A</sup> IPSCs that were completely blocked by the GABA<sup>A</sup> receptor antagonist, bicuculline. Paired GABA<sup>A</sup> IPSCs were evoked using a bipolar stainless steel stimulating electrode placed 200–500 mm rostral to the recording site in the VTA at 0.1 Hz (duration 100 µs, 50 ms inter-stimulation interval) and recorded using KCl containing electrodes and whole-cell voltage-clamp in neurons held at −70 mV. Pipettes were filled with (in mM): 125 KCl, 2.8 NaCl, 2 MgCl2, 2 ATP-Na+, 0.3 GTP-Na+, 0.6 EGTA and 10 HEPES (pH adjusted to 7.28 with KOH, osmolarity adjusted to 275–280 mOsm with sucrose). Stimulation intensity was adjusted to evoke baseline synaptic responses ranged between −200 pA and −800 pA (approximately 50% of maximal responses). The cell input resistance and series resistance were monitored through the experiment and if these values changed by more than 10%, data were not included. The appearance of an Ih current (≥50 pA) in response to stepping cells from −50 mV to

depression (LTDGABA) in response to a complex spiking spike-timing-dependent plasticity (STDP) protocol. Panel (A) represents sample bursts of the complex spiking protocol for induction of LTDGABA. (B,C) Single and average experiments showing induction of STDP recorded in Ih(+) (presumably DA) neurons in VTA (filled square symbols) or SNc (filled circle symbols). At the arrow, STDP was induced. Insets: averaged inhibitory postsynaptic currents (IPSCs) before and 25 min after STDP protocol. In this and all figures, 10 consecutive traces from each condition were averaged for illustration as inset. Calibration: 100 pA, 25 ms (VTA: 75 ± 1.1% of pre-STDP values, F(3,22) = 7.5, p < 0.0001, n = 8; SNc: 103 ± 2.4% of pre-STDP values, F(5,70) = 1.26, n = 14). Values shown throughout figure are the mean ± SEM.

−100 mV was used to identify putative SNc/VTA DA neurons. As a standard protocol in our lab including the present study we consistently record from a region of the VTA (in the dorsal and caudal VTA) that is shown to contain mostly nucleus accumbens (NAc)- projecting DA neurons with Ih positivity (Margolis et al., 2006a,b; Zhang et al., 2010). In addition we consistently consider other electrophysiological criterions for identification of DA neurons (AP characteristics and frequency) that are also commonly used to identify putative DA neurons (Johnson and North, 1992; Dacher et al., 2013; Kodangattil et al., 2013).

To induce STDP, we first obtained a stable baseline for 10 min and then DA cells were taken to current clamp and received trains of a sub-threshold presynaptic stimulation paired with back propagating action potentials (bAPs/postsynaptic spiking) at 5 Hz. To evoke bAPs, cells were injected with direct somatic currents of 1.5 nA for 5 ms through patch pipettes. STDP protocols consisted of 30 trains of five bursts repeated at 0.1 Hz. To induce LTDGABA using post-pre pairing, each burst was composed of three bAPs at 50 Hz followed by a single presynaptic stimulation (negative timing, −5 ms). A complex STDP protocol was used to induce LTDGABA where each burst composed of three bAPs preceded with three presynaptic stimulations at 50 Hz (the complex pairing protocol included both positive timing, +5 ms and negative timing, −15 ms, **Figure 1A**). Values are presented as means ± SEM. Statistical significance was assessed using repeated measures ANOVA with significance level of p < 0.05. Levels of STDP are reported as averaged IPSC amplitudes for 5 min just before STDP induction compared with averaged IPSC amplitudes during the 5 min period from 25 min to 30 min after protocol. Interleaved control experiments were performed with experiments in which drugs were bath applied. Salts and all drugs were purchased from Sigma-Research Biochemicals International or Tocris Bioscience.

# RESULTS

# VTA but Not SNc DA Neurons Exhibit LTDGABA in Response to a Complex Pairing STDP Protocol

Previously, we showed that the synaptic efficacy of GABAergic synapses onto VTA DA neurons can be bi-directionally modified by pre/post spike pairing in a narrow time window (Kodangattil et al., 2013; Authement et al., 2015). Here, we used a more natural and complex pattern of spiking (**Figure 1A**) to induce STDP at GABAergic synapses onto midbrain VTA/SNc DA neurons. Consistent with our previous results (Kodangattil et al., 2013), we were able to induce LTDGABA in VTA DA neurons in response to a combination of both positive (+5 ms) and negative (−15 ms) timing (**Figure 1B**). On the other hand, we found that SNc DA neurons did not exhibit any form of plasticity in response to the same complex STDP protocol (**Figure 1C**). Given this finding, we only examined the effects of D1R/D2R drugs on the induction of STDP (LTDGABA) in VTA DA neurons.

# Induction of LTDGABA in the VTA by STDP Protocols Requires D2R Activation

LTDGABA at GABAergic synapses onto VTA DA neurons can also be triggered in response to a traditional LTD pairing paradigm using LFS paired with modest depolarization (Dacher and Nugent, 2011). LTDGABA triggered in response to both LFS-pairing paradigm and STDP protocols is dependent on the postsynaptic scaffolding A-kinase anchoring protein 79/150 (AKAP79/150) signaling complex which selectively controls GABAergic synaptic strength and mediates the opposing effects of protein kinase A (PKA) and calcineurin (CaN) on GABA<sup>A</sup> receptor trafficking in VTA DA neurons (Dacher and Nugent, 2011; Dacher et al., 2013; Authement et al., 2015). Since LTDGABA in response to the LFS-pairing LTD paradigm is also D2R-dependent and modulated by morphine (Dacher and Nugent, 2011), we tested whether LTDGABA induced by STDP protocols also requires D2R activation. We attempted to induce LTDGABA by a complex STDP protocol or a post-pre STDP protocol as previously described (Kodangattil et al., 2013; Authement et al., 2015) while a D2R antagonist, sulpiride (10 µM), was present in the perfusate throughout the experiment. Sulpiride was able to completely block the induction of LTDGABA in response to the complex STDP protocol (**Figure 2A**) as well as post-pre STDP protocol (**Figure 2B**) suggesting that both the pairing and STDP protocols trigger this D2R-dependent LTDGABA. It should be mentioned that D2R activation by a D2R agonist results in a rundown in GABAergic IPSCs (a chemical form of LTD) in VTA DA neurons that is dependent on CaN activity upon inhibition of PKA-AKAP150 anchoring (Dacher et al., 2013).

# Induction of LTDGABA by a Complex STDP Protocol Also Requires D1R Activation

The roles of D1Rs and D2Rs have been implicated in induction and modulation of STDP (Pawlak and Kerr, 2008; Shen et al., 2008; Pawlak et al., 2010; Xu and Yao, 2010; Ruan

FIGURE 2 | DA D1-like receptor (D1R) and D2R activation is required for induction of LTDGABA in response to the complex STDP protocol. Panels (A–D) show average experiments of STDP with sample traces from Ih(+) neurons in response to the complex or post-pre STDP protocols in drug-free artificial cerebrospinal fluid (ACSF; controls, filled square symbols) or drug bath application (open square symbols) experiments. Control LTDGABA group in response to the complex STDP protocol is similar to Figure 1B in (A,C,D) Representing the interleaved control experiments conducted with drug-treated slice experiments. (A) The D2R antagonist blocked the induction of LTDGABA in response to the complex protocol (sulpiride: 104 ± 1.8% of pre-STDP values, F(4,14) = 0.33, n = 7). (B) Sulpiride prevented the induction of LTDGABA in response to the pre-post protocol while control slices showed robust LTDGABA (control: 69 ± 1.8% of pre-STDP values, F(4,25) = 7.01, p < 0.0001, n = 7; sulpiride: 96 ± 1.5% of pre-STDP values, F(4,32) = 0.63, n = 8). (C) The D1R agonist did not affect the induction of LTDGABA in response to the complex protocol (SKF81297: 75 ± 1.2% of pre-STDP values, F(3,12) = 3.7, p < 0.05, n = 7). (D) The D1R antagonist blocked the induction of LTDGABA in response to the complex protocol (SCH23390: 89 ± 1.04% of pre-STDP values, F(4,32) = 2.01, n = 10). Calibration: 100 pA, 25 ms.

et al., 2014). Although DA via both D1R and D2R may exert opposing effects on STDP where D1Rs enable the induction of LTP and D2Rs mainly favor LTD, the cooperative actions of both DA receptors in induction of LTD of STDP have been shown (Xu and Yao, 2010). Here we tested our hypothesis that D1R activation blocks LTDGABA and reverses the direction of plasticity to promote LTPGABA. To test this, we attempted to induce LTDGABA in response to the complex pairing STDP protocol in the presence of a D1R agonist (SKF81297, 10 µM) or a D1R antagonist (SCH233390, 10 µM). While LTDGABA was unaffected in the presence of the D1R agonist (**Figure 2C**), the D1R antagonist blocked the induction of LTDGABA (**Figure 2D**) suggesting that D1R activation also facilitates the induction of LTD at GABAergic synapses onto VTA DA neurons.

### DISCUSSION

Here, we have extended our previous studies on GABAergic STDP in VTA DA neurons to investigate how local VTA DA signaling through somatodendritic release of DA during STDP within the VTA would affect STDP. We found that while GABAergic synapses onto VTA DA neurons were able to express LTDGABA in response to a complex STDP protocol, GABAergic synapses onto SNc DA neurons did not show any form of plasticity in response to the same protocol. This suggests that this form of inhibitory plasticity may be selectively expressed at GABA<sup>A</sup> synapses in the VTA. LTDGABA triggered by our complex STDP protocol or a post-pre STDP protocol in VTA DA neurons was dependent on D2R activation similar to LTDGABA triggered in response to a typical LFS-pairing protocol shown by our group (Dacher and Nugent, 2011). D1R activation seems to commonly facilitate the induction of LTP in response to STDP protocols although it could also contribute to induction of LTD of STDP at the synapse (Pawlak and Kerr, 2008; Zhang et al., 2009; Ruan et al., 2014; Brzosko et al., 2015). Thus, we further examined whether boosting endogenous DA action on D1R activity during STDP induction through bath application of a D1R agonist could reverse the direction of plasticity towards LTP. Not only were we unable to trigger LTP in response to the complex protocol in the presence of D1R agonist, the induction of LTDGABA was completely unaffected ruling out the possibility of masking LTD by a simultaneous induction of a D1R-dependent LTP at these synapses. Interestingly, we observed a blockade of LTDGABA in response to the complex protocol in the presence of D1R antagonist suggesting that local endogenous DA could engage D1Rs to promote this plasticity. The selectivity of induction of GABAergic STDP at GABAergic synapses in the VTA vs. the SNc may not be surprising considering the distinct anatomical and functional populations of VTA and SNc DA neurons (Beier et al., 2015; Shin et al., 2017). Given that VTA DA neurons are found to be more heterogeneous than SNc DA neurons (Margolis et al., 2006b) and the postsynaptic nature of DA-dependent LTDGABA triggered at GABAergic synapses in the VTA (Dacher and Nugent, 2011; Dacher et al., 2013; Kodangattil et al., 2013; Authement et al., 2015), our data suggest that the expression of this STDP as a uniform property of GABA<sup>A</sup> synapses in the VTA may be due to distinct intrinsic characteristics of the VTA microcircuits that differ from the SNc local circuitry. Consistently, it has been shown that extracellular DA levels released within the VTA and SNc differ as the dendritic release of DA within the VTA is far greater than DA release within the SNc (Rice et al., 1997; Ford et al., 2010). Therefore, we propose that this regional difference in release properties of DA at dendritic locations may underlie the selective expression of STDP in the VTA. Given that activation of either D1Rs or D2Rs by endogenous DA was sufficient to trigger this plasticity and considering different localization of DA receptors in VTA neuronal populations and presynaptic terminals innervating VTA neurons, it remains to be known how and where the activation of either presynaptic or postsynaptic DA receptors mediate this plasticity in the VTA. GABA<sup>A</sup> synapses onto VTA DA neurons mainly originate from VTA GABAergic neurons (that comprise 30% of VTA neuronal populations) and rostromedial tegmental area (RMTg) neurons (Barrot et al., 2012). D2Rs are mainly expressed postsynaptically on VTA DA neurons to provide an auto feedback inhibition of DA neurons (Beckstead et al., 2007). D2Rs are assumed to be expressed on GABAergic terminals where their activation facilitates the induction of a presynaptic endocannabinoidmediated LTD at GABAergic synapses onto VTA DA neurons (Pan et al., 2008). We have also shown that inhibition of PKA activity (which is the main downstream signaling mechanism for both D1R and D2R) or disruption of AKAP150-PKA association promotes LTDGABA in response to STDP protocols by favoring CaN activity and endocytosis of GABA<sup>A</sup> receptors in VTA DA neurons (Authement et al., 2015). Given that D2Rs can inhibit PKA activity and the postsynaptic locus of LTDGABA expression, we assume that D2R activation acts through this signaling pathway to promote LTDGABA in response to complex STDP protocols.  DA increases presynaptic release of GABA in the midbrain through D1R activation (Cameron and Williams, 1993). It has been shown that NAc D1R-expressing medium spiny neuronal projections inhibiting VTA DA neurons preferentially make GABA<sup>B</sup> synapses onto VTA DA neurons while these D1R expressing neurons of NAc inhibit VTA GABAergic interneurons via activating GABA<sup>A</sup> synapses (Barrot et al., 2012). It will be interesting to test whether the postsynaptic action of DA on D1Rs in VTA DA neurons promotes LTD through modulation of NMDA receptor (Schilström et al., 2006) or GABA<sup>B</sup> receptor activity originating from D1 expressing NAc (Cameron and Williams, 1993; Kamikubo et al., 2007) in VTA DA neurons. It should be noted that there are several limitations to our study including the young age of the rats, the suboptimal cooler temperature for recordings at 28◦C, the use of electrical stimulation rather than optogenetic stimulation of specific afferent inputs to the VTA, and also the complexity of mimicking a natural pattern of neuronal firings in an in vitro STDP induction protocol. In fact, a recent study has shown the heterogeneity of distinct GABAergic inputs to VTA DA neurons where only GABAergic inputs arising from the VTA GABAergic interneurons show short-term plasticity (Polter et al., 2018). In sum, we demonstrated that synaptic actions of DA within the VTA is required for induction of GABAergic STDP. The selective expression of this DA-dependent STDP in the VTA presents an important synaptic learning mechanism that can be targeted by drugs of abuse or stress to alter DA signaling within VTA DA circuits and significantly impact reward-related behavioral learning.

#### DATA AVAILABILITY

All relevant data is contained within the manuscript: all datasets (GENERATED/ANALYZED) for this study are included in the manuscript.

#### AUTHOR CONTRIBUTIONS

FN designed the experiments. LL, MD and FN performed electrophysiology experiments, analyzed the data and prepared the figures. LL and FN wrote the manuscript.

#### REFERENCES

#### FUNDING

This work was supported by the National Institutes of Health (NIH)—National Institute on Drug Abuse (NIDA) Grant#R01 DA039533 to FN. The funding agency did not contribute to writing this article or decide in submitting the article.

#### ACKNOWLEDGMENTS

The opinions and assertions contained herein are the private opinions of the authors and are not to be construed as official or reflecting the views of the Uniformed Services University of the Health Sciences or the Department of Defense or the Government of the United States. We are grateful to Nugent Laboratory members, Ryan Shepard and Robert Hammack for their constructive comments on the manuscript.

**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.

Copyright © 2018 Langlois, Dacher and Nugent. 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.

# Data Driven Models of Short-Term Synaptic Plasticity

Elham Bayat Mokhtari <sup>1</sup> , J. Josh Lawrence<sup>2</sup> and Emily F. Stone<sup>1</sup> \*

<sup>1</sup> Department of Mathematical Sciences, The University of Montana, Missoula, MT, United States, <sup>2</sup> Pharmacology and Neuroscience, Texas Tech University Health Sciences Center, Lubbock, TX, United States

Simple models of short term synaptic plasticity that incorporate facilitation and/or depression have been created in abundance for different synapse types and circumstances. The analysis of these models has included computing mutual information between a stochastic input spike train and some sort of representation of the postsynaptic response. While this approach has proven useful in many contexts, for the purpose of determining the type of process underlying a stochastic output train, it ignores the ordering of the responses, leaving an important characterizing feature on the table. In this paper we use a broader class of information measures on output only, and specifically construct hidden Markov models (HMMs) (known as epsilon machines or causal state models) to differentiate between synapse type, and classify the complexity of the process. We find that the machines allow us to differentiate between processes in a way not possible by considering distributions alone. We are also able to understand these differences in terms of the dynamics of the model used to create the output response, bringing the analysis full circle. Hence this technique provides a complimentary description of the synaptic filtering process, and potentially expands the interpretation of future experimental results.

#### Edited by:

Vito Di Maio, Istituto di Scienze Applicate e Sistemi Intelligenti "Eduardo Caianiello" (ISASI), Italy

#### Reviewed by:

Gianluigi Mongillo, Université Paris Descartes, France Andrey Olypher, Georgia Gwinnett College, United States Katsunori Kitano, Ritsumeikan University, Japan

#### \*Correspondence:

Emily F. Stone<EMAIL_ADDRESS>
Received: 01 February 2018 Accepted: 27 April 2018 Published: 22 May 2018

#### Citation:

Bayat Mokhtari E, Lawrence JJ and Stone EF (2018) Data Driven Models of Short-Term Synaptic Plasticity. Front. Comput. Neurosci. 12:32. doi: 10.3389/fncom.2018.00032 Keywords: short term plasticity, epsilon machines, synaptic filtering, mutual information, interneuron-pyramidal cell synapses, causal state splitting reconstruction

# 1. INTRODUCTION

Short term plasticity at the synapse level can have profound effects on functional connectivity of neurons. Through repetitive activation, the strength, or efficacy, of synaptic release of neurotransmitter can be decreased, through depletion, or increased, through facilitation. A single synapse type can display different properties at different frequencies of stimulation.

The role of synaptic plasticity and computation has been analyzed and reported on in numerous papers over the past 30 years. A review of feed-forward synaptic mechanisms and their implications can be found in Abbott and Regehr (2004). In this paper Abbott and Regher state "The potential computational power of synapses is large because their basic signal transmission properties can be affected by the history of presynaptic and postsynaptic firing in so many different ways." They also outline the basic function of a synapse as a signal filter as follows: Synapses with an initial low probability of release act as high pass filters through facilitation, while synapses with an initially high probability of release exhibit depression and subsequently serve as low pass filters. Intermediate cases in which the synapse can act as a band-pass filter, exist. Identifying synapse-specific molecular mechanisms is currently an active area of research, involving subtle changes in expression of myriad calcium sensor isoforms (synaptotagmins), subtly configured to alter the microscopic rates of synaptic release, facilitation, depression, and vesicle replenishment (Fioravante and Regehr, 2011; Chen and Jonas, 2017; Jackman and Regehr, 2017).

The underlying mechanisms creating these effects may be inferred by fitting an a priori model to synaptic response data. We parameterize such a model combining the properties of facilitation and depression (FD) at the presynaptic neuron with experimental data from dual whole-cell recordings from a presynaptic parvalbumin-positive (PV) basket cell (BC) connected to a postsynaptic CA1 (Cornu Ammonis 1 subregion) pyramidal cell, for fixed frequency spike trains into the presynaptic PV BC (Stone et al., 2014; Lawrence et al., 2015). We later examine the response of the model to an in vivo-like Poisson spike train of input, where the inter-spike interval (ISI) follows an exponential distribution, in Bayat et al. (submitted). Here we investigate the information processing properties of the synapse in question, following (Markram et al., 1998) and using standard calculations of entropy and mutual information between the input spike train and output response. This, however, left us unsatisfied, as it did not indicate the history dependence of the response, which we believe is one of the more interesting features of plasticity models that involve presynaptic calcium concentration. We attempted using multivariate mutual information measures, but this very quickly collapses due to the "curse of dimensionality" (Bellman, 1957). In this paper we try to resolve the question using methods of Computational Mechanics, creating unifilar HMMs called epsilon machines, that represent the stochastic process that our synapse model creates. As a bonus we are using data itself (albeit synthetic data) to create models of plasticity that can be used to classify properties of different types of synapses.

As stated in the abstract, methods from Information theory rely on distribution measures which inherently ignore the ordering of the measured data stream. We seek to incorporate this important feature of plasticity, the dependence of the response of the synapse on the prior sequence of stimulation, directly through the construction of causal state machines. This can only add to the understanding of the process in cases where the input stimulus train is known. In experiments where only the output postsynaptic response is known, this technique is particularly useful. While the machines themselves cannot be interpreted in a physiological way, the information they provide can be used to classify synaptic dynamics and inform the construction of physiological models. The point of the analysis is to gain as much accurate information from experiments in short term synaptic plasticity as possible without imposing the bias of an assumed underlying physical model. To create synthetic data we use a very simple but otherwise complete model of short-term plasticity that incorporates a "memory" effect through the inclusion of calcium build-up and decay. This has roots in a real physiological process (the flooding of calcium into the presynaptic terminal can trigger the release of neurotransmitter), but we are not interested per se in creating a biophysically complete model here. The calcium dynamics simply introduces another time scale into the model, one that is physiologically relevant. We wish to explore the effect of this additional time scale on the complexity of the process.

Computational Mechanics is an area of study pioneered by Crutchfield and colleagues in the 1990s, (Crutchfield and Young, 1989; Crutchfield, 1994; Shalizi and Shalizi, 2002). Finding structure in time series with these techniques has been applied in such diverse arenas as layered solids (Varn et al., 2002), Geomagnetism (Clarke et al., 2003), climate modeling (Palmer et al., 2002), financial time series (Park et al., 2007), and more recently ecological models (Boschetti, 2008) and large scale multi-agent simulations (Parikh et al., 2016). In neuroscience, to name a few only, we note one application to spike train data (Haslinger et al., 2013), and a recent publication by Marzen et al. on the time resolution dependence of information measures of spike train data (Marzen et al., 2015).

We employ some of the simplest ideas from this body of work, namely decomposing a discrete time-discrete state data stream into causal states, which are made up of sequences of varying length that all give the same probability of predicting the same next data point in the stream. The data are discrete time by construction, and made into discrete symbols through a partition, so the process can be described by symbolic dynamics. We use the Causal State Splitting Reconstruction (CSSR) algorithm on the data to create the causal states and assemble a HMM that represents the transitions between the states, while emitting a certain symbol. This allows us to classify the synapse types and gives an idea of the differences in the history dependence of the processes as well.

Using an a priori model for short-term synaptic dynamics and fitting it to data, while a perfectly valid approach, allows only for the discovery of the parameters in the model and possibly a necessary model reduction to remove any parameter dependencies (too many parameters in the model for the data set to fit). The alternative approach is to allow the data itself to create the model. From these "data driven" models, conclusions can be drawn about the properties of the synapse that are explicitly discoverable from the experimental data. The ultimate goal is a categorization of the types of processes a synapse can create, and an assignment of those to different synapse types under varying conditions. Note that the complexity or level of biophysical detail of our model synapse is not important to this end. In fact, the best way to calibrate this method is using the simplest possible model of the dynamics that captures the history dependence of the plasticity. This is not consonant with the goal of incorporating as many physiological features as possible, whether they affect the dynamics significantly or not. In fact, in most cases the limited data in any electro-physiological experiment precludes identifying more than a handful of parameters in an a priori model, a point we discuss in Stone et al. (2014). Our goal is to classify the sort of filter the synapse creates under certain physiological conditions, rather than to identify specific detailed cellular level mechanisms.

We are motivated in this task by the work of Kohus et al. (2016), in which they present a comprehensive data set describing connectivity and synaptic dynamics of different

**Abbreviations:** BC, basket cell; CA1, Cornu Ammonis, early name for hippocampus; FD, facilitation and depression; IPSC, inhibitory postsynaptic current; ISI, interspike interval; KL, Kozachenko and Leonenko; KSG, Kraskov, Stögbauer, and Grassberger; mAChR, muscarinic acetylcholine receptors; MCMC, Monte Carlo Markov Chain; MI, Mutual Information; NT, neurotransmitter; PSR, postsynaptic response; PV, parvalbumin-positive; HMM, Hidden Markov Model; CSSR, Causal State Splitting Reconstruction.

interneuron (IN) subtypes in CA3 using paired cell recordings. They apply dynamic stimulation protocols to characterize the short-term synaptic plasticity of each synaptic connection across a wide range of presynaptic action potential frequencies. They discovered that while PV+ (parvalbumin positive) cells are depressing, CCK+ (Cholecystokinin positive) INs display a range of synaptic responses (facilitation, depression, mixed) depending upon postsynaptic target and firing rate. Classifying such a wide range of activity succinctly is clearly useful in this context. The discovery that the rate of particular observed oscillations in these cells (called sharp wave ripples) may be paced by the shortterm synaptic dynamics of the PV+BC in CA3 demonstrates the importance of these dynamics in explaining complex network phenomena.

The paper is organized as follows. The construct for an experimental paper with section 2 and section 3 is not an immediately obvious partition of our work, but we use it as best we can. In the section 2 we describe the background on the synaptic plasticity model, and some analysis of its properties. We also cover the necessary background from Computational Mechanics. Finally we show how the techniques are explicitly applied to our data. In the section 3 we present the epsilon machines created from data from three types of synapses from our FD model: depressing, facilitating, and mixed, at varying input frequencies. Here we also indicate similarities and differences in the actual machines. In the section 4 we speculate on the reasons for these features by referring back to the original model. In the last section we indicate directions for future work.

### 2. MATERIALS AND METHODS

#### 2.1. Model of Synaptic Plasticity

In Stone et al. (2014), we parameterize a simple model of presynaptic plasticity from work by Lee et al. (2008) with experimental data from cholinergic neuromodulation of GABAergic transmission in the hippocampus. The model is based upon calcium dependent enhancement of probability of release and recovery of signalling resources (For a review of these mechanisms see Khanin et al., 2006). It is one of a long sequence of models developed from 1998 to the present, with notable contributions by Markram et al. (1998) and Dittman et al. (2000). The latter is a good exposition of the model as it pertains to various types of short term plasticity seen in the central nervous system, and the underlying dependence of the plasticity is based on physiologically relevant dynamics of calcium influx and decay within the presynaptic terminal. In our work, we use the Lee model to create a two dimensional discrete dynamical system in variables for calcium concentration in the presynaptic area and the fraction of sites that are ready to release neurotransmitter into the synaptic cleft.

In the rest of this section we outline the model, which is used to generate synthetic data for our study of causal state models, or epsilon machines, of short-term plasticity.

In the model the probability of release (Prel) is the fraction of a pool of synapses that will release a vesicle upon the arrival of an action potential at the terminal. Following the work of Lee et al. (2008), we postulate that Prel increases monotonically as function of calcium concentration in a sigmoidal fashion to asymptote at some Pmax. The kinetics of the synaptotagmin-1 receptors that binds the incoming calcium suggests a Hill equation with coefficient 4 for this function. The half-height concentration value, K, and Pmax are parameters determined from the data.

After releasing vesicles upon stimulation, some portion of the pool of synapses will not be able to release vesicles again if stimulated within some time interval, i.e., they are in a refractory state. This causes "depression;" a monotonic decay of the amplitude of the response upon repeated stimulation. The rate of recovery from the refractory state is thought to depend on the calcium concentration in the presynaptic terminal (Dittman and Regehr, 1998; Wang and Kaczmarek, 1998). The model has a simple monotonic dependence of rate of recovery on calcium concentration, a Hill equation with coefficient of 1, starting at some kmin, increasing to kmax asymptotically as the concentration increases, with a half height of K<sup>r</sup> .

The presynaptic calcium concentration itself, [Ca], is assumed to follow first order decay kinetics to a base concentration, [Ca]base. At this point we choose that [Ca]base = 0, since locally (near the synaptotagmin-1 receptors) the concentration of calcium will be quite low in the absence of an action potential. The evolution equation for [Ca] then is simply τca d[Ca] dt = −[Ca] where τca is the calcium decay time constant, measured in milliseconds. Upon pulse stimulation, presynaptic voltage-gated calcium channels open, and the concentration of calcium at the terminal increases rapidly by an amount δ (measured in µm): [Ca] → [Ca] + δ at the time of the pulse. Note that calcium build-up is possible over a train of pulses if the decay time is long enough relative to the ISI.