A unified computational model for cortical post-synaptic plasticity

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 Ca²⁺/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.

Figure 1

Download asset Open asset

All pathways of Table 1B are Ca²⁺-dependent, but due to the variability in binding rates and quantities of different Ca²⁺-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 Ca²⁺ inputs and may serve as a basis for BCM-type rules of plasticity.

To examine the sensitivities of LTP- and LTD-inducing pathways to Ca²⁺, we simulated the injection of a prolonged square-pulse Ca²⁺ input of varying magnitude (illustrated in Figure 2A) into the post-synaptic spine and quantified the degree of activation of each of the Ca²⁺-binding molecules and the downstream signalling cascades. The simulations were carried out in the presence of mGluRs and $\beta$-adrenergic and cholinergic neuromodulation, which were modelled as prolonged square-pulse inputs as well.

Figure 2

Download asset Open asset

The injected Ca²⁺ quickly bound to Ca²⁺ 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 Ca²⁺ input. Consistent with experimental literature, a vast majority of Ca²⁺ 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 Ca²⁺ binding of these molecules in a steady state (5 min after the onset of Ca²⁺) and the overall activation/deactivation of downstream molecules as a function of the magnitude of the Ca²⁺ 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 Ca²⁺ flux is given, but their degrees of activity varied across the magnitude of the injected Ca²⁺ flux (Figure 2G). DGL was most completely activated throughout the Ca²⁺ amplitude, owing to the large equilibrium constant of its Ca²⁺ binding. At extremely large Ca²⁺ fluxes, CaM was more completely bound by Ca²⁺ than PLC and PLA2 (Figure 2G), but at lower Ca²⁺ 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 Ca²⁺ amplitudes. PKC, by contrast, became activated at relatively small Ca²⁺ 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 Ca²⁺ flux to the post-synaptic spine for the activation of the CaMKII pathway, relatively large Ca²⁺ flux for the activation of the PKA pathway, and relatively small Ca²⁺ flux for the activation of the PKC pathway.

The Ca²⁺ 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 ('Ca²⁺ 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 Ca²⁺ (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 Ca²⁺ entry (Figure 3D) caused a rapid increase in half-activated calmodulin (bound by two Ca²⁺ 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

Download asset Open asset

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 Ca²⁺ 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 Ca²⁺to the synapse model, which is a simple representation of inputs during synaptic plasticity induction protocols. Alternatively, Ca²⁺ 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 Ca²⁺ inputs and neuromodulatory inputs'), to determine the Ca²⁺ 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 Ca²⁺ 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, Ca²⁺, β-adrenergic and cholinergic ligands, and glutamate. All pathways leading to GluR1 and GluR2 phosphorylation and the consequent exocytosis and endocytosis are Ca²⁺-dependent: blocking Ca²⁺ 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

Download asset Open asset

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.

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 Ca²⁺ 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 Ca²⁺ 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 μm³ 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 Ca²⁺ flux, we could determine the number of Ca²⁺ 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 $N_{\mathrm{samp}}=200$ 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 Mg²⁺-gate opening in the NMDA receptor (Figure 5B–D, insets) (Hay and Segev, 2015; Markram et al., 2015). Thus, the Ca²⁺ flux time course varied across the pairing ISIs (Figure 5—figure supplement 1). These Ca²⁺ flux time series were imported into our biochemical model (Ca²⁺ 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 Ca²⁺ inputs) did not cause synaptic plasticity. Throughout the experiments, the activation of mGluRs was blocked.

Figure 5 with 1 supplement see all

Download asset Open asset

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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺-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 Ca²⁺ transient amplitude (Nevian and Sakmann, 2006), instead, our model suggests that the total Ca²⁺ is a better predictor of the plasticity outcome: the correlation coefficient between the post-STDP synaptic conductance and the peak Ca²⁺ transient amplitude (see Figure 5—figure supplement 1E) was 0.53, while that between the post-STDP synaptic conductance and the mean Ca²⁺ 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

Download asset Open asset

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 Ca²⁺ 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 'Ca²⁺ 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 Ca²⁺ 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., $\frac{{[\mathrm{GluR1}]}_{\mathrm{total}}}{{[\mathrm{GluR1}]}_{\mathrm{total}}+{[\mathrm{GluR2}]}_{\mathrm{total}}}$, from here on referred to as GluR1 ratio), the concentration of NCX (regulating the rate of Ca²⁺ 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 $f_{\mathrm{PKA}}$ and $f_{\mathrm{PKC}}$ 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 Ca²⁺ input.

We classified the parameter sets based on the total synaptic conductance 15 min after the onset of the stimulation with the high Ca²⁺ 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 Ca²⁺ injection resulted in LTD and the 250 particles/ms resulted in LTP.

Figure 7 with 1 supplement see all

Download asset Open asset

Our model with the standard protein concentrations (GluR1 ratio 0.5, [NCX] = 0.54 mM, $f_{\mathrm{PKA}}=f_{\mathrm{PKC}}=1.0$) 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 Ca²⁺ 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 $f_{\mathrm{PKA}}$, 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 (${\displaystyle ⪆0.6}$) GluR1 ratios, any NCX, PKA or PKC concentration could be used, but with smaller (${\displaystyle ⪅0.6}$) 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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ 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.

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 Ca²⁺ 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 Ca²⁺ dynamics, GluR1 activity, or the PKA-pathway signalling can have relatively large effects on the model output.

Figure 8

Download asset Open asset

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 Ca²⁺, β-adrenergic ligand and glutamate in addition to GluR1 fraction and factors for the protein concentrations of different pathways. Here, both the Ca²⁺ 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).

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

