A unified computational model for cortical post-synaptic plasticity
- Simula Research Laboratory, Norway
- Department of Biosciences, University of Oslo, Norway
- Faculty of Science and Technology, Norwegian University of Life Sciences, Norway
- Department of Physics, University of Oslo, Norway
- The Krasnow Institute for Advanced Study, George Mason University, United States
Abstract
Signalling pathways leading to post-synaptic plasticity have been examined in many types of experimental studies, but a unified picture on how multiple biochemical pathways collectively shape neocortical plasticity is missing. We built a biochemically detailed model of post-synaptic plasticity describing CaMKII, PKA, and PKC pathways and their contribution to synaptic potentiation or depression. We developed a statistical AMPA-receptor-tetramer model, which permits the estimation of the AMPA-receptor-mediated maximal synaptic conductance based on numbers of GluR1s and GluR2s predicted by the biochemical signalling model. We show that our model reproduces neuromodulator-gated spike-timing-dependent plasticity as observed in the visual cortex and can be fit to data from many cortical areas, uncovering the biochemical contributions of the pathways pinpointed by the underlying experimental studies. Our model explains the dependence of different forms of plasticity on the availability of different proteins and can be used for the study of mental disorder-associated impairments of cortical plasticity.
Introduction
Synaptic plasticity in the neocortex has been under intense research since the first observations of neocortical long-term potentiation (LTP) (Komatsu et al., 1981; Lee, 1982). Although most often studied in brain slices, synaptic plasticity in the neocortex is a key phenomenon underlying vital mammalian brain processes ranging from formation and storage of memories to attentional selection (Roelfsema and Holtmaat, 2018). These processes are impaired in heritable mental illnesses such as schizophrenia and fragile X syndrome, as well as neurodegenerative diseases such as Alzheimer’s disease, all of which have been associated with deficits in cortical plasticity (Kantrowitz et al., 2017; Martin and Huntsman, 2012; Koch et al., 2014). Improved understanding of neocortical synaptic plasticity all the way from molecular to circuit level is therefore needed to further our understanding of these yet incurable diseases.
Similar to hippocampal synaptic plasticity (Larkman and Jack, 1995), synaptic plasticity in the neocortex is highly variable — the outcomes of any plasticity-inducing protocol depends on the cortical area, neuron type as well as details of the stimulation protocol (Castro-Alamancos et al., 1995; Froc and Racine, 2005; Sjöström et al., 2008; Feldman, 2009). Computational models provide a tool for efficient hypothesis testing of mechanisms of neocortical plasticity, which helps to overcome the challenges posed by excessive variability. The foundations of our mechanistic understanding of neocortical synaptic plasticity lie upon the phenomenological Bienenstock-Cooper-Munro (BCM) theory, which predicts that small synaptic activity (later attributed to small Ca2+ transients [Bear et al., 1987; Lisman, 1989]) cause long-term depression (LTD) whereas large synaptic activity (large Ca2+ transients) give rise to LTP (Bienenstock et al., 1982). Simple BCM-based models and the closely related models of spike-timing-dependent plasticity (STDP) have been widely used to explain the emergence of input-specific cell assemblies mediating, e.g., orientation selectivity (Shouval et al., 1997) or memory traces (Klampfl and Maass, 2013) in the cortex. These models, however, typically fail to provide a mechanistic understanding of the biochemistry within the synapse — namely, they do not reveal how various molecules downstream of Ca2+ regulate the induction and maintenance of plasticity occurring in neuronal circuits, their composite neurons and synapses of the cortex. Moreover, current models often ignore the joint contributions of neuromodulators, which are critical for inducing some forms of cortical synaptic plasticity (Meunier et al., 2017; Brzosko et al., 2019). These shortcomings impede testing biochemical mechanisms of heritable mental illnesses associated with impaired cortical plasticity.
In this work, we aim at filling this gap of knowledge by introducing a biochemically detailed, mass-action law-based model of neocortical post-synaptic plasticity that can be used to study the induction of plasticity in different genetic conditions and neuromodulatory states, and under various stimulation protocols. Despite the lack of biochemically detailed models of synaptic plasticity in the neocortex, models of intracellular signalling have been used to study LTP and LTD in the hippocampus (Bhalla and Iyengar, 1999; Jȩdrzejewska-Szmek et al., 2017), cerebellum (Gallimore et al., 2018), and striatum (Blackwell et al., 2019). These models permit systematic studies on how patterns of Ca2+ inputs to the post-synaptic spine, either alone or in combination with neuromodulatory actions, activate different signalling pathways leading to post-synaptic plasticity in the form of, e.g., AMPA-receptor (AMPAR) phosphorylation and membrane insertion. We integrate quantitative descriptions of the intracellular signalling pathways underlying synaptic plasticity in the neocortex into a unified model that is capable of describing both stimulation protocol-dependent plasticity, as well as neocortically observed neuromodulator-gated forms of STDP. We show that our model can be tuned by alterations of protein expression to reproduce not only BCM-like forms of plasticity but also experimental observations on neocortical plasticity from various cortical areas. Our results help to quantify and explain the differences in molecular constituents of different forms of neocortical LTP and LTD, and the different, data-fitted versions of our model can be directly used to examine the effects of chemical inhibitors and genetic manipulations of signalling proteins on synaptic plasticity in different cortical cells.
Results
Model construction
We reviewed the literature of molecular signalling pathways that needed for neocortical LTP/LTD, in particular in the post-synaptic spine of pyramidal cells (Table 1A). Three main pathways were highlighted in the experimental studies, namely, the protein kinase A (PKA), protein kinase C (PKC), and Ca2+/calmodulin-dependent kinase II (CaMKII) pathways. To construct a computational model of cortical post-synaptic plasticity that describes these pathways, we adopted mass-action law-based descriptions of these pathways from biochemically detailed models of post-synaptic LTP/LTD in other brain areas, namely, hippocampus, basal ganglia and cerebellum (Table 1B). We prioritised the model components from hippocampal models due to the relatively small ontological differences between hippocampus and neocortex (Kirsch and Chechik, 2016). We focused on the effects of these pathways on AMPARs due to the better description of intracellular regulation of AMPAR dynamics in comparison to that of NMDA and kainate receptors or voltage-gated ion channels. In short, we based our model on that of Jȩdrzejewska-Szmek et al., 2017, which describes the PKA- and CaMKII-dependent phosphorylation of AMPAR subunit 1 (GluR1), and added the metabotropic glutamate receptor (mGluR) and muscarinic acetylcholine M1 receptor-mediated activation of PKC from Kim et al., 2013 and Blackwell et al., 2019, respectively. Other types of receptors that interact with these pathways, such as serotonin (5HT) and dopamine receptors (He et al., 2015), have been shown to underlie certain types of neocortical plasticity. Dopamine D1/D5 receptors as well as serotonin 5HT4, 5HT6 and 5HT7 receptors are coupled to Gs proteins whereas 5HT2 receptors are Gq-coupled. The effects of these neurotransmitters would therefore be similar to those of norepinephrine and acetylcholine in our model (depending on the receptor composition in the post-synaptic neuron), and thus they are omitted in the present work. We then adopted the reactions describing PKC-dependent phosphorylation and endocytosis of AMPAR subunit 2 (GluR2) and reinsertion to the membrane from Gallimore et al., 2018, which allowed the representation of post-synaptic depression with our model. The pathways included in the model are illustrated in Figure 1. A description of the model calibration is given in Materials and methods, section 'Construction and calibration of the biochemically detailed model of post-synaptic plasticity in the cortex', and the full set of model reactions and initial concentrations is provided in Tables 3 and 4, respectively.
Table 1
Pathways contributing to cortical synaptic plasticity.
(A) Experimental evidence on the requirement of various molecular species for specific types of synaptic regulation in different cortical areas. (B) Model components needed for describing the modes of plasticity listed in (A). References are made to previous computational models describing these pathways. The types of phosphorylation of AMPAR subunit that mediate the plasticity are printed in bold.
| (A) | ||||
|---|---|---|---|---|
| Pathway components | Type of neurons | Type of regulation | Pre-/post-synaptic | References |
| CaMKII | Cingulate cortex | Esophageal acid-induced sensitisation | post-syn. | Banerjee et al., 2013 |
| CaMKII | Prefrontal cortex, pyramidal neurons | 5-HT1-induced modulation of AMPA currents | post-syn. | Cai et al., 2002 |
| -adr. receptors, PKA | Visual cortex, layer 4 pyramidal cells | Potentiation of AMPA currents | post-syn. | Seol et al., 2007 |
| M1 receptors, PKC | Visual cortex, layer 4 pyramidal cells | Depression of AMPA currents | post-syn. | Seol et al., 2007 |
| D1–PKA | Prefrontal cortex, pyramidal neurons | Potentiation of AMPA currents | post-syn. | Sun et al., 2005 |
| -adr. receptors | Frontal cortex | Potentiation of field EPSPs | n/a | Sáez-Briones et al., 2015 |
| PKC | Cultured cortical neurons | Internalisation of AMPARs | post-syn. | Chung et al., 2000 |
| ERK | Visual cortex | Potentiation of field EPSPs | n/a | Di Cristo et al., 2001 |
| (B) | ||||
| Molecular pathway | Cell type and references | |||
| Ca2+ → CaM → CaMKII | Hippocampal CA1 neuron Bhalla and Iyengar, 1999; Jȩdrzejewska-Szmek et al., 2017, generic Hayer and Bhalla, 2005, | |||
| cerebellar Purkinje cells Gallimore et al., 2018, striatal spiny projection neuron Blackwell et al., 2019 | ||||
| CaMKII → GluR1 S831p | Hippocampal CA1 neuron Jȩdrzejewska-Szmek et al., 2017 | |||
| -adrenergic receptors → cAMP | Hippocampal CA1 neuron Jȩdrzejewska-Szmek et al., 2017 | |||
| cAMP → PKA | Hippocampal CA1 neuron Bhalla and Iyengar, 1999; Jȩdrzejewska-Szmek et al., 2017, cerebellar Purkinje | |||
| cells Gallimore et al., 2018 | ||||
| PKA → GluR1 S845p | Hippocampal CA1 neuron Jȩdrzejewska-Szmek et al., 2017 | |||
| M1 receptors → PLC | Cerebellar Purkinje cells Gallimore et al., 2018 | |||
| PLC → PKC | Hippocampal CA1 neuron Bhalla and Iyengar, 1999, striatal spiny projection neuron | |||
| Kim et al., 2013; Blackwell et al., 2019 | ||||
| cerebellar Purkinje cells Kotaleski et al., 2002; Gallimore et al., 2018 | ||||
| PKC → GluR2 S880p | Cerebellar Purkinje cells Gallimore et al., 2018 | |||
Figure 1
Signalling pathways included in the model.
The PKA-pathway-related proteins and signalling molecules are highlighted by blue, PKC-pathway molecules by yellow, and CaMKII-pathway molecules by green colours. Reactions associated with a molecular species in parenthesis indicate a dependency on the denoted species — for details, see Table 3. Acronyms: β-AR – β-adrenergic receptor; AC1 and AC8 – adenylyl cyclase type 1 or 8; CaM – calmodulin; CaMKII – calmodulin-dependent protein kinase II; cAMP – cyclic adenosine monophosphate; DAG – diacylglycerol; Epac1 – exchange factor directly activated by cAMP 1; Gi, Gq and Gs – G-protein type I, Q, or S; GluR1 and GluR2 – AMPAR subunit 1 or 2; mGluR – metabotropic glutamate receptor; M1R – cholinergic receptor M1; NCX – Na+-Ca2+ exchanger; Ng – neurogranin; NMDAR – NMDA receptor; PDE1 and PDE4 – phosphodiesterase type 1 or 4; PIP2 – phosphatidylinositol 4;5-bisphosphate; PKA – protein kinase A; PKCt and PKCp – transiently or persistently active protein kinase C; PLC – phospholipase C; PMCA – plasma membrane Ca2+ ATPase; PP1 – protein phosphatase 1; PP2A – protein phosphatase 2A; PP2B – protein phosphatase 2B (calcineurin). In this work, the NMDARs are considered only in section 'Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses': in the rest of the work, Ca2+ is directly injected as a square-pulse current into the spine.
Ca2+ activates multiple pathways that regulate the post-synaptic plasticity in cortical PCs
All pathways of Table 1B are Ca2+-dependent, but due to the variability in binding rates and quantities of different Ca2+-binding molecules, some pathways become more easily activated than others. This permits LTP or LTD to be induced in a way that is sensitive to the amount of Ca2+ inputs and may serve as a basis for BCM-type rules of plasticity.
To examine the sensitivities of LTP- and LTD-inducing pathways to Ca2+, we simulated the injection of a prolonged square-pulse Ca2+ input of varying magnitude (illustrated in Figure 2A) into the post-synaptic spine and quantified the degree of activation of each of the Ca2+-binding molecules and the downstream signalling cascades. The simulations were carried out in the presence of mGluRs and -adrenergic and cholinergic neuromodulation, which were modelled as prolonged square-pulse inputs as well.
Figure 2
Ca2+ activates CaMKII, PKA, and PKC pathways.
(A) Illustration of the stimulus protocols with Ca2+ flux amplitudes 150 (green), 200 (cyan), and 250 (purple) particles/ms. (B–F) Time courses of Ca2+ (in nM) bound to buffers (B), pumps (C), PKC-pathway proteins (D), or CaM (E), and the concentration of free Ca2+ ions (F), according to NeuroRD (solid; averaged across eight samples) or NEURON (dashed) simulations. Colours indicate the Ca2+ flux used (see A). (B) Number of Ca2+ ions bound to Ca2+ buffers, that is immobile buffer and calbindin. (C) Number of Ca2+ ions bound to Ca2+ pumps and exchangers, that is PMCA and NCX. (D) Number of Ca2+ ions bound to PKC-pathway proteins PLC and PLA2. (E) Number of Ca2+ ions bound to CaM, in all its forms. (F) Cytosolic Ca2+ concentration (mM) (G) Degrees of activation of different Ca2+-binding proteins in a steady state (5 min after onset of Ca2+ input) as a function of the magnitude of Ca2+ flux. The x-axis shows the amplitude of the Ca2+ input (see panel A), and the y-axis shows the ratio of the underlying species in a Ca2+-bound form over the total number of the proteins. For CaM, only the CaM molecules bound by four Ca2+ ions are considered activated — in PLC, PLA2, and DGL, binding of only one Ca2+ ion is needed for activation. Here, the measured quantity of active PLC includes both Gq-bound and non-Gq-bound CaPLC. Inset: zoomed-in view on the red area. (H) Ratio of the steady-state concentration of PKA catalytic subunit over the theoretical maximum where all PKA molecules were dissociated into residuals and catalytic subunits. Colour of the curve indicates the amplitude of the β-adrenergic ligand flux (particles/ms). (I) Fraction of phosphorylated CaMKII subunits. (J) Fraction of (transiently or persistently) activated PKC. Colour of the curve indicates the amplitude of the cholinergic and glutamatergic ligand flux (particles/ms). The grey area in panels (G–J) represents Ca2+ inputs that cause cytosolic Ca2+ concentration to reach extremely high levels (>1 mM) that are likely to lead to apoptosis.
The injected Ca2+ quickly bound to Ca2+ buffers (immobile buffer and calbindin, Figure 2B) and pumps (PMCA and NCX, Figure 2C) as well as to the proteins of the PKC pathway (phospholipase A2 (PLA2) and C (PLC), Figure 2D): a 95% saturation was reached in 1–2 s (Figure 2B–D). In contrast, the activation of calmodulin (CaM) was slower (Figure 2E): a 95% saturation was reached in 32–53 s, depending on the magnitude of the Ca2+ input. Consistent with experimental literature, a vast majority of Ca2+ was quickly bound and only a small fraction remained free in the cytosol (Figure 2F).
To further illustrate the differences between the activation patterns of these pathways, we quantified the degrees of Ca2+ binding of these molecules in a steady state (5 min after the onset of Ca2+) and the overall activation/deactivation of downstream molecules as a function of the magnitude of the Ca2+ injection. Both PKC pathway-mediating proteins PLC, diacylglycerol lipase (DGL), and PLA2, along with the PKA and CaMKII pathway-related protein CaM became completely activated if large enough Ca2+ flux is given, but their degrees of activity varied across the magnitude of the injected Ca2+ flux (Figure 2G). DGL was most completely activated throughout the Ca2+ amplitude, owing to the large equilibrium constant of its Ca2+ binding. At extremely large Ca2+ fluxes, CaM was more completely bound by Ca2+ than PLC and PLA2 (Figure 2G), but at lower Ca2+ amplitudes, CaM remained largely unbound (Figure 2G inset). This is reflected in the activation patterns of the catalytic subunit of PKA (Figure 2H) and CaMKII (Figure 2I), both of which are dependent on the activation of CaM and thus had a small response at low Ca2+ amplitudes. PKC, by contrast, became activated at relatively small Ca2+ amplitudes (Figure 2J). Of these three pathways, the PKC pathway was dependent on the cholinergic ligands or the activation of the mGluRs (Figure 2J), and the PKA pathway was dependent on the availability of β-adrenergic ligands (Figure 2H). Taken together, these results highlight the need for large Ca2+ flux to the post-synaptic spine for the activation of the CaMKII pathway, relatively large Ca2+ flux for the activation of the PKA pathway, and relatively small Ca2+ flux for the activation of the PKC pathway.
High-frequency stimulation (HFS) causes LTP and low-frequency stimulation (LFS) causes LTD in GluR1-GluR2-balanced synapses
The Ca2+ flux entering the post-synaptic spine is extremely large during and after synaptic transmission and low otherwise, which causes the signalling pathways to be activated and deactivated in a more dynamic way than described in the previous section ('Ca2+ activates multiple pathways that regulate the post-synaptic plasticity in cortical PCs'). The activation of these pathways and their dependence on the stimulus protocol are difficult to study experimentally due to methodological constraints (e.g., side effects of fluorescence indicators, lack of signal calibration, and poor temporal or spatial resolution), but biochemically detailed models, such as the one considered in this work, can provide insights into the transient molecular mechanisms behind LTP and LTD. Our model is particularly well suited to study the mechanisms behind CaMKII-, PKA- and PKC-mediated phosphorylation of AMPAR subunits, which are important mediators of long-term plasticity (Wang et al., 2005). Phosphorylation of GluR1 subunits at S845 increases the insertion rate of the AMPAR into the membrane, thus leading to post-synaptic LTP (Diering et al., 2016). Conversely, phosphorylation of GluR2 subunits at S880 increases the rate of receptor endocytosis from the membrane, and has thus been observed to lead to post-synaptic LTD (Xia et al., 2000). However, it is not the number of the membrane-expressed AMPAR subunits alone that determine the strength of the synapse, but different compositions of the subunits have different single-channel conductances, and phosphorylation at S831 of the GluR1 subunit also affects the conductance of the channel (Oh and Derkach, 2005).
To simulate the reaction dynamics in the post-synaptic spine under realistic input patterns, we applied the 4xHFS and LFS protocols. Each input contained transient (3 ms) influxes of Ca2+ (1900 particles/ms) into the cytosol and glutamate (20 particles/ms), acetylcholine (20 particles/ms) and β-adrenergic ligand (10 particles/ms) into the extracellular subspace near the spine membrane. We used a balanced ratio (1:1) of GluR1 and GluR2 subunits. We recorded the time courses of the concentrations of all CaMKII-, PKA-, and PKC-pathway molecules contributing to LTP or LTD to monitor their activity during and following the stimulation protocols. We also recorded the numbers of membrane-inserted GluR1 and GluR2 and their state of phosphorylation and used Equation 5 for determining the maximal synaptic conductance.
In the 4xHFS protocol, which typically causes LTP in plasticity experiments, our model predicts a large increase in total synaptic conductance (Figure 3A) due to a radical increase in membrane-inserted GluR1 subunits (Figure 3B) and a decrease in GluR2 subunits (Figure 3C). These changes in membrane-expression of AMPAR subunits were dependent on activations of many signalling proteins in the CaMKII (Figure 3D–H), PKA (Figure 3I–M), and PKC (Figure 3N–R) pathways. First, the Ca2+ entry (Figure 3D) caused a rapid increase in half-activated calmodulin (bound by two Ca2+ ions; Figure 3E), leading to a longer-lasting increase in active calmodulin (Figure 3F). Calmodulin activation led to an increase in the concentration of phosphorylated CaMKII (Figure 3G), which phosphorylated the GluR1-type receptors at S831 (Figure 3H). The β-adrenergic input (Figure 3I), in turn, bound to the β-adrenergic receptors and activated the Gs proteins (Figure 3J), which bound to the adenylyl cyclase AC1 to produce cyclic adenosine monophosphate (cAMP, Figure 3K). cAMP bound to PKA to release the catalytic subunits of PKA (Figure 3L), which led to a phosphorylation of the GluR1-type receptors at S845 (Figure 3M) and thus to increased membrane expression of GluR1 subunits and total synaptic conductance (Figure 3A–B). Due to the simultaneous activation of the CaMKII pathway, a significant proportion of double phosphorylated GluR1-type receptors was observed (Figure 3H,M). As for the PLC–PKC pathway, glutamate (Figure 3N, blue) bound to mGluRs and acetylcholine (Figure 3N, green) bound to muscarinic receptors (M1), and the activation of these receptors contributed to the activation of Gq proteins (Figure 3O). The activated Gq proteins bound with PLC and metabolised phosphatidylinositol 4,5-bisphosphate (Pip2) into diacylglycerol (DAG, Figure 3P), which activated PKC (Figure 3Q). This led to the phosphorylation of GluR2-type receptors at S880 (Figure 3R), which caused the decrease in membrane expression of GluR2 observed in Figure 3C.
Figure 3 with 4 supplements see all
4xHFS activates CaMKII, PKA, and PKC pathways and leads to LTP (A–R), while LFS activates the PKC pathway and leads to LTD (S–U).
(A) Total synaptic conductance in response to 4xHFS, determined by the numbers of membrane-inserted GluR1s and GluR2s — see Equation 5. The stimulation starts at 40 s and lasts until 53 s. (B–C) Concentration of membrane-inserted GluR1s (B) and GluR2s (C) in response to 4xHFS. (D–H) Concentration of different species in the CaMKII pathway, namely, intracellular unbound Ca2+ (D), CaM bound with two Ca2+ ions (E), CaM bound with four Ca2+ ions (active CaM; F), phosphorylated CaMKII, bound or unbound by CaMCa4 (G), and S831-phosphorylated and double-phosphorylated GluR1 subunits (H) in response to 4xHFS. (I–M) Concentration of different species in the cAMP-PKA pathway, namely, β-adrenergic ligand in all its forms (I), activated (GTP-bound but not bound to ATP) Gs and Gi proteins (J), intracellular cAMP (K), catalytic subunit of PKA (L), and S845-phosphorylated and double-phosphorylated GluR1 subunits (M) in response to 4xHFS. (N–R) Concentration of different species in the PLC-PKC pathway, namely, glutamate and acetylcholine in all their forms (N), activated (GTP-bound but not bound to DAG) Gq proteins (O), intracellular DAG (P), activated PKC (Q), and S880-phosphorylated GluR2 subunits (R) in response to 4xHFS. S: Total synaptic conductance in response to LFS. (T–U) Concentration of membrane-inserted GluR1s (T) and GluR2s (U) in response to LFS, which starts at 40 s and lasts until 220 s. The solid lines represent stochastic (NeuroRD) simulation results, while the dashed lines represent data from deterministic (NEURON RxD) simulations. β-adrenergic ligands, glutamate, and acetylcholine are measured in numbers of particles as they reside both at the membrane (when bound to receptors) and at the extracellular subspace near the spine membrane (when unbound); other species measured in concentration.
The differences between NEURON and NeuroRD simulation results in Figure 3M were due to the stochasticity in NeuroRD simulator — both smaller and larger GluR1 phosphorylation levels compared to NEURON simulation results (Figure 3M, dashed) were obtained when NeuroRD simulations were run with different random number seeds (not shown).
In the LFS protocol, which typically causes LTD in the experiments, our model predicts a prominent (20%) decrease in total synaptic conductance (Figure 3S) due to a decrease in GluR2 subunits. In this protocol, the Ca2+ inputs are insufficiently large to activate CaM, and the Gs proteins remain deactivated as well (data not shown). In consequence, CaMKII and PKA pathways remain deactivated, and the effect of the LFS protocol on GluR1 phosphorylation and membrane insertion is small (Figure 3T). By contrast, the PKC pathway is almost as strongly activated as in the 4xHFS protocol (data not shown), leading to prominent S880 phosphorylation of GluR2 (data not shown) and removal of GluR2 from the membrane (Figure 3U).
The expression of both LFS-induced LTD and 4xHFS-induced LTP of these types is dependent on the presence of both GluR1 and GluR2 subunits: GluR1-deficient synapses failed to show 4xHFS-induced LTP (Figure 3—figure supplement 1A) and GluR2-deficient synapses failed to show LFS-induced LTD (Figure 3—figure supplement 1B). To show that our results were not an artefact of the tetramer formation rule (Equation 1–5), we applied an alternative tetramer formation rule where GluR1 and GluR2 subunits randomly dimerised and the dimers paired with like dimers (which disallows the emergence of heterotetramers with 1:3 or 3:1 proportion of GluR1:GluR2 subunits; Gan et al., 2015). We reproduced the LFS-induced LTD and 4xHFS-induced LTP using this dimer-of-like-dimers rule with a modified (35:65) balance of GluR1 vs. GluR2 subunits (Figure 3—figure supplement 2A).
In the above analyses, we used brief square-pulse fluxes of Ca2+to the synapse model, which is a simple representation of inputs during synaptic plasticity induction protocols. Alternatively, Ca2+ current entering the post-synaptic spines can be estimated by using multicompartmental biophysically detailed neuron models. We simulated a model of layer 2/3 pyramidal cell, stimulated with synaptic inputs from a 6xHFSt or LFS-1Hz protocol (see Materials and methods, section 'Modelling the Ca2+ inputs and neuromodulatory inputs'), to determine the Ca2+ inputs entering the post-synaptic spine through NMDA receptors (NMDARs). In accordance with Figure 3 and experimental data from somatosensory cortex (Heusler et al., 2000), our model predicted that 6xHFSt induced LTP whereas LFS-1Hz induced LTD (Figure 3—figure supplement 3). Here, the 6xHFSt protocol was used instead of 4xHFS to model the same protocol as in Heusler et al., 2000; our model would also predict an LTP for 4xHFS (data not shown). The HFS-induced LTP and LFS-induced LTD of Figure 3 could also be reproduced with alternative durations of neuromodulator inputs, including 10 min bath applications (Figure 3—figure supplement 4). These results indicate that our model can reproduce HFS-induced LTP and LFS-induced LTD also when using realistic NMDAR-conducted Ca2+ transients and that these forms of plasticity are robust to the temporal profile of the neuromodulatory inputs.
The activations of the above pathways are dependent on the magnitude and dynamics of the inputs to the model, namely, Ca2+, β-adrenergic and cholinergic ligands, and glutamate. All pathways leading to GluR1 and GluR2 phosphorylation and the consequent exocytosis and endocytosis are Ca2+-dependent: blocking Ca2+ entry completely abolished 4xHFS-induced LTP (Figure 4A) that followed GluR1 insertion (Figure 4B) and GluR2 endocytosis (Figure 4C). Blocking β-adrenergic ligands abolished the 4xHFS-induced LTP (Figure 4A) by suppressing the membrane-insertion of GluR1 (Figure 4B), but had no effect on GluR2 endocytosis (Figure 4C). Likewise, blocking β-adrenergic ligands had no effect on LFS-induced LTD (not shown). In contrast, LFS-induced LTD (Figure 4E) that followed GluR2 endocytosis (Figure 4G) was reduced by blockade of mGluR activation while the number of GluR1 subunits at the membrane remained unaffected (Figure 4F). This reduction was strengthened by simultaneous blockade of cholinergic inputs (Figure 4E–G, yellow traces). Counterintuitively, blocking mGluR and M1-receptor activation also reduced the amplitude of the 4xHFS-induced LTP (Figure 4A) by disabling GluR2 endocytosis (Figure 4C) while it had no effect on GluR1 insertion (Figure 4B). The reason for this is that in the PKC pathway-blocked case there is a smaller post-4xHFS membrane-bound GluR1 ratio (fraction of GluR1 subunits over all GluR subunits at the membrane) than in the control case, and thus the probability of AMPARs being homomeric GluR1 tetramers (which had a very large conductance compared to other tetramers; Equation 5) is much smaller in the former case than in control (Figure 4D). Although qualitatively similar difference can be observed in post-LFS membrane-bound GluR1 ratios between PKC pathway-blocked case and control, the probability of homomeric GluR1 tetramers and their contribution to the synaptic conductance are very small in both cases (Figure 4H) and thus the LFS-induced LTD is not affected.
Figure 4
4xHFS-induced LTP is dependent on β-adrenergic ligands and LFS-induced LTD is dependent on activation of mGluRs or cholinergic receptors.
(A–D) 4xHFS-induced LTP in the control case (dark purple), without Ca2+ inputs (blue), without β-adrenergic ligands (green), and under blockade of PKC pathway-activation (mGluRs or cholinergic receptors; yellow). (E–H) LFS-induced LTD in the control case (dark purple), under the blockade of mGluR activation (blue), and under blockade of both mGluRs or cholinergic receptors (yellow). (A, E) Total synaptic conductance. (B, F) Membrane expression of GluR1. (C, G) Membrane expression of GluR2. (D, H) The fraction of membrane-inserted GluR1 over all membrane-inserted GluR subunits (left), the probability of an AMPAR tetramer being homomeric GluR1 (middle), and the relative contribution of homomeric GluR1 subunits to the total conductance (i.e., summed conductance of homomeric GluR1 tetramers divided by the summed conductance of all tetramers; right). The bars represent the values at the end of the 4xHFS (D) or LFS (H) simulation with (dark purple) and without (yellow) PLC-activating ligands.
Taken together, our results show that cortical synapses expressing both GluR1 and GluR2 subunits can express a frequency-dependent form of post-synaptic plasticity (LTP for high-frequency inputs, LTD for low-frequency inputs) that is gated by neuromodulators affecting the PKA and PKC pathways. Our findings also lend support to that GluR2 endocytosis may lead to either potentiation (Figure 4A) or depression (Figure 4E), depending on the prevalence of the GluR1 subunits at the membrane.
Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses
Cortical synapses typically exhibit a type of synaptic plasticity, namely STDP, that is dependent on both the pre- and post-synaptic activity. According to a classical model, the differences in the outcome of STDP for different pairing intervals of pre- and post-synaptic stimulus are explained by different amount of Ca2+ entering the post-synaptic spine, which is affected by both the pre-synaptically released glutamate and the elevation of post-synaptic membrane potential. Biophysically detailed neuron modelling offers a powerful tool for determining the size of these Ca2+ inputs as a function of the pairing interval.
We considered the LTP/LTD response to paired stimulation protocol using a multicompartmental model of a layer 2/3 pyramidal cell (Figure 5A; Markram et al., 2015). We placed a synaptic spine with volume 0.5 μm3 at a random location on the apical dendrite, 250–300 μm from the soma (Figure 5A, thick, black branches), and stimulated the head of the spine with glutamatergic synaptic currents (Hay and Segev, 2015; Markram et al., 2015; Figure 5A, black traces, top). In parallel, we stimulated the soma with a burst of four short (2 ms) supra-threshold square-pulse currents (Figure 5A, black traces, bottom). Given that approximately 10% of the NMDAR-mediated currents and none of the AMPAR-mediated currents are conducted by Ca2+ flux, we could determine the number of Ca2+ ions entering the spine at each time instant following the onset of the synaptic input (Figure 5A, grey traces). This experiment was repeated using different inter-stimulus intervals (ISI) between the synaptic and somatic stimuli and averaged across trials. The membrane potential dynamics at the post-synaptic spine depended on the ISI (Figure 5B–D), largest effects response being obtained with near-coincident stimuli (Figure 5C). The higher the membrane potential at the spine, the higher the value of the variable describing the Mg2+-gate opening in the NMDA receptor (Figure 5B–D, insets) (Hay and Segev, 2015; Markram et al., 2015). Thus, the Ca2+ flux time course varied across the pairing ISIs (Figure 5—figure supplement 1). These Ca2+ flux time series were imported into our biochemical model (Ca2+ transients in the spine model showed in Figure 5E–G), which allowed us to predict the magnitude of GluR subunit phosphorylation and membrane insertion for each pairing interval. When added as bath application, the β-adrenergic and cholinergic ligands were simulated by prolonged injections of 50 particles/s for 10 min, starting 8 min before the STDP protocol — these neuromodulators alone (without the electric stimulation-mediated Ca2+ inputs) did not cause synaptic plasticity. Throughout the experiments, the activation of mGluRs was blocked.
Figure 5 with 1 supplement see all
Layer 2/3 pyramidal cell plasticity in response to STDP protocol depends on neuromodulatory state and pairing interval.
(A) Layer 2/3 pyramidal cell morphology (grey, thin), locations of synaptic input highlighted (black, thick). Inset: Illustration of the inputs (black) and the recorded synaptic intracellular Ca2+ (grey). Scale bar 200 μm. (B–D) Membrane potential at the dendritic spine when the pre-synaptic stimulation onset is 50 ms after (B), at the same time as (C), or 50 ms prior to (D) the onset of the last somatic stimulus. Inset (red): Mg2+-gate variable as a function of time, ranging from −80 ms to 140 ms in a similar manner as the data in the main panel. (E–G) Concentration of free Ca2+ in the dendritic spine according to the biochemical spine model when the pre-synaptic stimulation onset is 50 ms after (B), at the same time as (C), or 50 ms prior to (D) the onset of the last somatic stimulus. (H–J) No LTD was induced by the stimulation protocol (1 Hz paired with post-synaptic stimulation for 2 min) in the absence of M1-receptor activation, but pairing-interval-dependent LTP was induced in presence of β-adrenergic inputs. (K–M) Pairing-interval-dependent LTD was induced when the synaptic input was coupled with cholinergic inputs, and STDP was induced when both cholinergic and β-adrenergic inputs were present. (H, K) Relative concentration of GluR1 at the membrane 16 min after the stimulation onset (normalised by concentration of membrane-inserted GluR1 at rest). (I, L) Relative concentration of GluR2 at the membrane 16 min after the stimulation onset (normalised by concentration of membrane-inserted GluR2 at rest). (J, M) Relative synaptic conductance (Equation 5) 16 min after the stimulation onset (normalised by synaptic conductance at rest).
We first confirmed that the membrane expression of the glutamate receptors was not strongly affected by paired synaptic and somatic stimulation in the absence of β-adrenergic (which activates the PKA pathway) and cholinergic (which activates the PKC pathway) neuromodulation. Our model predicted that there is little change in the membrane expression of GluR1 and GluR2 type receptor subunits in this stimulation protocol (Figure 5H–I, purple). Consequently, our model reproduced the observation (Seol et al., 2007) that this stimulation protocol led to little change in predicted synaptic conductance (Figure 5J, purple).
We next considered the paired synaptic-somatic stimulation in the presence of β-adrenergic ligand. Our model predicted a prominent (up to 70%) increase in GluR1 membrane expression with little effect on GluR2 membrane expression (Figure 5H–I, blue). The predicted increase in GluR1 membrane expression (Figure 5H) and the consequent increase in synaptic conductance (Figure 5J, blue) were most prominent when the ISI was around 20–80 ms, and modest for large ISIs. These predictions are consistent with the experiments where an ISI-dependent potentiation of the EPSCs in the presence of β-adrenergic receptor agonists and absence of cholinergic agonists was observed (Seol et al., 2007).
When β-adrenergic neurotransmission was blocked but M1 receptors were activated by cholinergic ligands, the model predicted a prominent (up to 60%) decrease in GluR2 membrane expression, with little effect on GluR1 membrane expression (Figure 5K–L, purple). Our model of synaptic conductance (Equation 5) predicted a decrease in total conductance in a GluR1-GluR2-balanced synapse for this condition (Figure 5M, purple), which is in line with the experimental data (Seol et al., 2007). The depression takes place throughout the tested ISIs, but the effect was smallest for ISIs very close to zero due to the counteracting effects of GluR1 membrane-insertion (Figure 5K, purple). Finally, when both β-adrenergic and cholinergic neurotransmission were active, our model predicted an increased GluR1 membrane expression and decreased GluR2 membrane expression, both of which were ISI dependent (Figure 5K–L, blue). In these simulations, the predicted synaptic conductance was increased for small and moderate pre-post intervals and decreased otherwise (Figure 5M, blue), which is qualitatively similar to experimental data (Seol et al., 2007). These results are dependent on the availability of both GluR1 and GluR2 subunits at the post-synaptic spine: in simulations where GluR1 or GluR2 subunits were absent, only LTD (Figure 3—figure supplement 1C) or LTP (Figure 3—figure supplement 1D), respectively, was induced by the STDP protocol. In a similar manner as the HFS- and LFS-induced plasticity in Figure 3—figure supplement 2A, we could reproduce the STDP using the dimer-of-like-dimers tetramer formation rule with a GluR1 fraction of 35% (Figure 3—figure supplement 2B). Taken together, our model with balanced numbers of GluR1 and GluR2 subunits reproduces the neuromodulator-gated STDP observed in layer 2/3 pyramidal cells of the visual cortex.
The combination of our biochemically detailed model with the biophysically detailed model of layer 2/3 pyramidal cell model provides a compelling means of hypothesis testing for cortical STDP in this cell type. We analyzed how the shape of the STDP curve of Figure 5M is affected by the number of spikes in each post-synaptic burst stimulus. Our simulations suggest that decreasing the number of spikes per burst decreases the amplitude of both LTP and LTD in the STDP protocol and, in particular, brings the LTD for large post-pre ISIs close to zero (Figure 6A). These alterations are mediated by changes in both the level of membrane-insertion of GluR1 and endocytosis of GluR2 subunits (Figure 6A, insets). For small and moderate pre-post ISIs, the effects of decreasing the number of post-synaptic stimuli on the STDP curve are expected: both GluR1 insertion and GluR2 endocytosis are of smaller amplitude, and hence the dampened LTP amplitude (Figure 6A). By contrast, for post-pre ISIs and large pre-post ISIs, decreasing the number of post-synaptic stimuli results in larger amplitude of GluR1 insertion and GluR2 endocytosis, which yields a dampened LTD for post-pre ISIs and strengthened LTP for large pre-post ISIs (Figure 6A). These counter-intuitive effects can be explained by the accumulation of small-conductance K+ (SK) conductance in the post-synaptic neuron: the larger the number of post-synaptic pulses in the pairing burst, the larger the SK currents (Figure 6B). The SK current decays slowly (matching the Ca2+ concentration decay), and remnants of the SK currents can be observed as long as 200 ms after the post-synaptic stimulus (Figure 6B inset). For large post-pre ISIs, the number of spikes per burst has little effect on the Ca2+ transients during the post-synaptic stimulation (Figure 6C inset), but the SK currents activated by a large number of spikes per burst contribute to significantly decrease the Ca2+ transient caused by the pre-synaptic stimulus during the decay period (Figure 6C). By contrast, for small pre-post intervals, the additional spikes in the post-synaptic stimulus significantly contribute to the Ca2+ transients (Figure 6D). To show that the effects of the number of spikes per post-synaptic burst are mediated by the SK current, we ran the simulation of Figure 5J using a partial to complete blockage of the SK currents in the biophysically detailed simulations of the layer 2/3 pyramidal cell. The paired-pulse protocol of Figure 5M (involving both β-adrenergic and cholinergic neuromodulation) caused an STDP in all cases, but decreasing the SK conductance shortened the post-pre LTD window and decreased the amplitude of LTD (Figure 6E). Similar effects were obtained with a decrease of Ca2+-channel conductances (not shown), which is in agreement with the data of Nevian and Sakmann, 2006. Our model predictions also agree with the observation that the plasticity outcome is not determined by Ca2+ transient amplitude (Nevian and Sakmann, 2006), instead, our model suggests that the total Ca2+ is a better predictor of the plasticity outcome: the correlation coefficient between the post-STDP synaptic conductance and the peak Ca2+ transient amplitude (see Figure 5—figure supplement 1E) was 0.53, while that between the post-STDP synaptic conductance and the mean Ca2+ input during the inter-stimulus interval (see Figure 5—figure supplement 1F) was 0.96 (Figure 6—figure supplement 1).
Figure 6 with 1 supplement see all
The STDP curve of layer 2/3 pyramidal cells is affected by the number of post-synaptic stimulus pulses associated with the pre-synaptic input.
(A) The STDP curves of Figure 5M when the number of spikes per post-synaptic burst was 1 (yellow), 2 (green), 3 (blue), or 4 (as in Figure 5; dark purple). Inset: relative concentrations of membrane-inserted GluR1 (top) or GluR2 (bottom) subunits — see Figure 5K–L for reference. (B) Top: somatic membrane potential time course (aligned according to the onset of the first stimulus) for different numbers of post-synaptic stimulus pulses. Bottom: somatic SK current-density time course in the four conditions. Inset: the SK current densities 200 ms after the onset of the first post-synaptic stimulus. (C–D) Ca2+ flux to the dendritic spine when the pre-synaptic stimulation onset is 200 ms after (C) or 30 ms before (D) the onset of the last post-synaptic stimulus. (E) The STDP curves of Figure 5M when the number of spikes per post-synaptic burst was four but the somatic SK conductance parameter was either normal (dark purple), 50% smaller (magenta), or 80% smaller (pink).
The model predicts multimodal, protein concentration- and neuromodulation-dependent rules of plasticity
Cortical neurons express a variety of forms of LTP/LTD depending on the brain region and cell type. In computational studies, neocortical plasticity is most typically described by simple rules according to which small-amplitude Ca2+ inputs lead to depression of the synapse whereas large-amplitude inputs lead to potentiation. Apart from a few examples Castellani et al., 2001; d’D'Alcantara et al., 2003; Castellani et al., 2005; Honda et al., 2013, these models typically do not describe the intracellular signalling machinery leading to the resulting plasticity Holthoff et al., 2002; Karmarkar et al., 2002; Badoual et al., 2006; Cornelisse et al., 2007; Kubota and Kitajima, 2008; Urakubo et al., 2008. Unlike biochemically detailed models, the simple models cannot be used to explore whether and how the prevalence of different plasticity-related proteins gives rise to various types of LTP/LTD or their impairments, which is an important question in the study of mental disorders with deficits in cortical plasticity. Here, we analysed the biochemical underpinnings of different types of plasticity rules using our unified model of cortical plasticity in order to predict the conditions for different forms of plasticity.
In a similar fashion to section 'Ca2+ activates multiple pathways that regulate the post-synaptic plasticity in cortical PCs', we simulated our model of the post-synaptic spine when stimulated with a prolonged (5 min) square-pulse influx of Ca2+ and neuromodulators. We randomly altered the model parameters controlling the initial concentrations of different proteins, namely, the ratio of GluR1 to all GluR subunits (i.e., , from here on referred to as GluR1 ratio), the concentration of NCX (regulating the rate of Ca2+ decay from the spine), and the concentrations of PKA-pathway and PKC-pathway proteins (upstream of PKA and PKC). Alterations of the initial concentration of CaMKII (the only molecule in our model that exclusively affects the CaMKII pathway) had little effect in most domains of plasticity considered here (not shown), and thus, we omitted it in this analysis. We sampled these parameters from the following intervals: GluR1 ratio from the interval from 0 to 1 (keeping the total concentration of GluR subunits fixed at 540 nM), NCX concentration from the interval from 0 to twice the original value (2 × 0.54 mM), and the PKA and PKC-pathway factors and from the interval from 0 to 2 (see Materials and methods, section 'Parameter alterations and model fitting'). We simulated the post-synaptic spine 150,000 times using different random values for these parameters under zero, low (50 particles/ms), medium (150 particles/ms), and high (250 particles/ms) levels of Ca2+ input.
We classified the parameter sets based on the total synaptic conductance 15 min after the onset of the stimulation with the high Ca2+ flux (250 particles/ms): the relative synaptic conductance varied between 0.16 and 5.92, and thus, we grouped the parameter sets to 16 classes using a bin size of 0.36 (Figure 7A). We then analysed the parameter distributions and their co-variations within these classes and how the different parameters affected the shape of the LTP/LTD curve within each class. A special subset of the LTP/LTD curves were the BCM-type plasticity curves, where either 50 or 150 particles/ms Ca2+ injection resulted in LTD and the 250 particles/ms resulted in LTP.
Figure 7 with 1 supplement see all
The fraction of GluR1s, number of Ca2+ extrusion proteins, and the concentrations of PKA and PKC-pathway proteins in the post-synaptic spine determine the type of LTP/LTD in the post-synaptic spine.
(A) The LTP/LTD curves for all 16 classes. Four values of Ca2+ input amplitude were considered: 0, 50, 150, and 250 particles/ms (x-axis; repeated and overlaid for space). The y-axis shows the relative synaptic conductance, that is, total synaptic conductance 15 min after the onset of the Ca2+ input divided by the total synaptic conductance before the Ca2+ input. 20 representative parameter sets are displayed from each class, coloured from purple (lowest relative synaptic conductance response for medium Ca2+ input) to green (highest conductance). The black, dashed trace in class six represents the model with the default concentration parameters. (B–E) Distribution of model parameters, that is, GluR1 ratio (B), NCX-concentration coefficient (C), PKA pathway-concentration coefficient (D), and PKC pathway-concentration coefficient in the 16 classes. Class 6 (purple) highlighted for further analysis. F–H: GluR1 ratio plotted against NCX-concentration coefficient (F), (G), and (H) in class 6. The contours represent the distribution of parameters (N = 5837) that produced class-6 plasticity. No parameters yielding class-6 plasticity were found beyond the purple contour, and the inner contours cover the parameter space where the distribution is higher than 0%, 20%, 40%, 60% or 80% of the maximal density value. The black and red markers represent parameter sets that produced two plasticity subclasses, namely, one where the total deviance (summed absolute difference) from the BCM-type plasticity produced by the default parameter set (black, N = 145) or from a linearly increasing LTP (red, N = 183) was less than 0.2 (a.u.). Inset: The LTP/LTD plasticity curves of the two subclasses. The thick lines represent the centre of the subclasses (black: relative conductances in response to 50, 150, and 250 Ca2+ ions/ms: 0.76, 0.96, 2.24; red: relative conductances in response to 50, 150, and 250 Ca2+ ions/ms: 1.41, 1.83, 2.24).
Our model with the standard protein concentrations (GluR1 ratio 0.5, [NCX] = 0.54 mM, ) produced a BCM-type curve in class 6 (Figure 7A, black dashed curve). Classes 11–16 exhibited the strongest LTP, with large synaptic conductance for both 150 and 250 particles/ms Ca2+ injection, whereas classes 1 and 2 only exhibited LTD (Figure 7A). Classes 3–12 exhibited BCM-type of plasticity but the majority of the LTP/LTD curves were of non-BCM type in each class (Figure 7A).
Three parameters — the GluR1 ratio, NCX concentration and , differed significantly across the 16 classes (Figure 7B–D). Low GluR1 ratio was needed for strong LTD and medium or high GluR1 ratio for strong LTP (Figure 7B). However, the strongest forms of LTP (classes 11–16) were induced only when GluR1 ratio was smaller than 1 (Figure 7B), because a very low number of GluR2 subunits implied that the synapse has many homomeric GluR1 tetramers at a basal state, and thus stimulation-induced GluR1 exocytosis and GluR2 endocytosis did not radically increase the number of homomeric GluR1 tetramers (Equation 5, see also Figure 4D). For LTD and moderate LTP (classes 1–5), any NCX concentration and PKA-pathway coefficient could be used, but very strong LTP (classes 10–16) required a small to medium NCX concentration (Figure 7C) and a medium to large PKA-pathway coefficient (Figure 7D). By contrast, PKC-pathway coefficient alone was not predictive of plasticity outcome (Figure 7E).
The model results for the large parameter distributions of Figure 7B–E imply that there are manifestly different combinations of parameters that lead to the same LTP/LTD outcome. To analyse this intrinsic variability, we studied the distributions of the model parameters within the class of moderate LTP (class 6, 132–168% LTP for 250 particles/ms; indicated by purple boxes in Figure 7B–E) in more detail. Dependencies among the four parameters could be observed in 2-dimensional contour plots of the parameter distributions (Figure 7F–H). With large () GluR1 ratios, any NCX, PKA or PKC concentration could be used, but with smaller () GluR1 ratios, smaller NCX concentration (Figure 7F) or larger PKA-pathway coefficients (Figure 7G) were needed to obtain class-6 type of plasticity. To illustrate how these parameters affect the shape of the plasticity curves within class 6, we plotted the parameter sets that produced a BCM-type LTP/LTD curve similar to the one produced by our default model (Figure 7F inset, black) or an LTP curve that was linear within this regime (Figure 7F inset, red). Moderate GluR1 ratios (0.40–0.57; Figure 7F, black) and moderate NCX concentrations (0.7–1.4 times the default value; Figure 7F, black) were needed for the default BCM-type plasticity, while for the linear LTP curve a larger GluR1 ratio (0.59–0.92; Figure 7F, red) was needed but the NCX concentration was more variable (values ranged from 0.4 to 1.7 times the default value; Figure 7G, red). The PKA pathway coefficients were generally larger in the default BCM-type plasticity parameter sets than in the parameter sets producing the linear LTP curve (Figure 7G). Figure 7H shows the distributions of a set of coefficients, i.e. the PKC pathway, which were not correlated with the plasticity outcome within this group.
Our previous analysis showed that PKC-pathway-mediated GluR2 endocytosis was important in lower stimulation frequency protocols (Figure 3) or in protocols with large separation between pre- and post-synaptic stimuli (Figure 5). To further analyze the contribution of PKC-pathway proteins to plasticity outcomes, we repeated the analysis of Figure 7 by clustering the plasticity outcome based on the relative synaptic strength after a steady-state Ca2+ input of low amplitude (50 particles/ms; Figure 7—figure supplement 1A). As observed with the previous clustering, the GluR1 ratio and NCX concentration differed across classes (Figure 7—figure supplement 1B and C). However, in this classification, the PKA-pathway coefficient was not predictive of the plasticity outcome (Figure 7—figure supplement 1D) whereas the PKC-pathway coefficient varied across the classes (Figure 7—figure supplement 1E). Separation between BCM-like plasticity and gradual LTD was also evident within class 6’, and due to the same GluR1, PKA and NCX parameters as with the original classification (Figure 7—figure supplement 1F–H). This shows our identification of critical parameters is robust to how the classification was performed.
Taken together, our results show that alterations of the concentrations of the proteins regulating Ca2+ efflux or PKA/PKC-pathway signalling and the numbers of GluR1 and GluR2 subunits, ranging from complete absence to moderate increase (±100%), have a large effect both on the type of plasticity (LTP or LTD) and on the sensitivity of the plasticity outcome to the amplitude of the Ca2+ flux. These data suggest that neocortical post-synaptic spines may exhibit a vast set of plasticity rules by downregulation or relatively mild upregulation of their protein expression.
A parametric analysis confirms the robustness of the model
We analysed the model responses to 4xHFS and LFS protocols (as in Figure 3) under small (±10%) changes in the parameters describing the initial concentrations and reaction rates (Figure 8). As expected, most parameter changes led to small deviations from the predicted magnitudes of LTP/LTD (Figure 8, grey bars). Alterations of the initial concentration of a number of species (10 out of 47) and reaction rates (12 out of 223) resulted in a notable (>15%) amplification or attenuation of LTD (Figure 8A) or LTP (Figure 8B). The parameters affecting the LFS-induced LTD were all related to GluR1 membrane insertion or total amount of GluR1 or GluR2 (Figure 8A), while the parameters affecting the 4xHFS-induced LTP were related to NCX-mediated Ca2+ extrusion, PP1 concentration, production of cAMP by AC1, PKA buffering/deactivation, or GluR1 membrane insertion (Figure 8B). Importantly, none of the parameter changes completely abolished the LTP or LTD. Taken together, our model is robust to small alterations in initial concentrations and reaction rates, but parameters influencing the Ca2+ dynamics, GluR1 activity, or the PKA-pathway signalling can have relatively large effects on the model output.
Figure 8
The model predictions of LTP and LTD are robust to small changes in model parameters.
Values of initial concentrations (47 parameters) or reaction rates (223 parameters) were changed one at the time by −10% or +10%, and the resulting synaptic conductance 16 min after LFS (A) or 4xHFS (B) protocol was measured (NEURON RxD simulations). The initial synaptic conductance is 33.4 pS (see Figure 3A,S), although some parameter changes mildly affected this value (data not shown). The x-axis shows the post-LFS (A) or post-HFS (B) synaptic conductance, and the y-axis shows the number of parameter alterations. Majority of the parameter changes had small effect on plasticity (grey bars), but changes in initial concentrations of 10 species and 12 reaction rates caused >15% change in the amplitude of LTP or LTD — these changes are represented by black (multi-pathway parameters), blue (PKA-pathway-related parameters), and green (PKC-pathway-related parameters) bars. The underlying parameter changes are printed above the corresponding bar.
The model flexibly reproduces data from various cortical LTP/LTD experiments
The richness of the intracellular signalling machinery behind LTP and LTD poses challenges for both qualitative and quantitative comparison between results from different cell types, obtained using different stimulation protocols, or even published by different laboratories Larkman and Jack, 1995. Computational biochemically detailed models have been proposed as an absolutely reproducible tool that is particularly suited for unifying our understanding of LTP and LTD across cell types and brain regions Manninen et al., 2010. Here, we show that our model for intracellular signalling in a cortical post-synaptic spine — through the use of varying concentrations of different proteins — can be flexibly tuned to reproduce data from the experimental literature of cortical LTP/LTD. This allows one to make predictions for the differences in intracellular machineries underlying each of the experiments, leading to a more complete view of the plasticity-related signalling pathways in different cell types in the cortex and the effects of the stimulation protocol on the plasticity outcome across cortical areas.
To show the flexibility of our model, we aimed to reproduce a large amount of data on cortical plasticity across cortical areas and stimulation paradigms. We reviewed the literature of cortical plasticity, and picked eight studies where one or more types of neurons were tested using one or more stimulation protocols and the outcome was quantified using electrophysiological measurements (Table 2). These studies comprised 11 data sets that described the response of a neuron population in entorhinal cortex (EC), prefrontal cortex (PFC), barrel cortex (BC), anterior cingulate cortex (ACC), visual cortex (VC), or auditory cortex (AuC) to plasticity inducing protocols (Table 2). For each experiment in each data set, we assigned an objective function that quantified the error between the predicted LTP/LTD outcome (measured in relative synaptic conductance) and the data (typically measured in fold change of field EPSP slope). The objective functions were averaged across different time instants (10, 15, and 20 min post-stimulus-onset). We then ran a multi-objective optimisation algorithm (see Materials and methods, section 'Parameter alterations and model fitting') that aimed at finding the values for model parameters that minimised these objective functions. The fitted parameters included the amplitudes of pre-synaptic stimulation-associated fluxes of Ca2+, β-adrenergic ligand and glutamate in addition to GluR1 fraction and factors for the protein concentrations of different pathways. Here, both the Ca2+ and the neuromodulatory inputs were square-pulse injections that followed every pre-synaptic stimulation although some of the studies of Table 2 used bath application of neuromodulatory agents; however, the temporal distribution of the neuromodulators has only a small effect in our model as shown earlier (Figure 3—figure supplement 4).We ran the optimiser for 20 generations. For data sets VC-1 and VC-2, we did not find parameter sets that would fulfil all four objective functions, and therefore, we re-fitted the model for these data sets excluding the CaMKII-blocked experiments (printed in grey in Table 2).
Table 2
List of LTP/LTD experiments in the cortex.
The first column labels the experimental data set and names the underlying study. The second column shows the considered synaptic pathway and the third column shows whether the observed LTP/LTD had a pre- or post-synaptic origin. The fourth and fifth columns show the frequency (in Hz) of stimulation and the number of pulses delivered, respectively: 10 × 4 means that 10 trains of 4 pulses with 10 ms interval (100 Hz) were delivered, and likewise, 25 × 5 means that 25 trains of 5 pulses with 10 ms interval were delivered. The sixth column tells whether the data were obtained in control conditions or under additional blockers or agonists. The seventh, eighth, ninth, and tenth columns show the relative change in synaptic strength 10, 15, and 20 min after the start of the stimulus protocol and an average SD of the relative synaptic strengths — these values were approximated from the LTP/LTD curves plotted in the underlying references. The rows correspond to experiments from a given reference that are divided to 11 different experimental data sets. Within each data set, the underlying system is assumed to be otherwise similar to the control except for the applied modifier: as an example, the chemical or genetic blockade of CaMKII activity (as performed in Ma et al., 2008 and Hardingham et al., 2003) is here expected to only affect the ability of CaMKII to autophosphorylate, and the rest of the model parameters are kept fixed. The experiments printed in grey were included in the underlying study, but were excluded from the main analyses of the present work (see main text). EC – entorhinal cortex; PFC – prefrontal cortex; BC – barrel cortex; ACC – anterior cingulate cortex; VC – visual cortex; AuC – auditory cortex; CC – corpus callosum. (*): The LFS of 900 3-ms pulses at 5 Hz in data sets VC-1 and VC-2 was replaced by 180 15-ms pulses at 1 Hz to decrease computational load in the optimisation.
| Data set | Reference | Pathway | Pre/post | Freq. | Npulses | Experiment | 10 min | 15 min | 20 min | SD |
|---|---|---|---|---|---|---|---|---|---|---|
| EC-1 | Ma et al., 2008 | horizontal | mostly | 100 | 100 | control | 1.3 | 1.4 | 1.3 | 0.1 |
| post | CaMKII blocked | 1.05 | 1.02 | 0.95 | 0.07 | |||||
| without post-syn. Ca2+ | 1.05 | 1.05 | 1.1 | 0.09 | ||||||
| EC-2 | Ma et al., 2008 | ascending | mostly | 100 | 100 | control | 1.6 | 1.6 | 1.6 | 0.11 |
| post | PKA blocked | 1.4 | 1.4 | 1.4 | 0.13 | |||||
| without post-syn. Ca2+ | 1.3 | 1.4 | 1.4 | 0.13 | ||||||
| PFC-1 | Sáez-Briones et al., 2015 | CC→PFC | n/a | 312 | 156 | control | 2.0 | 1.98 | 1.9 | 0.08 |
| without -adrenergic ligand | 1.34 | 1.4 | 1.36 | 0.09 | ||||||
| PFC-2 | Flores et al., 2011 | CC→PFC | n/a | 312 | 156 | control | 1.7 | 1.6 | 1.64 | 0.12 |
| without −1-receptor agonist | 1.43 | 1.45 | 1.43 | 0.1 | ||||||
| BC | Hardingham et al., 2003 | L4→L2/3 | n/a | 5 | 10 × 4 | control | 1.35 | 1.4 | 1.3 | 0.09 |
| CaMKII mutant | 1.25 | 1.2 | 1.1 | 0.09 | ||||||
| ACC | Song et al., 2017 | L5/6 → L2/3 | post | 5 | 10 × 4 | control | 1.55 | 1.4 | 1.4 | 0.05 |
| without s845 | 1.1 | 1.05 | 1.05 | 0.07 | ||||||
| without s831 | 1.35 | 1.4 | 1.3 | 0.1 | ||||||
| PFC-3 | Zhou et al., 2013 | L2/3 → L2/3 | mostly | 0.1 | 50 | control | 1.3 | 1.4 | 1.4 | 0.14 |
| post | without −1-receptor agonist | 1.1 | 1.2 | 1.2 | 0.13 | |||||
| VC-1 | Kirkwood et al., 1997 | L4 → L3 | n/a | 5 | 10 × 4 | (CTR, HFS) | 1.3 | 1.26 | 1.26 | 0.07 |
| (adult) | (without CaMKII, HFS) | 1.02 | 1.02 | 1.02 | 0.02 | |||||
| 5 | 900* | (CTR, LFS) | n/a | 0.95 | 0.95 | 0.05 | ||||
| (without CaMKII, LFS) | n/a | 0.88 | 0.93 | 0.03 | ||||||
| VC-2 | Kirkwood et al., 1997 | L4 → L3 | n/a | 5 | 10 × 4 | (CTR, HFS) | 1.2 | 1.18 | 1.18 | 0.05 |
| (4–5 w) | (without CaMKII, HFS) | 1.07 | 1.09 | 1.08 | 0.03 | |||||
| 5 | 900* | (CTR, LFS) | n/a | 0.79 | 0.82 | 0.03 | ||||
| (without CaMKII, LFS) | n/a | 0.82 | 0.89 | 0.03 | ||||||
| AuC-1 | Kotak et al., 2007 | L6 → L5 | n/a | 1 | 25 × 5 | LTP-expressing cells | 1.98 | 1.58 | 1.93 | 0.19 |
| AuC-2 | Kotak et al., 2007 | L6 → L5 | n/a | 1 | 25 × 5 | LTD-expressing cells | 0.77 | 0.68 | 0.67 | 0.09 |
We found groups of parameter sets that fit within one standard deviation (SD) on average from the target values of synaptic conductance for each data set of Table 2 (Figure 9A–C). There were differences in the numbers of acceptable parameter sets between the data sets due to differences in the postulated strength of the LTP/LTD, the number of experiments, and the SD of the post-stimulus synaptic conductance (Figure 9D). The data sets EC-1 and BC were particularly challenging to fit (<0.2% of the parameter sets tested by the optimiser gave an acceptable fit; Figure 9D). By contrast, the data set AuC-2 was the easiest to fit (3.9% of the parameter sets were acceptable; Figure 9D).
Figure 9 with 11 supplements see all
The model can be fit to LTP/LTD data from different cortical areas.
(A) The model could be fit to LTP/LTD data from data sets EC-1 (top), EC-2, PFC-1, PFC-2, BC, ACC, PFC-3, AuC-1, and AuC-2 (bottom). The curves represent the model predictions of the best-fit parameter sets, and the dots represent the experimental data from Table 2. For data sets other than AuC-1 and AuC-2, several experiments with various chemical agents or genetic mutations were performed for each neuron population: these are ordered as in Table 2 (e.g., in data set EC-1, purple (1 st experiment) corresponds to control, blue (2nd experiment) to CaMKII-blocked experiment, and green (3rd experiment) to the experiment where post-synaptic Ca2+ was blocked). (B) The model could not be fit to the complete LTP/LTD data from data sets VC-1 (top) and VC-2 (bottom). The best parameter sets correctly predicted the LTP/LTD in up to two experiments (e.g., the selected parameter sets reproduce the HFS data with and without CaMKII inhibitor, but failed to reproduce the LFS data). (C) The model could be fit to the LTP/LTD data from data sets VC-1 (top) and VC-2 (bottom) when CaMKII-blocked experiments were ignored. The vertical bars in (B) and (C) represent the SD from the experimental data. (D) Proportion of accepted parameter sets across the 20 generations of multi-objective optimisation (20’000 parameter sets in total) in each data set. (E) Box plots of selected parameters in the acceptable parameter sets of data sets EC-1 and EC-2 (three left-most pairs) and AuC-1 and AuC-2 (right-most pair). Values of and are linearly scaled such that the values 0 and 1 correspond to 0 and double the original value of the underlying parameters, respectively (CaM and CaMKII for , and R, Gs, AC1, and AC8 for , see Materials and methods, section 'Parameter alterations and model fitting'). The medians were significantly different in the compared data sets (U-test, p-value<0.001).
The obtained parameters reflect the pathways needed for the type of plasticity. For example, the LTP of synapses of the horizontal but not those of the ascending pathway to EC were blocked by CaMKII inhibition, while the LTP of synapses of the ascending pathway were blocked by PKA inhibition (Ma et al., 2008). This is reflected in the obtained parameter sets: the parameter controlling CaM and CaMKII concentrations () was significantly larger (U-test, p-value<0.001) in the horizontal-pathway (data set EC-1) synapse models, while the parameter controlling upstream PKA-pathway proteins () was significantly larger in the models reproducing the data from the ascending pathway (data set EC-2; Figure 9E). The GluR1-GluR2 ratio, in turn, was not significantly different between the two data sets (Figure 9E). As a contrasting example, our model predicts a large variety of parameters that reproduce the LTP and LTD of data sets AuC-1 and AuC-2 but, consistent with the results of Figure 7, the GluR1 ratio was significantly larger in the model parameters fitted to the data from LTP-expressing neurons than from LTD-expressing neurons (Figure 9E). The complete graphs of parameter value distributions in the 11 data sets and the parameter set producing the best fit (Figure 9A) for each data set are shown in Supplementary data, Figure 9—figure supplements 1–11. Taken together, our model-fitting experiment shows that the model can be fit to many types of multi-condition plasticity data — without altering the reaction rates — and that the resulting predictions of the underlying protein concentrations reflect the mechanisms proposed by the experimental studies.
The models obtained by fitting the initial concentrations to data provide an important tool for predicting the outcome of plasticity under various stimulus protocols and chemical agents. We carried out additional simulations with the obtained models using HFS protocol. The models fitted for data from EC (data sets EC-1 and EC-2, [Ma et al., 2008]), BC (Hardingham et al., 2003), ACC (Song et al., 2017), and LTP-expressing auditory cortical neurons (data set AuC-1, Kotak et al., 2007) predicted a steady increase in response to HFS, while models fitted for other cortical data predicted a mixture of LTP, LTD, and no change (Figure 10A). Furthermore, to obtain experimentally testable predictions for the dependency of the plasticity outcome in different cortical areas on the intracellular signalling, we simulated the models from each data set with the corresponding stimulus protocol with CaMKII, PKA, or PKC blockade. The inhibition of CaMKII impaired the LTP in data set EC-1 (horizontal pathway) but had little or no effect on the plasticity in other cortical areas (Figure 10B). The lack of effect of CaMKII blockade on the ascending pathway of EC — an experiment which was not included in the fitting of the model (Table 2) — validates the underlying models in this aspect since similar results were observed in Ma et al., 2008. Moreover, the similarities in the predicted effect of CaMKII-, PKA-, and PKC-pathway blockades between the two models of CC→PFC synapses (PFC-1 and PFC-2; Figure 10B–D) serve as an additional validation of these models. The inhibition of PKA impaired LTP in all cortical areas, except for LTP in the horizontal pathway of the EC (EC-1; Figure 10C). Our models also predicted that LTP of the CC→PFC pathway and LFS-induced LTP in VC can be effectively weakened or even be transformed to a mild LTD by PKA blockade (Figure 10C). Finally, our models predicted that PKC inhibition transformed all forms of LTD (LFS-induced LTD in VC-1 and VC-2; LTD in AuC-2) into LTP and impaired certain forms of LTP (LTP in EC-2; LTP in AuC-1) (Figure 10D). Taken together, our results suggest that almost all forms of post-synaptic plasticity in the cortex are likely to be PKA-dependent, and that many types of cortical plasticity are also influenced by CaMKII and PKC activity. Our results highlight the need for additional chemical or genetic manipulations to be done when experimenting on cortical plasticity in order to correctly reveal the intracellular signalling cascades in the post-synaptic spine.
Figure 10
The models describing plasticity in different cortical areas predict diverse responses to modified stimulation protocol and stimulation under chemical blockers.
(A) The predicted responses of the 20 best models in each data set to HFS (100 pulses at 100 Hz) stimulation. (B–D) The predicted responses of the 20 best models in each data set to the applied stimulation protocol (see Table 2) when CaMKII (B), PKA (C), or PKC (D) activity was blocked (red) or under control condition (black).
Discussion
We built a single-compartment model describing the major post-synaptic signalling pathways leading to LTP and LTD in the cortex. We showed that our model reproduced conventional types of LTP and LTD, where an HFS-induced increase in GluR1 can increase the synaptic conductance (LTP) and an LFS-induced endocytosis of GluR2 can decrease it (LTD; Figure 3) and reproduced STDP data from visual cortical layer 2/3 pyramidal cells (Figure 5). Our model explains how different forms of plasticity depend on the concentrations of PKA- and PKC-pathway proteins (Figure 7). We also showed that our model can be fit to explain the pathway dependencies of various types of neocortical LTP/LTD data published in the literature by altering the magnitude of Ca2+ and ligand fluxes and the concentrations of post-synaptic proteins regulating the Ca2+ efflux and PKA- and PKC-pathway dynamics (Figure 9). Our fitted models provide a powerful tool for testing hypotheses on the effects of chemical or genetic manipulations on the LTP and LTD in different cortical regions (Figure 10).
Role of GluR2 in synaptic plasticity in the neocortex
GluR2 subunits are highly expressed in neocortical neurons (Kondo et al., 1997), and their endocytosis mediates (or, at minimum, is correlated with) synaptic depression in many cortical regions such as ACC (Toyoda et al., 2007), VC (Heynen et al., 2003), and PFC (Van den Oever et al., 2008). Previous intracellular signalling-based models of neocortical LTP/LTD exist, but they do not take into account the contributions of GluR2 subunit. For example, in three previous models (D'Alcantara et al., 2003; Castellani et al., 2005; Honda et al., 2013), S831-phosphorylation-mediated LTP was described in a fashion similar to our model (although more approximations were made), but the phosphorylation site S845 was assumed to be basally phosphorylated and LTD was caused by modest Ca2+ inputs that led to PP1 or PP2B-mediated dephosphorylation of S845. Although there is support for this order of events (Lee and Kirkwood, 2011; Diering et al., 2016), newer findings have confirmed the low degrees of phosphorylation of both S831 and S845 at a basal state in cortical cells, especially in synaptic spines (Diering et al., 2016). A recent phenomenological model has also shed light on the dependency of cortical STDP on the pairing interval and location of the synapse (Ebner et al., 2019), but their mechanisms are not specific to any particular AMPAR or NMDAR subunit. To analyse the contributions of GluR2 subunits to neocortical LTP/LTD, we included in our model the signalling pathways leading to both phosphorylation-mediated exocytosis of GluR1 and endocytosis of GluR2. Our model could thus be used to study not only PKA-mediated LTP or PKC-mediated LTD but also their co-effects and co-dependencies.
Implications of the study
The modelling results of the present work give rise to experimentally testable predictions. For example, our STDP model, when stimulated without β-adrenergic ligands, suggests that at near-zero pairing-intervals the magnitude of the depression may be decreased or even switched to mild LTP (Figure 5J). In many experimental studies (including Seol et al., 2007), the type and magnitude of plasticity in this regime of STDP is not reported. We also predict that a mild LTP (24–60% LTP; class 4 in Figure 7) can be obtained through many differently weighted interactions of PKA and PKC pathways and Ca2+ extrusion strengths (Figure 7F–H, Figure 9—figure supplement 1 and 2, Figure 9—figure supplement 5–7).Importantly, this is the regime of a wealth of experimental LTP data (Tsumoto, 1990; Table 2), which is consistent with the great diversity of LTP mechanisms observed in the neocortex (Feldman, 2009). Based on our simulated data (Figure 10), we suggest that in order to correctly characterise the mechanisms behind LTP of especially this magnitude, both experiments that activate the PKA pathway and experiments that block or activate the PKC pathway should be carried out. It is also important to know whether and to what degree GluR1 and GluR2 subunits are present in the synapse, since the balance of GluR1 and GluR2 subunits seems to be a determinant parameter permitting certain types of plasticity while prohibiting others (Figures 7B and 9E, and Figure 3—figure supplement 1).
A key challenge in the study of synaptic plasticity is the diversity of LTP/LTD observed across the cell types in the brain (Granger and Nicoll, 2014). Differences in the transcriptome have been proposed as one of the sources for this variability (Lisachev and Shtark, 2018). We believe our model can be used to explain some of the discrepancies in the experimental data in this regard and expand the understanding of possible molecular contributors to LTP/LTD. For example, it is known that activation of PKA pathway by dopamine or noradrenaline in PFC pyramidal neurons increases the synaptic conductance through GluR1 membrane insertion (Sun et al., 2005; Xu et al., 2010). Our model is in agreement with this (Figure 3), but it also proposes that the LTP can be impaired by over- or underexpression of many involved proteins, such as AC1, I-1 (inhibitor of PP1), PDE4, PDE1, GluR1, and CaM, and even alterations in ATP concentration (see Figure 8B). Small differences in the concentrations of a number of such contributing proteins are likely to cause significant alterations to LTP observed in different brain areas and cell types.
Due to the inclusion of three major LTP/LTD pathways in the neocortex, our model provides a more accurate means than earlier models for exploring how the Ca2+ dynamics in the spine affects the plasticity outcome in many stimulation protocols, STDP in particular. Our model suggests that the plasticity outcome of the STDP protocol is strongly correlated with the total amount of Ca2+ entering the post-synaptic spine, and less so with the peak Ca2+ flux (Figure 6—figure supplement 1).The total amount of Ca2+ influx could thus provide a better biomarker for plasticity than the previously considered amplitude and duration of the Ca2+ transient (Evans and Blackwell, 2015).
Validity of the results and limitations of the study
Our model of total synaptic conductance of the post-synaptic spine is based upon a number of assumptions. First, the prediction of a large increase of conductance that follows the replacement of GluR2 subunits at the membrane by GluR1 subunits (e.g., Figure 3) is based upon the findings on differences in single-channel conductances of different types of AMPAR tetramers in hippocampal neurons (Oh and Derkach, 2005). Following Oh and Derkach, 2005, we assumed that CaMKII-phosphorylation of S831 only increases the conductance of GluR1 homomers and not that of GluR1/GluR2 heteromers, although also heteromers have been observed to increase their conductance in the presence of transmembrane AMPAR regulatory proteins (Kristensen et al., 2011). Second, we assumed a random tetramerization procedure in which each of the four subunits in the tetramer may be either GluR1 or GluR2 subunit. Traditionally, AMPARs were thought to assemble as dimers of like dimers, that is, that first GluR1s and GluR2s assemble into homomeric R1-R1 and R2-R2 dimers and R1-R2 heterodimers and that these three types of dimers only assemble into tetramers with a dimer of its own kind (Gan et al., 2015). However, recent findings of heterotetramers with only one GluR1 subunit (Zhao et al., 2019) challenge this model. To show that our results were not dependent on the details of this process, we reproduced our results using the alternative (dimer of like dimers) tetramer formation rule. Using a slightly modified GluR1-GluR2 balance (35:65), this model reproduced HFS-induced LTP and LFS-induced LTD (Figure 3—figure supplement 2A) as well as the neuromodulator-gated STDP (Figure 3—figure supplement 2A).In summary, our model predictions were not dependent on the assumptions on the tetramer formation rule.
Our model reproduced the qualitative results of STDP of layer 2/3 pyramidal cells in visual cortex being gated by neuromodulators, but there were quantitative differences. When acetylcholine was present, our model predicted a prominent decrease in GluR2 membrane-expression regardless the pairing interval (Figure 5I), which caused a notable LTD for very large pairing intervals (Figure 5J), whereas the experimental data showed attenuation of the depression for large inter-stimulus intervals (Seol et al., 2007). This discrepancy is likely caused by processes allowing slower time-scale (>50 ms) interaction between the pre- and post-synaptic stimulus that are either not included (e.g., Ca2+-induced Ca2+ release or cAMP-dependence of HCN channels) or not adequately strong in the multi-compartmental model. For example, the contributions of voltage-gated Ca2+ channels and SK channels to the neuron electrophysiology may be large (Mäki-Marttunen et al., 2017; Mäki-Marttunen et al., 2018) — to this end, we showed here that the SK currents are amplified by the subsequent pulses stimulating the post-synaptic neuron and that this is one factor increasing the LTD for large ISIs (Figure 6). Note that this prediction of lowered synaptic strength for large absolute ISIs is not to be considered a basal synaptic state under spontaneous activity since the amplitude of the LTD is significantly decreased both by the removal of cholinergic neuromodulation (Figure 6J) and a decrease of stimulating frequency (data not shown). On the other hand, mechanisms lacking from the biochemical model (e.g., voltage-dependence of the Ca2+-extrusion rate of NCX Weber et al., 2002) could also impede our results in this matter. Some aspects of cellular physiology could therefore be better represented if we incorporated both biochemical signalling and multicompartmental Hodgkin-Huxley-type modelling into the simulations, as done in modelling studies of persistent neuron firing (Neymotin et al., 2016) and astrocyte electrophysiology (Savtchenko et al., 2018).
The quality of the model fitting to experimental data in Figure 9 is restricted by the fact that not all of the LTP/LTD data in Table 2 were confirmed to have a post-synaptic origin. This may be the key source of discrepancy in the fitting of the model to the CaMKII-blocked data from Kirkwood et al., 1997 (Figure 9B), since CaMKII activation at the pre-synaptic spine may lead to EPSC potentiation through an increase in neurotransmitter release (Ninan and Arancio, 2004). This scenario is supported by Seol et al., 2007 where S831-deficient mice were observed to show normal post-synaptic LTP in the VC.
Outlook
Our results on interactions of the different pathways in post-synaptic spines including both GluR1 and GluR2 subunits provide valuable insights on the contributions of protein expression on the plasticity of the synapse. Previously, synaptic plasticity outcomes in the cortex have been conjectured to depend on the type of the post-synaptic cell type, in addition to the timing and frequency of the applied stimuli and dendritic filtering properties (Bi and Poo, 1998; Sjöström et al., 2001). Our model provides a way to analyse exactly which aspects of PKA-, PKC- and CaMKII-pathway signalling underlie these cell-type-dependent differences in synaptic plasticity. Combining our biochemically detailed model with biophysically detailed models from different cortical areas will provide models with better predictive power in the future. Moreover, our model can be used for initial testing of hypotheses concerning dysfunctions (including chemical and genetic manipulations) of many intracellular signalling proteins and their role in impairments of cortical synaptic plasticity. By altering the initial concentrations or reaction rates of various species according to disease-associated functional genetics data, the model can be used to provide insights into the disease mechanisms of mental disorders that express both genetic disposition of post-synaptic signalling pathways and plasticity-related phenotypes, such as schizophrenia (Devor et al., 2017).
Materials and methods
Construction and calibration of the biochemically detailed model of post-synaptic plasticity in the cortex
Request a detailed protocolWe created a model of pathways leading from Ca2+ inputs and activation of β-adrenergic receptors, metabotropic glutamate receptors, and muscarinic acetylcholine receptors to the phosphorylation and insertion of AMPARs into the membrane. We started by using the model of Jȩdrzejewska-Szmek et al., 2017 for GluR1 phosphorylation at sites S831 and S845, which are phosphorylated by PKA and CaMKII, respectively, as a basis for our unified model. We added the mGluR- and M1 receptor-dependent pathways leading to PKC activation from Kim et al., 2013 and Blackwell et al., 2019, respectively, and adopted the PKC-dependent endocytosis of GluR2 and reinsertion to the membrane from Gallimore et al., 2018 as these pathways are critical for neocortical plasticity (Seol et al., 2007). As we included molecular species from different models and as we omitted certain molecular species that affected the dynamics of the underlying species but were not imperative for the pathways we wanted to describe, calibration of the model reactions was necessary. Following Hayer and Bhalla, 2005, we allowed the insertion and removal of GluR1 subunits to and from the membrane that depended on their state of S845 phosphorylation. We also allowed spontaneous membrane insertion of non-S845-phosphorylated GluR1; we chose the rate of this reaction so that on average one fifth of the (non-S845-phosphorylated) GluR1s were membrane-expressed in steady state, as suggested by experimental data (Oh et al., 2006). We adjusted the forward rate of GluR2 insertion to the membrane to decrease the proportion of membrane-inserted vs. internalised GluR2s (Figure 11A), following experimental data according to which 45% of GluR2s were membrane-inserted at resting conditions (Ashby et al., 2004). We also adopted the three-step CaM activation of Gallimore et al., 2018 instead of the two-step activation of Jȩdrzejewska-Szmek et al., 2017 where the reaction rates of CaM binding two Ca2+ ions were linearly dependent on the number of Ca2+ ions. The response curve for CaM activation by Ca2+ was steeper in this model (Figure 11B), which was in better accordance with recent experimental data (Hoffman et al., 2014). As a result, our model predicted a more prominent activation of CaM in response to a large influx of Ca2+ but milder activation in respose to small Ca2+ influx than the model of Jȩdrzejewska-Szmek et al., 2017 (Figure 11C).
Figure 11 with 2 supplements see all
Calibration of the model.
Black curves represent the final model, while grey lines represent predictions of models where previous model components or tentative parameter values were used. (A) Concentration of membrane-inserted GluR2 in 4xHFS when the forward rate of the membrane insertion of non-phosphorylated GluR2 was 0.0055 1/ms Gallimore et al., 2018 (grey) or 0.00025 1/ms (black). The rate 0.00025 1/ms caused a resting-state concentration of 121 nM for the membrane-bound GluR2 subunits, which is 45% of the total GluR2 concentration (270 nM). (B) Steady-state concentration of activated (bound by four Ca2+ ions) CaM in response to a prolonged Ca2+ input amplitude when the two-step (grey) or three-step (black) activation of CaM by Ca2+ was used. The x-axis shows the corresponding steady-state concentration of free Ca2+. Here, the initial concentrations of molecular species were as in Li et al., 2020, namely, 50 μM for CaM, Ng, PP2B, and CaMKII and 0 μM for all other species. Red dots show experimental data from Hoffman et al., 2014. (C) Concentration time course of non-protein-bound activated CaM (inset) or total activated CaM (main figure) in response to 4xHFS when the two-step (grey) or three-step (black) activation of CaM by Ca2+ was used. (D) Percentage of S880-phosphorylated GluR2 15 min after LFS when different forward rates of the activation of persistent PKC (kf between 0.00005 and 0.005 1/(nM ms)) were used. The value kf = 0.0005 1/(nM ms) gave a percentage of 47%, in close agreement with Ashby et al., 2004. (E–G) The dynamics of transiently active PKC (E) were not strongly influenced by the forward rate of the activation of persistent PKC (reaction 140), but those of persistently active PKC (F) and S880-phosphorylated GluR2 (G) were significantly affected. Black curves show the data corresponding to kf = 0.0005 1/(nM ms), while the grey lines show the data corresponding to kf = 0.00015 1/(nM ms) (dashed) and kf = 0.0015 1/(nM ms) (dotted). (H) Predicted responses of an isolated PKA activation model (reactions 59 and 93) to a 16 s cAMP input (dim grey background) when different values of the forward rate of PKA binding with four cAMP molecules were used. The curves show the concentration of the catalytic PKA subunit when different forward rates of PKA–cAMP binding (from bottom to top: 0.4 ×109, 1.0×10 9, 1.6 ×109, 2.2×10 9, and 2.8 ×109 1/(nM4ms)) were used. The markers show the corresponding data when the two-step PKA–cAMP binding model of Jȩdrzejewska-Szmek et al., 2017 was used. Inset: summed absolute differences between the tentative data (curves) and simulated data from the previous model (markers). The model with the forward rate of kf = 1.6 × 109 1/(nM4ms) gave the closest correspondence to the model of Jȩdrzejewska-Szmek et al., 2017. (I) Concentration of S845-phosphorylated GluR1 in response to 4xHFS when the single-step (reaction 59, black) or two-step (from Jȩdrzejewska-Szmek et al., 2017, grey) PKA–cAMP binding was used. (J) Concentration of S831-phosphorylated GluR1 in response to 4xHFS when PKC did (black) or did not (grey, overlaid) phosphorylated S831 in GluR1s.
To allow long-term activation of PKC, we adopted a persistently activated form of PKC, mediated by arachidonic acid (AA), from Kotaleski et al., 2002. We calibrated the rates of this reaction as follows. The backward rate was chosen so that approximately 90% of PKC would be active after 10 min, inspired by experimental data of Shirai et al., 1998. The forward rate was chosen so that LFS with effective PLC activation led to approximately 50% of the GluR2s being phosphorylated (Figure 11D), following experimental data (Ahmadian et al., 2004). The implications of these adjustments on the dynamics of transiently activated PKC, persistently activated PKC and GluR2 S880 phosphorylation are illustrated in Figures 11E, F and G, respectively.
We adopted the simplified, mass-action law-based PKA activation model (reaction 59; Table 3) from Williamson et al., 2009 (where it was called model ‘C’) instead of the 2-stage, linearised cAMP-binding of PKA in Jȩdrzejewska-Szmek et al., 2017 and Blackwell et al., 2019. We fitted the forward rate to data simulated with the original model (Figure 11H) to produce a longer-lasting S845 phosphorylation (typically,>10 min duration of S845 phosphorylation was observed Seol et al., 2007; Xue et al., 2014) in the 4xHFS protocol (Figure 11I). To account for the experimental observation the PKC phosphorylates GluR1 at site S831 Roche et al., 1996, we added this reaction using the rates identical to those of GluR1-S831 phosphorylation by phosphorylated CaMKII (see reactions 71–72). However, the presence of this reaction did not have a large effect on the S831 phosphorylation of GluR1 under standard conditions (Figure 11J). Finally, we introduced an immobile Ca2+ buffer with a Ca2+ binding rate of 0.0004 1/(nM ms), a release rate 20.0 1/ms, and an initial concentration of 500 μM (these values are within the range of experimental observations and values used in models Matthews and Dietrich, 2015). The model reactions are described in Table 3 and the initial concentrations are listed in Table 4.
Table 3
List of model reactions.
(A) The reaction-rate units are in 1/ms, 1/(nMms), 1/(nM2ms), 1/(nM3ms), or 1/(nM4ms), depending on the number of reactants. Reactions are grouped by similar modes of action and identical forward and backward rates. The denominators , , and represent groups of species detailed below. †: backward reaction rate proportional to [PKAc], not to [PKAc]2. (B) Groups of species as used in panel (A).
| (A) | Forw. | Backw. | Forw. | Backw. | |||
|---|---|---|---|---|---|---|---|
| ID | Reaction | Rate | Rate | ID | Reaction | Rate | Rate |
| 1 | Ca + PMCA ⇌ PMCACa | 5e-05 | 0.007 | 71 | GluR1 + ⇌ GluR1 | 2.78e-08 | 0.002 |
| 2 | PMCACa ⇌ PMCA + CaOut | 0.0035 | 0.0 | 72 | GluR1 ⇌ GluR1 + | 0.0005 | 0 |
| 3 | Ca + NCX ⇌ NCXCa | 1.68e-05 | 0.0112 | 73 | GluR1 + PKAc ⇌ GluR1 | 4e-06 | 0.024 |
| 4 | NCXCa ⇌ NCX + CaOut | 0.0056 | 0.0 | 74 | GluR1 + PP1 ⇌ GluR1 | 8.7e-07 | 0.00068 |
| 5 | CaOut + Leak ⇌ CaOutLeak | 1.5e-06 | 0.0011 | 75 | GluR1 ⇌ GluR1 + PP1 | 0.00017 | 0 |
| 6 | CaOutLeak ⇌ Ca + Leak | 0.0011 | 0.0 | 76 | GluR1 + PP1 ⇌ GluR1 | 8.75e-07 | 0.0014 |
| 7 | Ca + Calbin ⇌ CalbinC | 2.8e-05 | 0.0196 | 77 | GluR1 ⇌ GluR1 + PP1 | 0.00035 | 0 |
| 8 | L ⇌ LOut | 0.0005 | 2e-09 | 78 | GluR1 + PP2BCaMCa4 ⇌ GluR1 | 2.01e-06 | 0.008 |
| 9 | L + R ⇌ LR | 5.555e-06 | 0.005 | 79 | GluR1 ⇌ GluR1 + PP2BCaMCa4 | 0.002 | 0 |
| 10 | LR + Gs ⇌ LRGs | 6e-07 | 1e-06 | 80 | GluR1 ⇌ GluR1_memb | 2e-07 | 8e-07 |
| 11 | Gs + R ⇌ GsR | 4e-08 | 3e-07 | 81 | GluR1_S845 | ||
| 12 | GsR + L ⇌ LRGs | 2.5e-06 | 0.0005 | ⇌ GluR1_memb_S845 | 3.28e-05 | 8e-06 | |
| 13 | LRGs ⇌ LRGsbg + GsaGTP | 0.02 | 0.0 | 82 | PDE1 + CaMCa4 ⇌ PDE1CaMCa4 | 0.0001 | 0.001 |
| 14 | LRGsbg ⇌ LR + Gsbg | 0.08 | 0.0 | 83 | PDE1CaMCa4 + cAMP ⇌ PDE1CaMCa4cAMP | 4.6e-06 | 0.044 |
| 15 | + PKAc ⇌ PKAc | 8e-07 | 0.00448 | 84 | PDE1CaMCa4cAMP ⇌ PDE1CaMCa4 + AMP | 0.011 | 0.0 |
| 16 | PKAc ⇌ p + PKAc | 0.001 | 0.0 | 85 | AMP ⇌ ATP | 0.001 | 0.0 |
| 17 | ppLR + PKAc ⇌ PKAcppLR | 1.712e-05 | 0.00448 | 86 | PDE4 + cAMP ⇌ PDE4cAMP | 2.166e-05 | 0.0034656 |
| 18 | pppLR + PKAc ⇌ PKAcpppLR | 0.001712 | 0.00448 | 87 | PDE4cAMP ⇌ PDE4 + AMP | 0.017233 | 0.0 |
| 19 | ppppLR + Gi ⇌ ppppLRGi | 0.00015 | 0.00025 | 88 | + ⇌ PKAc | 2.5e-07 | 8e-05 |
| 20 | ppppLRGi ⇌ ppppLRGibg + GiaGTP | 0.000125 | 0.0 | 89 | PKAc ⇌ pPDE4 + PKAc | 2e-05 | 0.0 |
| 21 | pppp ⇌ pppp + Gibg | 0.001 | 0.0 | 90 | pPDE4 ⇌ PDE4 | 2.5e-06 | 0.0 |
| 22 | p ⇌ | 2.5e-06 | 0.0 | 91 | pPDE4 + cAMP ⇌ pPDE4cAMP | 0.000433175 | 0.069308 |
| 23 | pp ⇌ p | 2.5e-06 | 0.0 | 92 | pPDE4cAMP ⇌ pPDE4 + AMP | 0.3446674 | 0.0 |
| 24 | R + PKAc ⇌ PKAcR | 4e-08 | 0.00448 | 93 | PKAcAMP4 ⇌ PKAr + 2*PKAc | 0.00024 | 2.55e-05 |
| 25 | pR + PKAc ⇌ PKAcpR | 4e-07 | 0.00448 | 94 | Ca + fixedbuffer ⇌ fixedbufferCa | 0.0004 | 20.0 |
| 26 | ppR + PKAc ⇌ PKAcppR | 4e-06 | 0.00448 | 95 | Glu ⇌ GluOut | 0.0005 | 2e-10 |
| 27 | pppR + PKAc ⇌ PKAcpppR | 0.0004 | 0.00448 | 96 | Ca + PLC ⇌ PLCCa | 4e-07 | 0.001 |
| 28 | ppppR + Gi ⇌ ppppRGi | 7.5e-05 | 0.000125 | 97 | GqaGTP + PLC ⇌ PLCGqaGTP | 7e-07 | 0.0007 |
| 29 | ppppRGi ⇌ ppppRGibg + GiaGTP | 6.25e-05 | 0.0 | 98 | Ca + PLCGqaGTP ⇌ PLCCaGqaGTP | 8e-05 | 0.04 |
| 30 | GsaGTP ⇌ GsaGDP | 0.01 | 0.0 | 99 | GqaGTP + PLCCa ⇌ PLCCaGqaGTP | 0.0001 | 0.01 |
| 31 | GsaGDP + Gsbg ⇌ Gs | 0.1 | 0.0 | 100 | PLCCa + Pip2 ⇌ PLCCaPip2 | 3e-08 | 0.01 |
| 32 | GiaGTP ⇌ GiaGDP | 0.000125 | 0.0 | 101 | PLCCaPip2 ⇌ PLCCaDAG + Ip3 | 0.0003 | 0.0 |
| 33 | GiaGDP + Gibg ⇌ Gi | 0.00125 | 0.0 | 102 | PLCCaDAG ⇌ PLCCa + DAG | 0.2 | 0.0 |
| 34 | GsaGTP + AC1 ⇌ AC1GsaGTP | 3.85e-05 | 0.01 | 103 | PLCCaGqaGTP + Pip2 ⇌ PLCCaGqaGTPPip2 | 1.5e-05 | 0.075 |
| 35 | AC1 + CaMCa4 ⇌ AC1 | 6e-06 | 0.0009 | 104 | PLCCaGqaGTPPip2 ⇌ PLCCaGqaGTPDAG + Ip3 | 0.25 | 0.0 |
| 36 | + ATP ⇌ | 1e-05 | 2.273 | 105 | PLCCaGqaGTPDAG ⇌ PLCCaGqaGTP + DAG | 1.0 | 0.0 |
| 37 | AC1GsaGTPCaMCa4ATP | 106 | Ip3degrad + PIkinase ⇌ Ip3degPIk | 2e-06 | 0.001 | ||
| ⇌ cAMP + AC1GsaGTPCaMCa4 | 0.02842 | 0.0 | 107 | Ip3degPIk ⇌ PIkinase + Pip2 | 0.001 | 0.0 | |
| 38 | + ⇌ AC1Gsa | 6.25e-05 | 0.01 | 108 | PLC ⇌ PLC + GqaGDP | 0.012 | 0.0 |
| 39 | ⇌ cAMP + | 0.002842 | 0.0 | 109 | GqaGTP ⇌ GqaGDP | 0.001 | 0.0 |
| 40 | AC1GiaGTPCaMCa4ATP | 110 | GqaGDP ⇌ Gqabg | 0.01 | 0.0 | ||
| ⇌ cAMP + AC1GiaGTPCaMCa4 | 0.0005684 | 0.0 | 111 | Ca + DGL ⇌ CaDGL | 0.000125 | 0.05 | |
| 41 | AC1CaMCa4ATP ⇌ cAMP + AC1CaMCa4 | 0.005684 | 0.0 | 112 | DAG + CaDGL ⇌ DAGCaDGL | 5e-07 | 0.001 |
| 42 | AC8 + CaMCa4 ⇌ AC8CaMCa4 | 1.25e-06 | 0.001 | 113 | DAGCaDGL ⇌ CaDGL + 2AG | 0.00025 | 0.0 |
| 43 | CaM + 2*Ca ⇌ CaMCa2 | 1.7e-08 | 0.035 | 114 | Ip3 ⇌ Ip3degrad | 0.01 | 0.0 |
| 44 | + Ca ⇌ | 1.4e-05 | 0.228 | 115 | 2AG ⇌ 2AGdegrad | 0.005 | 0.0 |
| 45 | + Ca ⇌ | 2.6e-05 | 0.064 | 116 | DAG + DAGK ⇌ DAGKdag | 7e-08 | 0.0008 |
| 46 | CaM + Ng ⇌ NgCaM | 2.8e-05 | 0.036 | 117 | DAGKdag ⇌ DAGK + PA | 0.0002 | 0.0 |
| 47 | CaM + PP2B ⇌ PP2BCaM | 4.6e-06 | 1.2e-06 | 118 | Ca + PKC ⇌ PKCCa | 1.33e-05 | 0.05 |
| 48 | CaMCa + PP2B ⇌ PP2B | 4.6e-05 | 1.2e-06 | 119 | PKCCa + DAG ⇌ PKCt | 1.5e-08 | 0.00015 |
| 49 | PP2BCaM + 2*Ca ⇌ PP2BCaMCa2 | 1.7e-07 | 0.35 | 120 | Glu + MGluR ⇌ MGluR_Glu | 1.68e-08 | 0.0001 |
| 50 | CaMCa4 + CK ⇌ CKCaMCa4 | 1e-05 | 0.003 | 121 | MGluR_Glu ⇌ MGluR_Glu_desens | 6.25e-05 | 1e-06 |
| 51 | 2*CKCaMCa4 ⇌ Complex | 1e-07 | 0.01 | 122 | Gqabg + MGluR_Glu ⇌ MGluR_Gqabg_Glu | 9e-06 | 0.00136 |
| 52 | CKpCaMCa4 + CKCaMCa4 ⇌ pComplex | 1e-07 | 0.01 | 123 | MGluR_Gqabg_Glu ⇌ GqaGTP + MGluR_Glu | 0.0015 | 0.0 |
| 53 | CK + Complex ⇌ CK + pComplex | 1e-07 | 0.0 | 124 | GluR2 + PKC ⇌ GluR2 | 4e-07 | 0.0008 |
| 54 | 2*Complex ⇌ Complex + pComplex | 1e-05 | 0.0 | 125 | GluR2 ⇌ GluR2 + PKC | 0.0047 | 0 |
| 55 | Complex + pComplex ⇌ 2*pComplex | 3e-05 | 0.0 | 126 | GluR2 + PP2A ⇌ GluR2 | 5e-07 | 0.005 |
| 56 | CKpCaMCa4 ⇌ CaMCa4 + CKp | 8e-07 | 1e-05 | 127 | GluR2 ⇌ GluR2 + PP2A | 0.00015 | 0 |
| 57 | CKp + PP1 ⇌ CKp | 4e-09 | 0.00034 | 128 | GluR2 ⇌ GluR2_memb | 0.00024545 | 0.0003 |
| 58 | CKp ⇌ PP1 + CK | 8.6e-05 | 0.0 | 129 | GluR2_S880 ⇌ GluR2_memb_S880 | 0.0055 | 0.07 |
| 59 | PKA + 4*cAMP ⇌ PKAcAMP4 | 1.6e-15 | 6e-05 | 130 | ACh + M1R ⇌ AChM1R | 9.5e-08 | 0.0025 |
| 60 | Epac1 + cAMP ⇌ Epac1cAMP | 3.1e-08 | 6.51e-05 | 131 | Gqabg + AChM1R ⇌ AChM1RGq | 2.4e-05 | 0.00042 |
| 61 | I1 + PKAc ⇌ I1PKAc | 1.4e-06 | 0.0056 | 132 | Gqabg + M1R ⇌ M1RGq | 5.76e-07 | 0.00042 |
| 62 | I1PKAc ⇌ Ip35 + PKAc | 0.0014 | 0.0 | 133 | ACh + M1RGq ⇌ AChM1RGq | 3.96e-06 | 0.0025 |
| 63 | Ip35 + PP1 ⇌ Ip35PP1 | 1e-06 | 1.1e-06 | 134 | AChM1RGq ⇌ GqaGTP + AChM1R | 0.0005 | 0.0 |
| 64 | Ip35 + PP2BCaMCa4 ⇌ Ip35PP2B | 9.625e-05 | 0.33 | 135 | ACh ⇌ | 0.006 | 0 |
| 65 | Ip35PP2B ⇌ I1 + PP2B | 0.055 | 0.0 | 136 | Ca + PLA2 ⇌ CaPLA2 | 6e-07 | 0.003 |
| 66 | PP1PP2BCaMCa4 ⇌ PP1 + PP2BCaMCa4 | 0.0015 | 0.0 | 137 | CaPLA2 + Pip2 ⇌ CaPLA2Pip2 | 2.2e-05 | 0.444 |
| 67 | GluR1 + PKAc ⇌ GluR1 | 4.02e-06 | 0.024 | 138 | CaPLA2Pip2 ⇌ CaPLA2 + AA | 0.111 | 0.0 |
| 68 | GluR1 ⇌ GluR1 + PKAc | 0.006 | 0 | 139 | AA ⇌ Pip2 | 0.001 | 0.0 |
| 69 | GluR1 + CK ⇌ GluR1 | 2.224e-08 | 0.0016 | 140 | PKCt + AA ⇌ PKCp | 5e-09 | 1.76e-07 |
| 70 | GluR1 ⇌ GluR1 + CK | 0.0004 | 0 |
| (B) | |||||||
|---|---|---|---|---|---|---|---|
| ∈ {LR, pLR} | (, , ) ∈ { (_CKpCam, _S831, CKpCaMCa4), (_PKCt, | ||||||
| ∈ {LR, pLR, ppLR, pppLR, R, pR, ppR, pppR} | _S831, PKCt), (_PKCp, _S831, PKCp), (_S845_CKpCam, _S845_S831, | ||||||
| (, ) ∈ { (LRGibg, LR), (RGibg, R) } | CKpCaMCa4), (_S845_PKCt, _S845_S831, PKCt), (_S845_PKCp, | ||||||
| ∈ {LR, R, pR} | _S845_S831, PKCp), (_memb_CKpCam, _memb_S831, CKpCaMCa4), | ||||||
| ∈ {LR, pLR, ppLR, pR, ppR} | (_memb_PKCt, _memb_S831, PKCt), (_memb_PKCp, _memb_S831, PKCp), | ||||||
| (, ) ∈ { (GsaGTP, GsaGTPCaMCa4), (GsaGTPGiaGTP, | (_memb_S845_CKpCam, _memb_S845_S831, CKpCaMCa4), | ||||||
| GsaGTPGiaGTPCaMCa4), (, CaMCa4) } | (_memb_S845_PKCt, _memb_S845_S831, PKCt), (_memb_S845_PKCp, | ||||||
| (, ) ∈ { (AC1GsaGTPCaMCa4, AC1GsaGTPCaMCa4ATP), | _memb_S845_S831, PKCp) } | ||||||
| (AC1GsaGTPGiaGTPCaMCa4, AC1GsGiCaMCa4ATP), (AC1GiaGTPCaMCa4, | (, ) ∈ { (_S831, _S831_PKAc), (_memb_S831, | ||||||
| AC1GiaGTPCaMCa4ATP), (AC1CaMCa4, AC1CaMCa4ATP), (AC8CaMCa4, | _memb_S831_PKAc) } | ||||||
| AC8CaMCa4ATP) } | (, ) ∈ { (_S845, _S845_PP1), (_memb_S845, | ||||||
| (, , ) ∈ { (GiaGTP, AC1GsaGTP, GTPGiaGTP), (GiaGTP, | _memb_S845_PP1) } | ||||||
| AC1CaMCa4, GTPCaMCa4), (AC1GiaGTP, GsaGTP, GTPGiaGTP) } | (, ) ∈ { (_S845_PP1, ), (_memb_S845_PP1, | ||||||
| (, ) ∈ { (AC1GsGiCaMCa4ATP, AC1GsaGTPGiaGTPCaMCa4), | _memb) } | ||||||
| (AC8CaMCa4ATP, AC8CaMCa4) } | (, ) ∈ { (_S845_S831, _S845_S831_PP1), (_S831, | ||||||
| (, ) ∈ { (CaMCa2, CaMCa3), (PP2BCaMCa2, | _S831_PP1), (_memb_S845_S831, _memb_S845_S831_PP1), | ||||||
| PP2BCaMCa3) } | (_memb_S831, _memb_S831_PP1) } | ||||||
| (, ) ∈ { (CaMCa3, CaMCa4), (PP2BCaMCa3, | (, ) ∈ { (_S845_S831_PP1, _S845), (_S845_S831_PP1, | ||||||
| PP2BCaMCa4) } | _S831), (_S831_PP1, ), (_memb_S845_S831_PP1, _memb_S845), | ||||||
| (, ) ∈ { (2, CaMCa2), (4, CaMCa4) } | (_memb_S845_S831_PP1, _memb_S831), (_memb_S831_PP1, | ||||||
| ∈ {pCaMCa4, CaMCa4} | _memb) } | ||||||
| (, ) ∈ { (, PP1), (CaMCa4, CaMCa4PP1) } | (, ) ∈ { (_S845, _S845_PP2B), (_S845_S831, | ||||||
| (, ) ∈ { (PP1, ), (CaMCa4PP1, CaMCa4) } | _S845_S831_PP2B), (_memb_S845, _memb_S845_PP2B), | ||||||
| (, ) ∈ { (, CaMCa4), (PP1, P2BCaMCa4) } | (_memb_S845_S831, _memb_S845_S831_PP2B) } | ||||||
| ∈ {CaMCa4, P2BCaMCa4} | (, ) ∈ { (_S845_PP2B, ), (_S845_S831_PP2B, | ||||||
| (, ) ∈ { (, _PKAc), (_memb, _memb_PKAc) } | _S831), (_memb_S845_PP2B, _memb), (_memb_S845_S831_PP2B, | ||||||
| (, ) ∈ { (_PKAc, _S845), (_S831_PKAc, _S845_S831), | _memb_S831) } | ||||||
| (_memb_PKAc, _memb_S845), (_memb_S831_PKAc, | ∈ {, _PKAc, _CKCam, _CKpCam, _CKp, _PKCt, _PKCp, | ||||||
| _memb_S845_S831) } | _S831, _S831_PKAc, _S831_PP1} | ||||||
| (, , ) ∈ { (, CaMCa4, _CKCam), (, p, | ∈ {, _CKCam, _CKpCam, _CKp, _PKCt, _PKCp, _S831, | ||||||
| _CKp), (_S845, CaMCa4, _S845_CKCam), (_S845, p, _S845_CKp), | _PP1, _S831_PP1, _PP2B, _S831_PP2B} | ||||||
| (_memb, CaMCa4, _memb_CKCam), (_memb, p, _memb_CKp), | (, , ) ∈ { (PKAc, PDE4, PDE4), (PDE4cAMP, PKAc, | ||||||
| (_memb_S845, CaMCa4, _memb_S845_CKCam), (_memb_S845, p, | _PDE4_cAMP) } | ||||||
| _memb_S845_CKp) } | (, ) ∈ { (PDE4, ), (_PDE4_cAMP, cAMP) } | ||||||
| (, , ) ∈ { (_CKCam, _S831, CaMCa4), (_CKp, _S831, | (, ) ∈ { (GqaGTP, ), (CaGqaGTP, Ca) } | ||||||
| p), (_S845_CKCam, _S845_S831, CaMCa4), (_S845_CKp, _S845_S831, | (, , ) ∈ { (, t, _PKCt), (, p, _PKCp), | ||||||
| p), (_memb_CKCam, _memb_S831, CaMCa4), (_memb_CKp, _memb_S831, | (_memb, t, _memb_PKCt), (_memb, p, _memb_PKCp) } | ||||||
| p), (_memb_S845_CKCam, _memb_S845_S831, CaMCa4), | (, , ) ∈ { (_PKCt, _S880, t), (_PKCp, _S880, p), | ||||||
| (_memb_S845_CKp, _memb_S845_S831, p) } | (_memb_PKCt, _memb_S880, t), (_memb_PKCp, _memb_S880, p) } | ||||||
| (, , ) ∈ { (, CKpCaMCa4, _CKpCam), (, | (, ) ∈ { (_S880, _S880_PP2A), (_memb_S880, | ||||||
| PKCt, _PKCt), (, PKCp, _PKCp), (_S845, CKpCaMCa4, | _memb_S880_PP2A) } | ||||||
| _S845_CKpCam), (_S845, PKCt, _S845_PKCt), (_S845, PKCp, | (, ) ∈ { (_S880_PP2A, ), (_memb_S880_PP2A, | ||||||
| _S845_PKCp), (_memb, CKpCaMCa4, _memb_CKpCam), (_memb, PKCt, | _memb) } | ||||||
| _memb_PKCt), (_memb, PKCp, _memb_PKCp), (_memb_S845, CKpCaMCa4, | ∈ {, _PKCt, _PKCp} | ||||||
| _memb_S845_CKpCam), (_memb_S845, PKCt, _memb_S845_PKCt), | ∈ {, _PP2A} | ||||||
| (_memb_S845, PKCp, _memb_S845_PKCp) } | |||||||
Table 4
List of initial concentrations of molecular species.
All non-mentioned species have an initial concentration of 0 nM.
| Species | Conc. (nM) | Species | Conc. (nM) | Species | Conc. (nM) | |||
|---|---|---|---|---|---|---|---|---|
| CaOut | extracell. Ca2+ | 1900000 | AMP | adenosine monophosphate | 980 | Pip2 | phosphatidylinositol 4,5-bisphosphate | 24000 |
| Leak | leak channels | 2000 | Ng | neurogranin | 20000 | PIkinase | phosphatidylinositol kinase | 290 |
| Calbin | calbindin | 150000 | CaM | calmodulin | 60000 | Ip3degPIk | Ip3-bound PI kinase | 400 |
| CalbinC | Ca2+-bound calbindin | 15000 | PP2B | protein phosphatase 2B | 2300 | PKC | protein kinase C | 15000 |
| LOut | extracell. β-adr. ligand | 2500000 | CK | CaMKII | 23000 | DAG | diacylglycerol | 90 |
| Epac1 | Epac1 | 500 | PKA | protein kinase A | 6400 | DAGK | DAG kinase | 300 |
| PMCA | Ca2+ pump | 22000 | I1 | inhibitor-1 | 2200 | DGL | DAG lipase | 1600 |
| NCX | Ca2+ exchanger | 540000 | PP1 | protein phosphatase 1 | 1600 | CaDGL | Ca2+-bound DAG lipase | 250 |
| L | β-adrenergic ligand | 10 | GluR1 | AMPAR subunit type 1 | 180 | DAGCaDGL | Ca2+-and DAG-bound DAG lipase | 90 |
| R | β-adrenergic receptor | 1600 | GluR1_memb | membrane-inserted GluR1 | 90 | Ip3degrad | degraded Ip3 | 600 |
| Gs | S-type G-protein | 13000 | PDE4 | phosphodiesterase type 4 | 670 | GluR2 | AMPAR subunit type 2 | 14 |
| Gi | I-type G-protein | 2600 | fixedbuffer | immobile buffer | 500000 | GluR2_memb | membrane-inserted GluR2 | 256 |
| AC1 | adenylyl cyclase type 1 | 430.0 | mGluR | metab. glutamate receptor | 800 | PP2A | protein phosphatase 2A | 500 |
| ATP | adenosine triphosphate | 2000000 | GluOut | extracell. glutamate | 1000000 | M1R | acetylcholine receptor M1 | 450 |
| AC8 | adenylyl cyclase type 8 | 370 | Gqabg | Q-type G-protein | 1400 | PLA2 | phospholipase A2 | 1000 |
| PDE1 | phosphodiesterase type 1 | 12000 | PLC | phospholipase C | 250 |
Throughout the work, we simulated the signalling pathways in a single compartment representing a dendritic spine of size 0.5 μm3. In reality, some of the molecular species are prevalently present in the cytosol, some attached to the membrane, some in the extracellular medium in an immediate vicinity to the membrane, and others outside the cell further away from the synaptic cleft (free in the extracellular medium or sequestered to other cells). As commonly done in the field, we solved this problem by introducing species that represent a molecular species confined in a particular location: reactions 1–6 describe the extrusion of Ca2+ from the cytosol into the extracellular medium, reactions 8, 95, and 135 describe the escape of ligands from the vicinity of the synapse, and reactions 80–81 and 128–129 for the translocation of the AMPARs to/from the membrane (Table 3). All stimulations start after 4040 s of simulation without inputs, which is sufficient for attaining a steady state for all species (Figure 11—figure supplement 1).
Statistical model for numbers of AMPAR tetramers at the membrane and the total synaptic conductance
Request a detailed protocolAMPARs have different conductances depending on their subunit composition and phosphorylation state (Oh and Derkach, 2005), but it is challenging to take this into account in models that include a large number of receptor subunits. In our model, AMPAR subunits GluR1 and GluR2 can be in one of 21 or five states, respectively, when counting all the different phosphorylation states and bonds with other molecules (Table 3), which leads to 284 possible types of tetramers. This makes it virtually impossible to model the dynamics of AMPAR tetramer assembly using the mass-action law-based approach where the concentration of each type of species is monitored (Michalski and Loew, 2012). To avoid this problem, we used a statistical model that estimated the numbers and types of different types of AMPAR tetramers given the numbers of GluR1 and GluR2 subunits located at the membrane.
We assumed that the composition of AMPAR tetramers is random such that there is no preference of one type of subunit being more likely to bind with any other type of subunit. Thus, the probability of a tetramer being a GluR1 homomer without any S831-phosphorylated subunits is approximately
(1)
where N1 and N2 are the numbers of GluR1 and GluR2 subunits bound to the membrane, respectively, and is the number of S831-phosphorylated GluR1 subunits at the membrane (note that the subunits are included in all GluR1 subunits, that is, ). Accordingly, the probabilities of a tetramer being a GluR1 homomer with at least one S831-phosphorylated subunit, a GluR2 homomer, or a heteromer, are approximately:
The number of membrane-bound tetramers that the N1 GluR1 subunits and N2 GluR2 subunits at the membrane can form is . Here, we ignore the unpaired subunits by estimating that — we also disregard the states of the non-membrane-bound subunits as they are not assumed to contribute to the synaptic conductance. This gives us approximate values for expected numbers of different types of tetramers on the membrane:
These estimates allow us to determine the total maximal synaptic conductance as the sum of the numbers of these tetramers multiplied with the corresponding single-channel conductances:
(5)
The single-channel conductance values 12.4 pS, 18.9 pS, 2.2 pS, and 2.5 pS are taken from experimental data (Oh and Derkach, 2005).
Modelling the Ca2+ inputs and neuromodulatory inputs
Request a detailed protocolWe modelled the neurotransmission to the post-synaptic spine as fluxes of Ca2+ ions, β-adrenergic ligand, glutamate, and acetylcholine (labelled as Ca, L, Glu, and ACh, respectively, in Table 3). We used various stimulation paradigms: In sections 'Ca2+ activates multiple pathways that regulate the post-synaptic plasticity in cortical PCs' and 'The model predicts multimodal, protein concentration- and neuromodulation-dependent rules of plasticity', long-lasting, single pulses of input species were applied. In sections 'High-frequency stimulation (HFS) causes LTP and low-frequency stimulation (LFS) causes LTD in GluR1-GluR2-balanced synapses' and 'A parametric analysis confirms the robustness of the model', we used the following repeated stimulus protocols: HFS — 100 pulses of Ca2+ (3 ms), repeated at 100 Hz; 4xHFS — 4 trains of HFS, separated by 3 s of quiescence; LFS — 900 pulses of Ca2+ (3 ms), repeated at 5 Hz. Unless otherwise stated, each Ca2+ pulse was accompanied by a 3 ms pulse of β-adrenergic ligand, glutamate, and acetylcholine. The activation of cholinergic and noradrenergic terminals by electrical stimulation is supported by experimental data in, e.g., slices of mouse prefrontal cortex (Mundorf et al., 2001). In section 'The model flexibly reproduces data from various cortical LTP/LTD experiments', the same approach was used, but the frequencies and numbers of repetitions of the inputs were taken from the experiments (see Table 2).
In section 'Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses' (and Figure 3—figure supplement 3 of section 'High-frequency stimulation (HFS) causes LTP and low-frequency stimulation (LFS) causes LTD in GluR1-GluR2-balanced synapses'), we used a multicompartmental model of a layer 2/3 pyramidal cell (Markram et al., 2015) (L23_PC_cADpyr229_1) to determine the amplitudes and time courses of the Ca2+ inputs conducted by NMDARs when different stimulus patterns were applied. This model included the fast Na+ current , Shaw-related K+ current , muscarinic K+ current , and hyperpolarization-activated cyclic nucleotide-gated current in the apical dendrite. The axo-somatic region contained all these (except for ) as well as the low and high-voltage activated Ca2+ currents and , the small-conductance Ca2+-dependent K+ current , the transient K+ current , and the persistent Na+ and K+ currents and in the axo-somatic region (Markram et al., 2015). We placed 10 post-synaptic spines, consisting of a 0.5 μm long and 0.1 μm thick neck and a 0.4 μm long and 0.4 μm thick head, to the proximal apical dendrite (250–300 μm from the soma). For an analysis of the effects of synapse location on Ca2+ flux and plasticity in cortical pyramidal cells, see Ebner et al., 2019. Each spine was equipped with the AMPA–NMDA synapse model of Hay and Segev, 2015 using the NMDA gating mechanism of Spruston et al., 1995 and an adjustment in the pre-synaptic resource update (Mäki-Marttunen et al., 2019). The ten synapses were synchronously stimulated but the glutamate release was probabilistic, the release events at different synapses being independent of each other as in Hay and Segev, 2015). We set the maximal AMPAR conductance to 0.001 μS, a value typically used in computational neuron models. We set the maximal conductance of NMDARs to 0.0032 μS to compensate for the lack of the slow component in the model of the NMDA current (Markram et al., 2015) — a value significantly smaller or larger than this abolished the LTP or LTD, respectively, in our STDP model, while the AMPAR conductance was a less crucial parameter (Figure 11—figure supplement 2).We estimated the numbers of Ca2+ ions entering into the post-synaptic spine across time, and used these numbers as the input to the biochemical model. In section 'Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses', following the experiments of Seol et al., 2007, we used extracellular [Mg2+] of 1.0 mM and a 1 Hz pre-synaptic stimulation, where each pulse was paired with a burst of post-synaptic stimulus currents (Seol et al., 2007). For Figure 3—figure supplement 3 of section 'High-frequency stimulation (HFS) causes LTP and low-frequency stimulation (LFS) causes LTD in GluR1-GluR2-balanced synapses', following the experiments of Heusler et al., 2000 we used [Mg2+] of 1.3 mM and solely pre-synaptic stimulation of one of the two stimulus protocols: 6xHFSt — 10 bursts of 4 pulses (at 100 Hz), repeated every 100 ms, and the whole train repeated 6 times every 10 s; LFS-1Hz — 1800 pulses delivered at a frequency of 1 Hz.
The effects of LTP/LTD on the size of Ca2+ inputs were not considered in this work.
Parameter alterations and model fitting
Request a detailed protocolIn sections 'The model predicts multimodal, protein concentration- and neuromodulation-dependent rules of plasticity' and 'The model flexibly reproduces data from various cortical LTP/LTD experiments', we altered the initial concentrations of many proteins to explore the parameter space or to perform model fitting. We chose to fit protein concentrations instead of reaction rates, since the reaction rates can be considered to be the same across cell types while the protein expression is known to be cell-type and age-dependent. This is analogous to fitting maximal conductances that correlate with ion-channel densities in Hodgkin-Huxley-type models instead of the ion-channel activation and inactivation curve parameters as is usually done in the fitting of biophysically detailed neuron models. The concentrations of upstream PKA-pathway proteins R (β-adrenergic receptor), Gs, AC1, and AC8 were varied in proportion using a factor parameter , and, likewise, the concentrations of upstream PKC-pathway proteins mGluR, M1, Gq, and PLC using a factor parameter . Furthermore, in section 'The model flexibly reproduces data from various cortical LTP/LTD experiments', CaM and CaMKII were altered in proportion by a factor , phosphatases PP1 and PP2B by a factor , and phosphodiesterases PDE1 and PDE4 by a factor . In both sections, the rapidity of Ca2+ extrusion was varied by altering the concentration of NCX, and in section 'The model flexibly reproduces data from various cortical LTP/LTD experiments', the concentrations of PKA and PKC were varied in addition to the upstream proteins — these concentrations were varied within the interval from 0 to double the original value. For the multi-objective optimisation in section 'The model flexibly reproduces data from various cortical LTP/LTD experiments', we used the Python implementation (published by the authors of Bahl et al., 2012) of the non-dominated sorting genetic algorithm II (NSGA-II) (Deb et al., 2002) with population size 1000. To restrict to physiologically realistic Ca2+ dynamics, we disregarded the data where free Ca2+ concentrations rose above 2 μM for one or more levels of Ca2+ input in section 'The model predicts multimodal, protein concentration- and neuromodulation-dependent rules of plasticity'. In a similar manner, in section 'The model flexibly reproduces data from various cortical LTP/LTD experiments', we introduced an objective function that penalised parameter sets that produced Ca2+ transients larger than 2 μM.
Simulation software and code accessibility
Request a detailed protocolFor deterministic simulations of intracellular signalling, we used the NEURON simulator with the reaction-diffusion (RxD) extension (McDougal et al., 2013). For stochastic simulations, we used NeuroRD software (https://github.com/neurord). In both types of simulations, we used adaptive time-step integration methods. The NEURON simulator was also used for simulating the multicompartmental model of layer 2/3 pyramidal cell in section 'Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses'. The full model along with the fitting and data-analysis algorithms (Python scripts) that were used in this study are publicly available in ModelDB at http://modeldb.yale.edu/260971.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files. Simulation scripts can be found at http://modeldb.yale.edu/260971.
-
ModelDBID 260971. Biochemically detailed model of LTP and LTD in a cortical spine (Maki-Marttunen et al 2020).
References
-
Biophysical and phenomenological models of multiple spike interactions in spike-timing dependent plasticityInternational Journal of Neural Systems 16:79–97.https://doi.org/10.1142/S0129065706000524
-
Automated optimization of a reduced layer 5 pyramidal cell model based on experimental dataJournal of Neuroscience Methods 210:22–34.https://doi.org/10.1016/j.jneumeth.2012.04.006
-
AMPA receptor subunits expression and phosphorylation in cingulate cortex in rats following esophageal acid exposureNeurogastroenterology & Motility 25:973–e776.https://doi.org/10.1111/nmo.12233
-
Molecular mechanisms underlying striatal synaptic plasticity: relevance to chronic alcohol consumption and seekingEuropean Journal of Neuroscience 49:768–783.https://doi.org/10.1111/ejn.13919
-
Serotonin 5-HT1A receptors regulate AMPA receptor channels through inhibiting Ca2+/calmodulin-dependent kinase II in prefrontal cortical pyramidal neuronsThe Journal of Biological Chemistry 277:36553–36562.https://doi.org/10.1074/jbc.M203752200
-
Different forms of synaptic plasticity in somatosensory and motor Areas of the neocortexThe Journal of Neuroscience 15:5324–5333.https://doi.org/10.1523/JNEUROSCI.15-07-05324.1995
-
Bidirectional synaptic plasticity as a consequence of interdependent Ca2+-controlled phosphorylation and dephosphorylation pathwaysEuropean Journal of Neuroscience 17:2521–2528.https://doi.org/10.1046/j.1460-9568.2003.02693.x
-
A fast and elitist multiobjective genetic algorithm: nsga-iiIEEE Transactions on Evolutionary Computation 6:182–197.https://doi.org/10.1109/4235.996017
-
Calcium: amplitude, duration, or location?The Biological Bulletin 228:75–83.https://doi.org/10.1086/BBLv228n1p75
-
Synaptic mechanisms for plasticity in neocortexAnnual Review of Neuroscience 32:33–55.https://doi.org/10.1146/annurev.neuro.051508.135516
-
Assembly of AMPA receptors: mechanisms and regulationThe Journal of Physiology 593:39–48.https://doi.org/10.1113/jphysiol.2014.273755
-
Expression mechanisms underlying long-term potentiation: a postsynaptic view, 10 years onPhilosophical Transactions of the Royal Society B: Biological Sciences 369:20130136.https://doi.org/10.1098/rstb.2013.0136
-
Dendritic excitability and gain control in recurrent cortical microcircuitsCerebral Cortex 25:3561–3571.https://doi.org/10.1093/cercor/bhu200
-
Molecular switches at the synapse emerge from receptor and kinase trafficPLOS Computational Biology 1:e20.https://doi.org/10.1371/journal.pcbi.0010020
-
Neurogranin alters the structure and calcium binding properties of calmodulinJournal of Biological Chemistry 289:14644–14655.https://doi.org/10.1074/jbc.M114.560656
-
β-adrenergic signaling broadly contributes to LTP inductionPLOS Computational Biology 13:e1005657.https://doi.org/10.1371/journal.pcbi.1005657
-
Neurophysiological effects of bitopertin in schizophreniaJournal of Clinical Psychopharmacology 37:447–451.https://doi.org/10.1097/JCP.0000000000000722
-
Mechanisms and significance of spike-timing dependent plasticityBiological Cybernetics 87:373–382.https://doi.org/10.1007/s00422-002-0351-0
-
On expression patterns and developmental origin of human brain regionsPLOS Computational Biology 12:e1005064.https://doi.org/10.1371/journal.pcbi.1005064
-
Emergence of dynamic memory traces in cortical microcircuit models through STDPJournal of Neuroscience 33:11515–11529.https://doi.org/10.1523/JNEUROSCI.5044-12.2013
-
Dopaminergic modulation of cortical plasticity in Alzheimer's disease patientsNeuropsychopharmacology 39:2654–2661.https://doi.org/10.1038/npp.2014.119
-
Subcellular interactions between parallel fibre and climbing fibre signals in purkinje cells predict sensitivity of classical conditioning to interstimulus intervalIntegrative Physiological & Behavioral Science 37:265–292.https://doi.org/10.1007/BF02734249
-
Mechanism of Ca2+/calmodulin-dependent kinase II regulation of AMPA receptor gatingNature Neuroscience 14:727–735.https://doi.org/10.1038/nn.2804
-
A model for synaptic development regulated by NMDA receptor subunit expressionJournal of Computational Neuroscience 24:1–20.https://doi.org/10.1007/s10827-007-0036-8
-
Synaptic plasticity: hippocampal LTPCurrent Opinion in Neurobiology 5:324–334.https://doi.org/10.1016/0959-4388(95)80045-X
-
AMPA receptor regulation during synaptic plasticity in Hippocampus and neocortexSeminars in Cell & Developmental Biology 22:514–520.https://doi.org/10.1016/j.semcdb.2011.06.007
-
BookLong-term potentiation-associated gene expression: involvement of the tumour proteinIn: Stuchlik A, editors. The Hippocampus-Plasticity and Functions. IntechOpen. pp. 1–53.https://doi.org/10.5772/intechopen.73219
-
Postsynaptic signal transduction models for long-term potentiation and depressionFrontiers in Computational Neuroscience 4:152.https://doi.org/10.3389/fncom.2010.00152
-
Pathological plasticity in fragile X syndromeNeural Plasticity 2012:1–12.https://doi.org/10.1155/2012/275630
-
Buffer mobility and the regulation of neuronal calcium domainsFrontiers in Cellular Neuroscience 9:48.https://doi.org/10.3389/fncel.2015.00048
-
Reaction-diffusion in the NEURON simulatorFrontiers in Neuroinformatics 7:28.https://doi.org/10.3389/fninf.2013.00028
-
Modulation of synaptic plasticity in the cortex needs to understand all the playersFrontiers in Synaptic Neuroscience 9:2.https://doi.org/10.3389/fnsyn.2017.00002
-
Catecholamine release and uptake in the mouse prefrontal cortexJournal of Neurochemistry 79:130–142.https://doi.org/10.1046/j.1471-4159.2001.00554.x
-
Spine Ca2+ signaling in spike-timing-dependent plasticityJournal of Neuroscience 26:11001–11013.https://doi.org/10.1523/JNEUROSCI.1749-06.2006
-
Extrasynaptic membrane trafficking regulated by glur1 serine 845 phosphorylation primes ampa receptors for ltpThe Journal of Biological Chemistry 281:752–758.https://doi.org/10.1074/jbc.M509677200
-
Control of synaptic plasticity in deep cortical networksNature Reviews Neuroscience 19:166–180.https://doi.org/10.1038/nrn.2018.6
-
Distinct effects of fatty acids on translocation of gamma- and epsilon-subspecies of protein kinase CJournal of Cell Biology 143:511–521.https://doi.org/10.1083/jcb.143.2.511
-
Dendritic excitability and synaptic plasticityPhysiological Reviews 88:769–840.https://doi.org/10.1152/physrev.00016.2007
-
Dendritic glutamate receptor channels in rat hippocampal CA3 and CA1 pyramidal neuronsThe Journal of Physiology 482 ( Pt 2:325–352.https://doi.org/10.1113/jphysiol.1995.sp020521
-
Long-term depression requires postsynaptic AMPA GluR2 receptor in adult mouse cingulate cortexJournal of Cellular Physiology 211:336–343.https://doi.org/10.1002/jcp.20940
-
Long-term potentiation and depression in the cerebral neocortexThe Japanese Journal of Physiology 40:573–593.https://doi.org/10.2170/jjphysiol.40.573
-
Requirement of an allosteric kinetics of NMDA receptors for spike timing-dependent plasticityJournal of Neuroscience 28:3310–3323.https://doi.org/10.1523/JNEUROSCI.0303-08.2008
-
Prefrontal cortex AMPA receptor plasticity is crucial for cue-induced relapse to heroin-seekingNature Neuroscience 11:1053–1058.https://doi.org/10.1038/nn.2165
-
Phosphorylation of AMPA receptors: mechanisms and synaptic plasticityMolecular Neurobiology 32:237–250.https://doi.org/10.1385/MN:32:3:237
Decision letter
-
Harel Z ShouvalReviewing Editor; University of Texas Medical School at Houston, United States
-
Michael J FrankSenior Editor; Brown University, United States
-
Harel Z ShouvalReviewer; University of Texas Medical School at Houston, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This is a very ambitious project, that models many of the different molecular interactions that can affect synaptic plasticity. The inclusion of neuromodulators is important, novel and can account for recent data about the role of neuromodulators. The statistical model of AMPA receptors is interesting novel and practical. It is approximation that should be validated experimentally. The authors take seriously the role of data. Due to the models high dimensionality and complexity it is hard to currently validate many assumptions. The significance of this paper is not necessarily the specific assumptions being made, rather it is that many pathway previously ignored are now included, and are shown to significantly contribute.
Decision letter after peer review:
Thank you for submitting your article "A unified computational model for cortical post-synaptic plasticity" for consideration by eLife. Your article has been reviewed by two peer reviewers, including Harel Z Shouval as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Michael Frank as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
As the editors have judged that your manuscript is of interest, but as described below that additional experiments are required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is 'in revision at eLife'. Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)
Summary:
This is a very ambitious project, that models many of the different molecular interactions that can affect synaptic plasticity. The inclusion of neuromodulators is important and relatively novel. The statistical model of AMPA receptors is interesting, novel and practical. The authors take seriously the role of data. It is quite a heroic effort, but on the other hand the model is hard to verify, it is hard to extract more general principles, and the predictions of the model are limited. This limits the overall impact of the work.
The model is very complex, and almost none of the reaction coefficients are known. There are so many parameter combinations that can yield the same results. What do we actually learn from this? This is an essential part of this model, not much can be done to address this.
Spine model should be used in all simulations of plasticity protocol. It is not clear and how the magnitude of the calcium currents through the NMDAR, and of efflux through diffusion and pumps are calibrated.
It is not clear if all initial conditions are at the fixed points.
The dynamics of the neuromodulators are odd and sometimes unclear. They should be justified, explained better or corrected if these dynamics cannot be justified. Does the application of neuromodulator to the bath on its own change the efficacies?
Essential revisions:
1) All simulations of plasticity protocols should use the spine model, and it should be explained how the NMDA currents in the spine model are calibrated.
2) All simulations should start from a fixed point of the model.
3) The neuromodulator dynamics should be better explained and justified. Does it make sense to have square neuromodulator pulses? What are the neuromodulator dynamics in the STDP experiments? Especially must take into account that in the Seol et al., 2007 paper neuromodulators are bath applied. Does bath application of neuromodulators, alone, without stimulation change the efficacies?
Reviewer #1:
This is a very ambitious project, almost heroic. The attempt is to carry out true quantitative modeling of many of the signal transduction pathways involved in LTP/LTD in neocortex. The work here is very detailed and extensive, and there are many interesting and novel components to this work. However, I am not yet convinced that this goal is feasible. The model is very complex, there is a huge number of unknown parameters and I am not sure that I understand all of the methods and therefore do not know if they are appropriate.
The good:
1) The inclusion of neuromodulators is important and relatively novel.
2) The statistical model of AMPA receptors is interesting novel and practical.
3) The authors take seriously the role of data.
The bad:
1) The model is very complex, and almost none of the reaction coefficients are known. There are so many parameter combinations that can yield the same results. What do we actually learn from this?
2) All simulation of plasticity including LFS and HFS should be carried out on the basis if the spine model in Neuron. It is not clear to me though how that was calibrated.
3) The neuromodulator dynamics, which are pulses for every Ca pulse, are problematic, and not justified. It is unclear to me what the NM dynamics are in the subsection “Paired-pulse stimulus protocol induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses”
Detailed points – Major:
The good:
4) The inclusion of neuromodulators and their effects on the PKC and PKA pathways is important. Most previous models have ignored these effects, but results in adult neocortical slices have proven that these pathways are essential for LTP and LTD.
Most previous models have ignored the role of the different subunits in the AMPA receptor heteromers. Modeling of these receptors directly as independent species in a mass action approach is nearly impossible due to the huge number of possible combinations. The statistical model proposed is a novel idea that has significant practical advantages. There are assumptions that go into this, which the authors acknowledge. Nevertheless, it is a novel and sensible approach. I think a separate work could simply be based on this model and the statistical testing of its validity.
These authors take seriously the role of data and test the model (with different parameters) against many different experimental results.
The major problems:
5) This is a very complex model with 140 reactions and 47 types of elementary molecules which could be in many states. A tiny fraction of the reaction coefficients have been measured, even if they have it might not be relevant for in vivo conditions. The determination of these coefficients in the paper does not seem to depend on measured coefficients, but instead seems to depend on these coefficients being appropriate for matching experimental results. As we know, and as is shown for example in Figure 7, there is a very large space of parameter combinations that can produce very similar experimental results. What do we actually learn from all of this then? This problem is actually made clear in Figure 7, what are the authors trying to clarify in this very complex figure?
6) The real complete model is the one that uses the Calcium transients from a Neuron simulation, subsection “Paired-pulse stimulus protocol induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses” and Figure 6. Only here there is a possibility that realistic calcium transients are used, and in principle this should have been used from all comparisons to data. Not only STDP protocols have real synapses, HFS and LFS stimuli also use real synapses. However, it is not clear to me how this synapse model was calibrated. Peak Ca influx rate here is much larger than for the other cases, more than an order of magnitude larger than in Figure 3 and much larger than in Figures 4, 5 as well. How was this calibrated? Was this based on the number of NMDAR and this influx through each? If so, what are these numbers? An alternative way of calibrating the spine model is to use estimates of calcium influx from Ca imaging, for example from Sabatini and Svoboda, 2002. A third option is to use whatever works, however then even this elementary component of the model is not based on biophysical realism. Several details and references are given in the subsection “Modelling the Ca2+ inputs and neuromodulatory inputs”, but this is still not clear to me. It would also be useful to see what the spine voltages are and how they affect the Mg block. Why is there more calcium influx at -30 ms vs. only presynaptic stimulation? Is the -30 ms measured from the first or last spike in the 4-spike train? What is it about a spike that occurred 30 ms prior to the presynaptic stimulus that affects calcium influx?
7) Neuromodulator dynamics. It is important that neuromodulators have been included here. However, the assumptions about their dynamics do not make sense to me. In several sections neuromodulators pulses are assumed to follow the Ca pulses? What is the logic here, that stimulation of axons also causes neuromodulator release? This clearly does not seem to match experiments like Seol et al., 2007, where neuromodulators are bath applied and should just be at a constant level. Is there any evidence that in other slice experiments neuromodulators are indeed release at every pulse? What about the culture experiments? It is also not clear to me what is done with neuromodulators in the subsection “Paired-pulse stimulus protocol induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses”, where Ca transients were taken from the spine model. Here it is clear to me that constant Neuromodulator levels should be used as in the experiment.
8) Are initial conditions steady at steady state for each parameter combination? It is not clear to me if all plasticity simulations are started at the steady levels of the system for the given parameter set? Are they?
Reviewer #2:
1) The authors develop a complicated model of the biochemical pathways underlying LTP and LTD. It is quite a heroic effort, but on the other hand the model is hard to verify, it is hard to extract more general principles, and the predictions of the model are limited. This limits the overall impact of the work.
2) Some of the design decisions are hard to follow.
For instance, why was the CaM activation made steeper?
3) I also did not understand how the model fitting was done by changing initial concentrations (Materials and methods). Changing the reaction rates would be a more conventional way. I wonder how these concentrations develop in the absence of stimulation. Do they stay the same, or do they have to be clamped to certain values? Yet in the subsection “The model flexibly reproduces data from various cortical LTP/LTD experiments” other variables are changed to fit the data (what are 'factors for the protein concentrations'?).
4) The model's complexity make it difficult to understand it's properties. For instance, does CaMKII act as a switch, and is the expression essentially binary (O'Connor and Wang)? Does it fit the observations of Nevian and Sakmann?
5) The STDP curves look odd, with no below baseline LTD for short negative intervals.
6) Does the last sentence of the subsection “High-frequency stimulation causes LTP and low-frequency stimulation causes LTD in GluR1-GluR2-balanced synapses” really imply a causal relation, so that GuR2 endocytosis leads to potentiation or depression? If so, the mechanism was not clear to me.
7) The y-axis labels on the plots are odd. In Figure 2 they put the quantity as the plot label, and the units as axis label. The authors do it correctly on the x-axis.
In other figures other conventions are followed.
https://doi.org/10.7554/eLife.55714.sa1Author response
Essential revisions:
1) All simulations of plasticity protocols should use the spine model, and it should be explained how the NMDA currents in the spine model are calibrated.
Our STDP simulations do indeed use the electrical model to determine calcium input, and now we have added simulations of HFS and LFS using calcium inputs from the multicompartmental neuron model. We have explained the calibration of the model in more detail, and we have added a figure that justifies the chosen parameter values.
2) All simulations should start from a fixed point of the model.
We have confirmed that the simulations start from a steady state, and added a figure (Figure 11—figure supplement 1) to show this.
3) The neuromodulator dynamics should be better explained and justified. Does it make sense to have square neuromodulator pulses? What are the Neuromodulator dynamics in the STDP experiments? Especially must take into account that in the Seol et al., 2007 paper neuromodulators are bath applied. Does bath application of neuromodulators, alone, without stimulation change the efficacies?
We have added justification to the dynamics of the neuromodulators, as explained in more detail below. Furthermore, we have now simulated the same bath-application protocol as used in Seol et al., 2007. We also showed that the type of dynamics (pulse-associated or bath-applied) of the neuromodulators does not make a large difference to the results of the model. The bath application of neuromodulators alone did not cause plasticity, as we now mention in the text.
Reviewer #1:
[…] 1) The model is very complex, and almost none of the reaction coefficients are known. There are so many parameter combinations that can yield the same results. What do we actually learn from this?
We thank the reviewer for this comment. We agree with the reviewer that the model is complex when measured in numbers of reactions or equations. Nevertheless, all dynamics of the modeled molecular species follow the law of mass action, which gives the model an intrinsic modular property: if future modelers want to simplify parts of it or (as we believe to be more likely) extend it to take into account even more details, they will only have to replace a set of targeted species and their interactions with the rest of the model species without refitting the whole model. This is an important difference compared to simpler models with few parameters and non-mass-action-law type of interactions that have been particularly fitted to describe the plasticity-related phenomena.
It would be ideal to have all reaction rates explicitly fitted to data obtained from carefully constrained experiments, but such data are scarcely published. Nonetheless, rate constants for many of the reactions were constrained by data. We used off-the-shelf models of the biochemistry of PKA, PKC and CaMKII pathways where the reaction rates and initial concentrations were usually fitted directly to enzyme kinetic or FRET imaging data.
What we can learn from such a complex model is the interplay of different molecular species in the induction of post-synaptic plasticity in the neocortex. This is an important research topic in polygenic mental disorders, schizophrenia in particular, where subtle changes (due to common gene variants) in plasticity-related proteins can lead to pathological cellular and network physiology. Our model can provide biochemically detailed model predictions that are urgently needed in discovering mechanisms for these disorders. A step toward this direction is Figure 7, which, although it shows that various combinations of model parameters can result in similar plasticity curve, also pinpoints the most sensitive molecular species in different types of plasticity. Many of these sensitive species are indeed in the groups of proteins encoded by risk genes of schizophrenia (see Figure 1 in Devor et al., 2017). Like any model predictions, such biochemistry-based predictions on disease pathophysiology may be biased due to under-constrained model parameters (e.g., rate constants) and thus should be tested experimentally – until that happens, they can serve as useful and potentially ground-breaking proof-of-principle mechanisms.
2) All simulation of plasticity including LFS and HFS should be carried out on the basis if the spine model in Neuron. It is not clear to me though how that was calibrated.
Using our biophysically detailed NEURON simulations, we have now calculated the Ca2+ transients resulting from HFS and LFS stimulation protocols. In particular, we showed that our model reproduced the LTP caused by a theta-HFS protocol (6xHFSt) and the LTD caused by a 1-Hz LFS protocol (LFS-1Hz) as described in Heusler et al., 2000. We have extended the description of the model calibration for this. We now write:
“In the above analyses, we used brief square-pulse fluxes of Ca2+to the synapse model, which is a simple representation of inputs during synaptic plasticity induction protocols. […] In accordance with Figure 3 and experimental data from somatosensory cortex [Heusler et al., 2000], our model predicted that 6xHFSt induced LTP whereas LFS-1Hz induced LTD (Figure 3—figure supplement 3).”
As both reviewers noted, the biochemical model is in itself rather complex, and we try to avoid bringing extra complexity to the manuscript when not absolutely necessary. Moreover, when we talk about the proposed “unified” model that can provide explanations for a range of forms of plasticity in the neocortex, we mainly refer to the biochemical network of the synapse, not the combination of the biochemical network and the multicompartmental biophysical model (the latter of which is known to vary between pyramidal cell types and cortical areas). For these reasons, we consider some of the simulation experiments added in response to the reviewer’s comments supplementary, and thus included these results in the supplementary material. We also believe it helps the reader that we start from generic HFS-induced LTP and LFS-induced LTD rather than a specific LTP and LTD in the somatosensory cortex.
3) The Neuromodulator dynamics, which are pulses for every Ca pulse, are problematic, and not justified. It is unclear to me what the NM dynamics are in the subsection “Paired-pulse stimulus protocol induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced”.
Little is known about the dynamics of neuromodulatory inputs during different stimulation protocols. We know that norepinephrinic projections span all layers of cortex (Simpson KL, Waterhouse BD, Lin RC. Anat Rec A Discov Mol Cell Evol Biol. 2006;288(2):166‐173), and that cholinergic projections, depending on the source, may span either all layers or only deep layers of the cortex (Obermayer J, Verhoog MB, Luchicchi A, Mansvelder HD. Front Neural Circuits. 2017;11:100). Thus, it is reasonable to hypothesize that HFS and LFS, when applied to the cortex, also cause rhythmic neurotransmitter release. There indeed exists evidence that noradrenaline and acetylcholine are released in response to electrical stimulation in cortical slices (see e.g. Mundorf et al., 2001; Vizi ES, Zsilla G, Caron MG, Kiss JP. J Neurosci. 2004;24(36):7888‐7894). Thus, our model of short pulses that accompany the Ca2+ pulses is justified experimentally. Nevertheless, the amplitude of this release may vary between cortical areas and slices prepared in different ways, as some experiments (such as Seol et al., 2007) indicate that bath application of neuromodulators (or corresponding agents) is required for activation of PKA or PKC pathways – this is why we allowed the amplitudes of the stimulus-associated neuromodulatory fluxes to vary in our model fitting to experimental data (Figure 9).
To justify these choices, we now write:
“Unless otherwise stated, each Ca2+ pulse was accompanied by a 3-ms pulse of β-adrenergic ligand, glutamate, and acetylcholine. The activation of cholinergic and noradrenergic terminals by electrical stimulation is supported by experimental data in, e.g., slices of mouse prefrontal cortex [Mundorf et al., 2001].”
To show that our modelling results are not dependent on the type of neuromodulator fluxes (pulse-associated or bath-applied), we have performed the HFS and LFS simulation experiments using the same amounts of neuromodulators but with different temporal distributions. We now write:
“The HFS-induced LTP and LFS-induced LTD of Figure 3 could also be reproduced with alternative durations of neuromodulator inputs, including 10-minute bath applications (Figure 3—figure supplement 4). These results indicate that our model can reproduce HFS-induced LTP and LFS-induced LTD also when using realistic NMDAR-conducted Ca2+ transients and that these forms of plasticity are robust to the temporal profile of the neuromodulatory inputs.”
We have now clarified the neuromodulator dynamics in the subsection “Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses”. We now write:
“When added as bath application, the β-adrenergic and cholinergic ligands were simulated by prolonged injections of 50 particles/s for 10 min, starting 8 min before the STDP protocol – these neuromodulators alone (without the electric stimulation-mediated Ca2+ inputs) did not cause synaptic plasticity.”
4) The inclusion of neuromodulators and their effects on the PKC and PKA pathways is important. Most previous models have ignored these effects, but results in adult neocortical slices have proven that these pathways are essential for LTP and LTD.
Most previous models have ignored the role of the different subunits in the AMPA receptor heteromers. Modeling of these receptors directly as independent species in a mass action approach is nearly impossible due to the huge number of possible combinations. The statistical model proposed is a novel idea that has significant practical advantages. There are assumptions that go into this, which the authors acknowledge. Nevertheless, it is a novel and sensible approach. I think a separate work could simply be based on this model and the statistical testing of its validity.
We thank the reviewer for these comments. We agree that the validity of the statistical model for AMPAR conductance could be tested using detailed models of AMPAR assembly and transport, and that it should be a separate paper.
5) This is a very complex model with 140 reactions and 47 types of elementary molecules which could be in many states. A tiny fraction of the reaction coefficients have been measured, even if they have it might not be relevant for in vivo conditions. The determination of these coefficients in the paper does not seem to depend on measured coefficients, but instead seems to depend on these coefficients being appropriate for matching experimental results. As we know, and as is shown for example in Figure 7, there is a very large space of parameter combinations that can produce very similar experimental results. What do we actually learn from all of this then? This problem is actually made clear in Figure 7, what are the authors trying to clarify in this very complex figure?
We have now simplified Figure 7 to transmit the main idea – that, despite the large number of model parameters (initial concentrations of proteins), the outcome of plasticity in the model is largely determined by 3-4 groups of species, namely, the GluR1 vs. GluR2 ratio, the main Ca2+ extrusion proteins, and the concentrations of PKA and PKC pathway-related proteins. Understanding the contributions of different proteins to synaptic plasticity is important since there is a large diversity of reported types of plasticity in the cortex. Previously, we illustrated three classes of LTP/LTD curves in panels F-L, in part because not all experiments demonstrate both LTP and LTD; now, we have mainly focused on a single class, namely, the one (class 6) that produced the same LTP outcome as the model with the default parameters for a large Ca2+ input. We have also removed panels I-L as they are not imperative for the main message of our manuscript. Instead, we added a supplementary figure (Figure 7—figure supplement 1) where the classification is based on plasticity outcome for small Ca2+ inputs, which better highlights the contribution of the PKC pathway.
We now write:
“The model results for the large parameter distributions of Figure 7B–E imply that there are manifestly different combinations of parameters that lead to the same LTP/LTD outcome. […] This shows our identification of critical parameters is robust to how the classification was performed.”
To highlight the model with the default parameters instead of any arbitrary parameter set, we write:
“Our model with the standard protein concentrations (GluR1 ratio 0.5, [NCX]=0.54 mM, fPKA = fPKC = 1.0) produced a BCM-type curve in class 6 (Figure 7A, black dashed curve).”
To further clarify the take-home message of this section, we now write:
“These data suggest that neocortical post-synaptic spines may exhibit a vast set of plasticity rules by downregulation or relatively mild upregulation of their protein expression.”
This conclusion is important, as the last Results section shows the implications of this hypothesis in synapse models fitted to real data – namely, we show that the particular types of post-synaptic plasticity can indeed be obtained by concentration changes of a similar magnitude as shown here.
6) The real complete model is the one that uses the Calcium transients from a Neuron simulation subsection “Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses” and Figure 6. Only here there is a possibility that realistic calcium transients are used, and in principle this should have been used from all comparisons to data. Not only STDP protocols have real synapses, HFS and LFS stimuli also use real synapses. However, it is not clear to me how this synapse model was calibrated. Peak Ca influx rate here is much larger than for the other cases, more than an order of magnitude larger than in Figure 3 and much larger than in Figures 4, 5 as well. How was this calibrated? Was this based on the number of NMDAR and this influx through each? If so, what are these numbers? An alternative way of calibrating the spine model is to use estimates of calcium influx from Ca imaging, for example from Sabatini and Svoboda, 2002. A third option is to use whatever works, however then even this elementary component of the model is not based on biophysical realism. Several details and references are given in the subsection “Modelling the Ca2+ inputs and neuromodulatory inputs”, but this is still not clear to me. It would also be useful to see what the spine voltages are and how they affect the Mg block. Why is there more calcium influx at -30 ms vs. only presynaptic stimulation? Is the -30 ms measured from the first or last spike in the 4-spike train? What is it about a spike that occurred 30 ms prior to the presynaptic stimulus that affects calcium influx?
We agree with the reviewer that it is important that our model is responsive to realistic levels of Ca2+ inputs. To complement the simple input protocols of Figure 3 (previously Figure 4), we have now implemented the HFS and LFS protocols using real synapses in addition to the previously shown STDP protocol. We now write:
“In the above analyses, we used brief square-pulse fluxes of Ca2+ to the synapse model, which is a simple representation of inputs during synaptic plasticity induction protocols. […] In accordance with Figure 3 and experimental data from somatosensory cortex [Heusler et al., 2000], our model predicted that 6xHFSt induced LTP whereas LFS-1Hz induced LTD (Figure 3—figure supplement 3).”
To introduce this extension of the study, we added details into the Materials and methods, where we now write:
“For Figure 3—figure supplement 3 of the subsection “High-frequency stimulation (HFS) causes LTP and low-frequency stimulation (LFS) causes LTD in GluR1-GluR2-balanced synapses”, following the experiments of [Heusler et al., 2000] we used [Mg2+] of 1.3 mM and solely pre-synaptic stimulation of one of the two stimulus protocols: 6xHFSt – 10 bursts of 4 pulses (at 100 Hz), repeated every 100 ms, and the whole train repeated 6 times every 10 sec; LFS-1Hz – 1800 pulses delivered at a frequency of 1 Hz.”
In the same section, we have now also added details on the calibration of the biophysical model. We write:
“In the subsection “Paired pre- and post-synaptic stimulation induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses” (and Figure 3—figure supplement 3 of the subsection “High-frequency stimulation (HFS) causes LTP and low-frequency stimulation (LFS) causes LTD in GluR1-GluR2-balanced synapses”), we used a multicompartmental model of a layer 2/3 pyramidal cell [Markram et al., 2015] (L23 PC cADpyr229 1) to determine the amplitudes and time courses of the Ca2+ inputs conducted by NMDA receptors (NMDARs) when different stimulus patterns were applied. […] We set the maximal conductance of NMDARs to 0.0032 μS to compensate for the lack of the slow component in the model of the NMDA current [Markram et al., 2015] – a value significantly smaller or larger than this abolished the LTP or LTD, respectively, in our STDP model, while the AMPAR conductance was a less crucial parameter (Figure 11—figure supplement 2).”
The calcium in Figure 2 (earlier Figure 3) was injected as a constant influx for long duration, and therefore, the amplitudes of the intracellular calcium responses were much smaller than in HFS and LFS protocols (where only 3-ms pulses of calcium were given). In Figure 3 (earlier Figure 4), the rates of the square-pulse calcium inputs in HFS and LFS experiments were indeed in line with the observations from one of the only publications of calcium in response to 100 Hz input (Regehr WG, Tank DW. J Neurosci 1992;12(11):4202-23). Note that the calcium response to 100 Hz is not 100 times the response to a single PSP, to avoid producing unrealistic peak responses. The calcium inputs of Figure 5 (earlier Figure 6) are of the same order of magnitude as in Figures 3-4. We have now plotted the Ca2+ concentration in the spine (Figure 5E-G, Figure 6C-D) instead of the Ca2+ influx to make the comparison with experimental data easier. The Ca2+ influx and its dependence on the ISI in the STDP protocol are now illustrated in Figure 5—figure supplement 1.
We have also added panels Figure 5B-D showing the membrane potential time courses at the spine heads and the corresponding [Mg2+]-block factors. The inter-stimulus interval was indeed determined from the onset of the pre-synaptic stimulus until the onset of the first post-synaptic stimulus. We have now changed it to the interval from the onset of the pre-synaptic stimulus until the onset of the last post-synaptic stimulus, to use the same convention as experimentalists. This shifts the STDP curve 30 ms to the right. In the legend of Figure 5, we now write:
“Membrane potential at the dendritic spine when the pre-synaptic stimulation onset is 50 ms after (B), at the same time as (C), or 50 ms prior to (D) the onset of the last somatic stimulus. Inset (red): Mg2+-gate variable as a function of time, ranging from -80 ms to 140 ms in a similar manner as the data in the main panel.”
7) Neuromodulator dynamics. It is important that Neuromodulators have been included here. However, the assumptions about their dynamics do not make sense to me. In several sections Neuromodulators pulses are assumed to follow the Ca pulses? What is the logic here, that stimulation of axons also causes neuromodulator release? This clearly does not seem to match experiments like Seol et al., 2007, where neuromodulators are bath applied and should just be at a constant level. Is there any evidence that in other slice experiments neuromodulators are indeed release at every pulse? What about the culture experiments? It is also not clear to me what is done with neuromodulators in the subsection “Paired-pulse stimulus protocol induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses” where Ca transients were taken from the spine model. Here it is clear to me that constant Neuromodulator levels should be used as in the experiment.
We agree that in the subsection “Paired-pulse stimulus protocol induces PKA- and PKC-dependent spike-timing-dependent plasticity (STDP) in GluR1-GluR2-balanced synapses” it makes sense to model the same bath application of neuromodulatory agonists as done in Seol et al. We reran the simulations for Figure 5 using the same protocol of neuromodulation (10 min bath application, starting 8 min before the onset of the pairing) and updated Figure 5 with those new results (which indeed are the same as before). We have also shown that the temporal dynamics of the neuromodulation do not play a crucial role in determining the outcome of the plasticity (see our response to point #2). As it is not known what the temporal dynamics are like in vivo, we keep the pulse-like neuromodulation previously implemented in Figure 4 (now Figure 3).
We now write:
“When added as bath application, the β-adrenergic and cholinergic ligands were simulated by prolonged injections of 50 particles/s for 10 min, starting 8 min before the STDP protocol”.
To clarify that M1 receptors were activated by cholinergic ligands and that the induced LTD was dependent on the pairing interval, we rephrased the sentence as follows:
“When β-adrenergic neurotransmission was blocked but M1 receptors were activated by cholinergic ligands, the model predicted a prominent (up to 60%) decrease in GluR2 membrane expression, with little effect on GluR1 membrane expression (Figure 5K–L, purple).”
On the question on neuromodulatory inputs in slices, please see our response to point #2. We are not aware of evidence of norepinephrine or cholinergic neurotransmission in cultured cortical neurons. However, the data we previously mentioned as being from cultured neurons from auditory cortex (Kotak et al., 2007) were in fact from slices, we apologize for this error in the text. Thus, all experimental data that we reproduce in the present work were from cortical slices, except from Flores et al., 2011, which recorded field responses in prefrontal cortex in vivo.
8) Are initial conditions steady at steady state for each parameter combination? It is not clear to me if all plasticity simulations are started at the steady levels of the system for the given parameter set? Are they?
We always simulated the model for more than 1 hour before stimulation onset. This is now clarified, where we write:
“All stimulations start after 4040 sec of simulation without inputs, which is sufficient for attaining a steady state for all species (Figure 11—figure supplement 1).”
The added supplementary figure (Figure 11—figure supplement 1) shows that after one hour the concentrations change by much less than 1 molecule/sec.
Reviewer #2:
1) The authors develop a complicated model of the biochemical pathways underlying LTP and LTD. It is quite a heroic effort, but on the other hand the model is hard to verify, it is hard to extract more general principles, and the predictions of the model are limited. This limits the overall impact of the work.
We thank the reviewer for this comment. We would like to point out that by mapping model entities to biological entities, the biochemical network model is straightforward to verify, because experimentalists can measure plasticity outcome while inhibiting specific molecules to test some of our predictions. A major goal of the proposed work is to make sense of the vast and confusing literature on synaptic plasticity, in which a diversity of induction protocols in a diversity of brain regions yields sometimes contradictory plasticity outcomes. Thus, a general principle extracted from our modeling study is that a relatively small number of key molecules can explain some of this diversity (as analyzed in Figure 7 and Figure 9). Another general principle is that the balance of GluR1/GluR2 subunits is important for determining the type of plasticity (see also our response to point #6 below). On the scarcity of model predictions: we have now performed additional simulations using our combination of biophysically detailed simulations of layer 2/3 pyramidal cells with the biochemically detailed post-synaptic spine model. Our results show that the shape of the STDP curve is radically affected by the number of post-synaptic stimuli per burst (Figure 6), fewer stimuli resulting in both weaker LTP and weaker LTD in an SK-channel dependent manner. This may explain why large number of stimuli per post-synaptic burst are often used in experimental settings. Our results also suggest that STDP is more correlated with the total Ca2+ entering the post-synaptic spine than with the peak Ca2+ flux (Figure 6—figure supplement 1).
In the Results section, we now added a paragraph to present and discuss the abovementioned predictions:
“The combination of our biochemically detailed model with the biophysically detailed model of layer 2/3 pyramidal cell model provides a compelling means of hypothesis testing for cortical STDP in this cell type. […] Our model predictions also agree with the observation that the plasticity outcome is not determined by Ca2+ transient amplitude [Nevian and Sakmann, 2006], instead, our model suggests that the total Ca2+ is a better predictor of the plasticity outcome: the correlation coefficient between the post-STDP synaptic conductance and the peak Ca2+ transient amplitude (see Figure 5—figure supplement 1E) was 0.53, while that between the post-STDP synaptic conductance and the mean Ca2+ input during the inter-stimulus interval (see Figure 5—figure supplement 1F) was 0.96 (Figure 6—figure supplement 1).”
We have also rewritten parts of the Discussion to better highlight the above-mentioned take-home messages. We now write:
“We also showed that our model can be fit to explain the pathway dependencies of various types of neocortical LTP/LTD data published in the literature by altering the magnitude of Ca2+ and ligand fluxes and the concentrations of post-synaptic proteins regulating the Ca2+ efflux and PKA- and PKC-pathway dynamics (Figure 9).”
Furthermore, we write:
“It is also important to know whether and to what degree GluR1 and GluR2 subunits are present in the synapse, since the balance of GluR1 and GluR2 subunits seems to be a determinant parameter permitting certain types of plasticity while prohibiting others (Figures 7B, 9E, and Figure 3—figure supplement 1).”
We added a Discussion paragraph on the total amount of Ca2+ vs. peak Ca2+ flux in determining the outcome of plasticity:
“Due to the inclusion of three major LTP/LTD pathways in the neocortex, our model provides a more accurate means than earlier models for exploring how the Ca2+ dynamics in the spine affects the plasticity outcome in many stimulation protocols, STDP in particular. […] The total amount of Ca2+ influx could thus provide a better biomarker for plasticity than the previously considered amplitude and duration of the Ca2+ transient [Evans and Blackwell, 2015].”
To argue for the overall impact of the present work, we would like to point out that our study is the first one (as far as we are aware) to model the contributions of both GluR1 and GluR2 to synaptic plasticity in the neocortex. This is significant since both subunits are strongly expressed in the neocortex. Our statistical plasticity rule based on the numbers GluR1 and GluR2 at the membrane is an important step in this prediction pipeline. To further verify our model according to higher-level plasticity data, we have now implemented the HFS and LFS using Ca2+ transients from a biophysically detailed model (Figure 3—figure supplement 3) as previously done for the STDP experiment. On the lack of general principles, we have now simplified Figure 7 to better highlight the most central model parameters affecting the plasticity outcomes in our model.
2) Some of the design decisions are hard to follow.
For instance, why was the CaM activation made steeper?
A steeper Ca2+ activation curve has been observed in experiments. We are now showing the calcium-sensitivity of our CaM model and that of the old CaM model, overlaid with experimental data of Hoffman et al., 2014, in Figure 11B. We believe that this property of steeper CaM activation may be important for bidirectional LTP/LTD so that small amounts of calcium activate PKC but not PKA pathway, while larger calcium inputs also activate the PKA pathway – however, analyzing this particular question is out of the scope of the present work. Note that the data used to adjust the reactions was independent of the synaptic plasticity experiments we were trying to emulate. We have now addressed this:
“The response curve for CaM activation by Ca2+ was steeper in this model (Figure 11B), which was in better accordance with recent experimental data [Hoffman et al., 2014]”.
In Figure 11B, we replaced the amplitude of the Ca2+ flux in a 4xHFS protocol on the x-axis by the steady-state Ca2+ concentration. This is explained in the legend of Figure 11:
“Steady-state concentration of activated (bound by four Ca2+ ions) CaM in response to a prolonged Ca2+ input amplitude when the two-step (grey) or three-step (black) activation of CaM by Ca2+ was used. […] Red dots show experimental data from [Hoffman et al., 2014].”
3) I also did not understand how the model fitting was done by changing initial concentrations (Materials and methods). Changing the reaction rates would be a more conventional way. I wonder how these concentrations develop in the absence of stimulation. Do they stay the same, or do they have to be clamped to certain values?
Yet in the subsection “The model flexibly reproduces data from various cortical LTP/LTD experiments” other variables are changed to fit the data (what are 'factors for the protein concentrations'?).
For the model construction (Materials and methods), we changed reaction rates because experimental measurements of protein concentration are less common. On the other hand, for fitting data from various cortical regions our reasoning is that the reaction rates should be relatively constant across cells of the same type, and even across cell types. However, we know that the numbers of proteins vary dramatically depending on the metabolic state of the cell, and there are also cell-type-specific and age-dependent differences in the levels of expression of different proteins. We therefore kept the rate constants fixed and varied the concentrations of different species in the subsection “The model flexibly reproduces data from various cortical LTP/LTD experiments” to show that our model can reproduce a number of experimental data sets with the same reaction rates. This could be compared to fitting of Hodgkin-Huxley-type models, where usually the same activation and inactivation parameters are assumed for an ion-channel species in many cell types, but the maximal conductance term, relating to the number of ion channels per membrane area, is varied.
We now write:
“We chose to fit protein concentrations instead of reaction rates, since the reaction rates can be considered to be the same across cell types while the protein expression is known to be cell-type and age-dependent. This is analogous to fitting maximal conductances that correlate with ion-channel densities in Hodgkin-Huxley-type models instead of the ion-channel activation and inactivation curve parameters as is usually done in the fitting of biophysically detailed neuron models.”
The concentrations of the species are in a steady state when we start the stimulation, this is now shown in Figure 11—figure supplement 1.
By factors for protein concentrations we mean parameters that control the number of a handful of proteins. We have now clarified this:
“The concentrations of upstream PKA-pathway proteins R (β-adrenergic receptor), Gs, AC1, and AC8 were varied in proportion using a factor parameter fPKA ∈ [0,2], and, likewise, the concentrations of upstream PKC-pathway proteins mGluR, M1, Gq, and PLC using a factor parameter fPKC ∈ [0,2].”
4) The model's complexity make it difficult to understand it's properties. For instance, does CaMKII act as a switch, and is the expression essentially binary (O'Connor and Wang)? Does it fit the observations of Nevian and Sakmann?
Our CaMKII model is not bistable, but once activated, CaMKII slowly becomes dephosphorylated (see reactions 52 and 56). The rate of this deactivation is dependent on the rate of CaMCa4 unbinding from phosphorylated CaMKII (reaction 56, forward reaction rate). These rate constants for CaMKII were tuned in a previous model to match the data of (De Koninck P and Schulman H, Science 1998 Jan 9;279(5348):227-30). However, due to the autophosphorylation of CaMKII, CaMKII activation is very steep (Figure 2I).
Our STDP model agrees with the data from Nevian and Sakmann, 2006, on the effects of decreased Ca2+ channel conductivity on LTD, namely, that a blockade decreases the amplitude of LTD. Our model suggests that these effects are due to decreased SK currents for post-pre intervals. Our model also agrees with their observation that the peak of the Ca2+ transient is a poor predictor of the plasticity outcome. As for the molecular dependence of the STDP, Nevian and Sakmann, 2006, showed that the STDP they observed was dependent on mGluRs – our predictions agree that it is PLC-pathway-dependent, but in our model and in Seol et al., 2007, the activation of the PLC pathway was driven by cholinergic receptors instead of mGluRs.
We now write:
“The paired-pulse protocol of Figure 5M (involving both β-adrenergic and cholinergic neuromodulation) caused an STDP in all cases, but decreasing the SK conductance shortened the post-pre LTD window and decreased the amplitude of LTD (Figure 5E). […] Our model predictions also agree with the observation that the plasticity outcome is not determined by Ca2+ transient amplitude [Nevian and Sakmann, 2006], instead, our model suggests that the total Ca2+ is a better predictor of the plasticity outcome: the correlation coefficient between the post-STDP synaptic conductance and the peak Ca2+ transient amplitude (see Figure 5—figure supplement 1E) was 0.53, while that between the post-STDP synaptic conductance and the mean Ca2+ input during the inter-stimulus interval (see Figure 5—figure supplement 1F) was 0.96 (Figure 6—figure supplement 1).”
5) The STDP curves look odd, with no below baseline LTD for short negative intervals.
In the original report, we measured the ISI from the onset of the pre-synaptic stimulus to the onset of the first post-synaptic stimulus. We have now changed the measurement of the ISI to be consistent with experimental conventions – it is now the interval from the pre-synaptic stimulus onset to the onset of the last post-synaptic stimulus (i.e., fourth somatic pulse). This means that our STDP curve is now shifted 30 ms to the right when using the four somatic stimulus pulses, separated by 10 ms. This has the consequence of now exhibiting below baseline LTD for short negative intervals. In the legend of Figure 5 we now write:
“Membrane potential at the dendritic spine when the pre-synaptic stimulation onset is 50 ms after (B), at the same time as (C), or 50 ms prior to (D) the onset of the last somatic stimulus.”
6) Does the last sentence of the subsection “High-frequency stimulation causes LTP and low-frequency stimulation causes LTD in GluR1-GluR2-balanced synapses” really imply a causal relation, so that GluR2 endocytosis leads to potentiation or depression? If so, the mechanism was not clear to me.
Yes. Our model uses the data of Oh et al., 2005, which shows that homomeric GluR1 receptors have much larger conductance than heteromeric GluR1-GluR2 tetramers (and that the conductance of homomeric GluR1 receptors, unlike that of heteromeric GluR1-GluR2 tetramers, is further increased through S831 phosphorylation). As decreasing GluR2 subunits from the spine membrane through S880 phosphorylation can significantly increase the proportion of homomeric GluR1 tetramers, our model predicts that this can lead to increase in the total conductance although the total number of AMPAR tetramers is decreased. We have now illustrated this more clearly in panels D and H of Figure 4 (earlier Figure 5) by adding a third pair of bars showing the contribution of homomeric GluR1s to the total synaptic conductance. In the legend of Figure 4, we write:
“D, H: The fraction of membrane-inserted GluR1 over all membrane-inserted GluR subunits (left), the probability of an AMPAR tetramer being homomeric GluR1 (middle), and the relative contribution of homomeric GluR1 subunits to the total conductance (i.e., summed conductance of homomeric GluR1 tetramers divided by the summed conductance of all tetramers; right).”
7) The y-axis labels on the plots are odd. In Figure 2 they put the quantity as the plot label, and the units as axis label. The authors do it correctly on the x-axis.
In other figures other conventions are followed.
We have now written the labels so that units are always in parentheses. The quantity is written, as before, either before the units on the axis label or in the main panel.
https://doi.org/10.7554/eLife.55714.sa2Article and author information
Author details
Tuomo Mäki-Marttunen
Funding
Research Council of Norway (248828)
- Tuomo Mäki-Marttunen
- Andrew G Edwards
- Gaute T Einevoll
Horizon 2020 - Research and Innovation Framework Programme (785907)
- Gaute T Einevoll
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
UNINETT Sigma2 resources (project NN9529K) were used for simulations. Funding: Research Council of Norway (248828), European Union Horizon 2020 Research and Innovation Programme under Grant Agreement No. 785907 [Human Brain Project (HBP) SGA2].
Senior Editor
- Michael J Frank, Brown University, United States
Reviewing Editor
- Harel Z Shouval, University of Texas Medical School at Houston, United States
Reviewer
- Harel Z Shouval, University of Texas Medical School at Houston, United States
Publication history
- Received: February 3, 2020
- Accepted: July 29, 2020
- Accepted Manuscript published: July 30, 2020 (version 1)
- Version of Record published: August 13, 2020 (version 2)
- Version of Record updated: January 28, 2021 (version 3)
Copyright
© 2020, Mäki-Marttunen et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 2,234
- Page views
-
- 335
- Downloads
-
- 14
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
A unified computational model for cortical post-synaptic plasticity
eLife 9:e55714.
https://doi.org/10.7554/eLife.55714
Further reading
-
- Computational and Systems Biology
Different strains of a microorganism growing in the same environment display a wide variety of growth rates and growth yields. We developed a coarse-grained model to test the hypothesis that different resource allocation strategies, corresponding to different compositions of the proteome, can account for the observed rate-yield variability. The model predictions were verified by means of a database of hundreds of published rate-yield and uptake-secretion phenotypes of Escherichia coli strains grown in standard laboratory conditions. We found a very good quantitative agreement between the range of predicted and observed growth rates, growth yields, and glucose uptake and acetate secretion rates. These results support the hypothesis that resource allocation is a major explanatory factor of the observed variability of growth rates and growth yields across different bacterial strains. An interesting prediction of our model, supported by the experimental data, is that high growth rates are not necessarily accompanied by low growth yields. The resource allocation strategies enabling high-rate, high-yield growth of E. coli lead to a higher saturation of enzymes and ribosomes, and thus to a more efficient utilization of proteomic resources. Our model thus contributes to a fundamental understanding of the quantitative relationship between rate and yield in E. coli and other microorganisms. It may also be useful for the rapid screening of strains in metabolic engineering and synthetic biology.
-
- Computational and Systems Biology
- Neuroscience
Biological motor control is versatile, efficient, and depends on proprioceptive feedback. Muscles are flexible and undergo continuous changes, requiring distributed adaptive control mechanisms that continuously account for the body's state. The canonical role of proprioception is representing the body state. We hypothesize that the proprioceptive system could also be critical for high-level tasks such as action recognition. To test this theory, we pursued a task-driven modeling approach, which allowed us to isolate the study of proprioception. We generated a large synthetic dataset of human arm trajectories tracing characters of the Latin alphabet in 3D space, together with muscle activities obtained from a musculoskeletal model and model-based muscle spindle activity. Next, we compared two classes of tasks: trajectory decoding and action recognition, which allowed us to train hierarchical models to decode either the position and velocity of the end-effector of one's posture or the character (action) identity from the spindle firing patterns. We found that artificial neural networks could robustly solve both tasks, and the networks'units show tuning properties similar to neurons in the primate somatosensory cortex and the brainstem. Remarkably, we found uniformly distributed directional selective units only with the action-recognition-trained models and not the trajectory-decoding-trained models. This suggests that proprioceptive encoding is additionally associated with higher-level functions such as action recognition and therefore provides new, experimentally testable hypotheses of how proprioception aids in adaptive motor control.