Download asset Open asset

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 ($f_{\mathrm{CaMKII}}$) 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 ($f_{\mathrm{PKA}}$) 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

Download asset Open asset

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 Ca²⁺ 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.

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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ entering the post-synaptic spine, and less so with the peak Ca²⁺ flux (Figure 6—figure supplement 1).The total amount of Ca²⁺ influx could thus provide a better biomarker for plasticity than the previously considered amplitude and duration of the Ca²⁺ transient (Evans and Blackwell, 2015).

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., Ca²⁺-induced Ca²⁺ release or cAMP-dependence of HCN channels) or not adequately strong in the multi-compartmental model. For example, the contributions of voltage-gated Ca²⁺ 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 Ca²⁺-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.

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

We created a model of pathways leading from Ca²⁺ 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 Ca²⁺ ions were linearly dependent on the number of Ca²⁺ ions. The response curve for CaM activation by Ca²⁺ 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 Ca²⁺ but milder activation in respose to small Ca²⁺ influx than the model of Jȩdrzejewska-Szmek et al., 2017 (Figure 11C).

Figure 11 with 2 supplements see all

Download asset Open asset

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 Ca²⁺ buffer with a Ca²⁺ 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.

Throughout the work, we simulated the signalling pathways in a single compartment representing a dendritic spine of size 0.5 μm³. 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 Ca²⁺ 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).

AMPARs 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 28⁴ 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

where N₁ and N₂ are the numbers of GluR1 and GluR2 subunits bound to the membrane, respectively, and $N_{1,\mathrm{phos}.}$ is the number of S831-phosphorylated GluR1 subunits at the membrane (note that the $N_{1,\mathrm{phos}.}$ subunits are included in all GluR1 subunits, that is, $N_{1,\mathrm{phos}.}\le N_1$). 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 N₁ GluR1 subunits and N₂ GluR2 subunits at the membrane can form is $\frac{N_1+N_2}{4}$. Here, we ignore the unpaired subunits by estimating that $\frac{N_1+N_2}{4}\approx \lfloor \frac{N_1+N_2}{4}\rfloor$ — 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:

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

We modelled the neurotransmission to the post-synaptic spine as fluxes of Ca²⁺ ions, β-adrenergic ligand, glutamate, and acetylcholine (labelled as Ca, L, Glu, and ACh, respectively, in Table 3). We used various stimulation paradigms: In sections 'Ca²⁺ 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 Ca²⁺ (3 ms), repeated at 100 Hz; 4xHFS — 4 trains of HFS, separated by 3 s of quiescence; LFS — 900 pulses of Ca²⁺ (3 ms), repeated at 5 Hz. Unless otherwise stated, each Ca²⁺ 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 Ca²⁺ inputs conducted by NMDARs when different stimulus patterns were applied. This model included the fast Na⁺ current $I_{\mathrm{Na},\mathrm{t}}$, Shaw-related K⁺ current $I_{\mathrm{Kv3}.1}$, muscarinic K⁺ current $I_{\mathrm{m}}$, and hyperpolarization-activated cyclic nucleotide-gated current $I_{\mathrm{HCN}}$ in the apical dendrite. The axo-somatic region contained all these (except for $I_{\mathrm{m}}$) as well as the low and high-voltage activated Ca²⁺ currents $I_{\mathrm{CaLVA}}$ and $I_{\mathrm{CaLVA}}$, the small-conductance Ca²⁺-dependent K⁺ current $I_{\mathrm{SK}}$, the transient K⁺ current $I_{\mathrm{K},\mathrm{t}}$, and the persistent Na⁺ and K⁺ currents $I_{\mathrm{Na},\mathrm{p}}$ and $I_{\mathrm{K},\mathrm{p}}$ 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 Ca²⁺ 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 Ca²⁺ 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 [Mg²⁺] 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 [Mg²⁺] 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 Ca²⁺ inputs were not considered in this work.

In 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 $f_{\mathrm{PKA}}\in [0,2]$, and, likewise, the concentrations of upstream PKC-pathway proteins mGluR, M1, Gq, and PLC using a factor parameter $f_{\mathrm{PKC}}\in [0,2]$. Furthermore, in section 'The model flexibly reproduces data from various cortical LTP/LTD experiments', CaM and CaMKII were altered in proportion by a factor $f_{\mathrm{CaMKII}}\in [0,2]$, phosphatases PP1 and PP2B by a factor $f_{\mathrm{PP}}\in [0,2]$, and phosphodiesterases PDE1 and PDE4 by a factor $f_{\mathrm{PDE}}\in [0,2]$. In both sections, the rapidity of Ca²⁺ 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 Ca²⁺ dynamics, we disregarded the data where free Ca²⁺ concentrations rose above 2 μM for one or more levels of Ca²⁺ 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 Ca²⁺ transients larger than 2 μM.

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

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

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

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