**1. Introduction**

Ionotropic glutamatergic receptors are a class of membrane receptors divided into three main subtypes, classified according to their activation to the selective agonists: NMDA (N-Methyl-D-aspartic acid), AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid), and Kainato. They play a key role in the process of synaptic transmission, which takes place in excitatory glutamatergic synapses, and dysregulations in their normal activities have been widely linked to numerous neurological disorders and synaptopathies [1–5]. Particularly, NMDA and AMPA receptors have been identified as crucial in the molecular mechanism underlying the process of synaptic plasticity, a process that leads to the modulation in the strength of the neuronal response to stimulation, linked to learning and memory [6–8].

Complex cognitive functions such as learning and multiple forms of memory are carried out by the hippocampal formation, which can dynamically sample, encode, store, and retrieve information coming from the sensory experience [9–11]. The constant encoding and integration of new information is possible thanks to the ability of a neural circuit to continuously reshape its topology and modulate the strength of its connections. In the hippocampal circuits, synaptic plasticity events that individual cells may undergo during synaptic transmissions occur in the form of Long Term Potentiation (LTP) and Long Term Depression (LTD). The trisynaptic circuit, particularly, has been extensively studied because

**Citation:** Micheli, P.; Ribeiro, R.; Giorgetti, A. A Mechanistic Model of NMDA and AMPA Receptor-Mediated Synaptic Transmission in Individual Hippocampal CA3-CA1 Synapses: A Computational Multiscale Approach. *Int. J. Mol. Sci.* **2021**, *22*, 1536. https://doi.org/10.3390/ ijms22041536

Academic Editor: Masoud Jelokhani-Niaraki

Received: 1 December 2020 Accepted: 1 February 2021 Published: 3 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of its apparently simple connectivity and the experimental accessibility of its structures. Inside this pathway, CA3 Shaffer collateral axons innervate CA1 pyramidal cells, forming excitatory glutamatergic synapses. The high density of NMDA receptors expressed on the surface of the dendritic CA1 spines confers to this synapse the ability to easily undergo NMDA receptor-mediated LTP and LTD, which has been substantially evidenced to be essential for some forms of explicit learning in mammals [12,13].

In Schaffer collateral-CA1 synapses, AMPA and NMDA receptors populate the membrane of the CA1 spine, actively participating in synaptic transmission. AMPA receptors are GluR1-GluR4 containing homo/hetero-tetrameric receptors that mediate fast excitatory neurotransmission in glutamatergic synapses. The early phase of synaptic plasticity events that occur in Schaffer collateral-CA1 synapses are associated with alterations in the number of AMPA receptors expressed on the spine membrane through activation of exocytosis or endocytosis mechanisms, as well as changes in AMPA receptors conductance through phosphorylation modifications [14,15]. Together, these molecular mechanisms lead to fine modulations in the strength of the synaptic transmission. The reactions underlying such modulation are controlled by the transient variations in the Ca2+ concentration that occur in the post-synaptic spine, especially, due to the activation of NMDA receptors. NMDA receptors are hetero-tetrameric glutamatergic ionotropic receptors permeable to Na2+, K+, Ca2+, and Mg2+ ions [16,17]. The permeability to Mg2+ ions gives to NMDA receptors a pronounced voltage-dependent behavior. At resting membrane potentials, external Mg2+ ions enter into the receptor's pore but, unlike the other permeating ions, they bind tightly to the pore, blocking it and impairing further ion permeation [18,19]. One of the most accepted physiological mechanisms needed to efficiently unblock NMDA receptors, thus generating an inward Ca2+ flux, is a temporal coincidence between the release of pre-synaptic neurotransmitter and a depolarization of the post-synaptic spine (of sufficient amplitude and duration) elicited by post-synaptic activity. This synchronicity is taken into account in the *Spike Timing Dependent Plasticity (STDP)* paradigm that also includes the post-synaptic dendritic activity expressed in the form of *back-propagating action potentials* (bAPs) [20,21]. The transient post-synaptic Ca2+ inward current generated by the activation and unblocking of NMDA receptors critically acts on the kinetic equilibrium of the different calcium-binding proteins involved in LTP/LTD-inducing pathways, such as Ca2+/Calmodulin-dependent Kinase II (CaMKII) [22–24].

Dysfunctions on LTP/LTD-mediated synaptic plasticity have been associated with many neurological disorders such as epilepsy and Alzheimers, Hughtington, and Parkinson's diseases [4,25–30].

A comprehensive and detailed understanding of the molecular mechanisms underlying synaptic transmission and neuroplasticity is then crucial for the physio-pathological characterization of many cognitive functions. However, even if LTP/LTD-mediated synaptic plasticity has been extensively studied, providing a substantial description of a full integration of the interaction networks underlying the whole synaptic transmission, deeply characterized at the molecular level, is currently a major challenge. This could be the starting point for the identification of new therapeutic strategies, aimed at re-tuning the global behavior of the intricate network of molecular interactions underlying synaptic plasticity, thus restoring its functional integrity.

Systems biology models have been shown to be key in approaching the complexity of this type of interaction networks. These models use a holistic approach to unveil the complexity of the molecular pathways and to catalogue all the biological complexes and the relationships between them [31]. They have evolved from empirical descriptions to fundamental mathematical equations applied by computational methods, allowing us to envision how such systems change over time under different conditions. In this way, one can infer qualitative features of the whole system, such as the downstream consequences of a single altered interaction, and consequently identify, for example, pharmacological targets or even predict the severity of a structural variant of a molecular species.

Here, we present, and render available to the scientific community (see Data Availability section), a mathematical model of the CA3 Schaffer collateral-CA1 transmission. Although other integrated and detailed models of glutamatergic synapses have been proposed recently [32,33], a clinical-oriented application of such models, able to also take into account the molecular characterization of particular disease-associated variants, is lacking. The rationale of our work was to provide a synaptic model that can be easily reproduced, run, and be integrated into larger analytical pipelines, proposing a novel viewpoint on the possible applications of comprehensive and detailed system biology models.

Our model allows us to simulate several features of the CA3-CA1 synaptic transmission process. These include (1) glutamate release inside the synaptic cleft as a result of a pre-synaptic stimulation, (2) bAP in the post-synaptic dendritic spine, (3) kinetic description of the gating mechanism of both NMDA and AMPA receptors, (4) estimation of the excitatory post-synaptic currents (EPSCs) and excitatory post-synaptic potentials (EPSPs), including the explicit calculation of the NMDA-mediated inward Ca2+ current, and (5) kinetic descriptions of the Ca2+-dependent molecular reactions that take place inside the post-synaptic spine and lead to the activation of CaMKII. Here we report some of the qualitative features observed in the receptors-specific contributions to synaptic transmission, as well as in the timing of pre/post-synaptic stimulation. Finally, we offer a further integration of our systems biology approach with a molecular level modeling of disease associated variants. This approach may pave the way to novel multiscale approaches to be used in the pharmacology or structural systems biology field. Because complex biological systems do not rely on individual metabolic networks, having a fully integrated description of metabolic networks allows us to envision the system as a whole instead of a sum of its parts [34]. It follows that the combination of integrated pathways with molecular detail observations, as the one we are presenting here, may bring to light new therapeutic strategies and bring us closer to the new era of personalized medicine.

#### **2. Results and Discussions**

This section is divided in two main subsections. In the first part, we present the implementation of the mechanistic model, providing an overview on the structure of the pipeline through the description of the individual modules, implemented to describe different fragments of the system. The second part contains the simulation of the model under different parameter configurations. This allows us to infer some qualitative features of the system, with a particular focus on the timing between pre and post-synaptic stimuli, and finally to assess shifts in the global system behavior given by the introduction of rare variants in the NMDA receptors associated with diseases.

#### *2.1. An Integrative, Python-Based Pipeline for Simulating Glutamatergic Synaptic Transmission*

We developed an integrative mathematical pipeline for the easy running of numerical simulations of synaptic transmission in individual CA3 Schaffer collateral-CA1 synapses, driven by both pre- and post-synaptic stimulation. The pipeline is composed of four different main modules, each one aimed at modelling a different part of the whole transmission process. Starting from the definition of a stimulation pattern, new modules were progressively implemented and added on top of each other, defining a linear pipeline for simulating the synaptic transmission (Figure 1).

**Figure 1.** Conceptual scheme of the pipeline to simulate our synaptic transmission model. This scheme illustrates all of the four modules of our framework: (i) Stimulation Pattern Design (SPD module); (ii) Receptors Gating Simulation (RGS module); (iii) Excitatory Post-Synaptic Currents and Potentials (EPSCs/EPSPs) Calculation (CPC module); (iv) Ca2+/Calmodulin-dependent Kinase II (CaMKII) Activation Simulation (CAS module).

> 2.1.1. Stimulation Pattern Design (SPD)

This module implements a series of functions that easily allow us to define a stimulation pattern that will drive the synaptic transmission. Such stimulation patterns can be composed of both pre- and post-synaptic stimuli, organized as trains of bursts. Here, highly customizable patterns can be designed by setting the number of stimuli composing each burst, the intra-burst, and the inter-burst frequencies for both pre- and post-synaptic stimuli (Figure 2).

**Figure 2.** Example scheme of a stimulation pattern. Pre- and post-synaptic stimuli are organized as trains or bursts. Each burst is composed of a sequence of stimuli, delivered at an intra-burst frequency. Inter-burst frequency defines the interval between each burst. Number of stimuli per burst, intra-burst, and inter-burst frequencies can be defined during the stimulation pattern design, for both pre- and post-synaptic patterns.

Pre-synaptic stimuli are idealized and modeled as the instantaneous rise and fall of the free glutamate concentration in the synaptic cleft, assuming a square pulse-like shape. In this article, we will refer to a pre-synaptic stimulus as a "glutamate pulse". The quantity of released glutamate (i.e., the pulse amplitude, expressed in μM) and the glutamate exposure time inside the cleft (i.e., the pulse width, expressed in ms) of each presynaptic stimulus can be independently parametrized. Post-synaptic stimuli are modeled as dendritic back-propagating action potentials, consisting of transient depolarization potentials of the post-synaptic spine membrane. The shape of such stimuli has been defined using a two-component exponential function (see methods, Section 3.1.1 for further details), as proposed by Shouval et al. [35]. The stimulation pattern defined in this module will constitute the input of the following modules.

#### 2.1.2. Receptors Gating Simulation (RGS)

Pre-synaptic stimuli, defined during the design of the stimulation pattern, are used as input to a second module, which is used to simulate the interactions between the neurotransmitter and the AMPA and NMDA receptors. This module contains the compartmental kinetic description of both the receptor-neurotransmitter binding reactions and the gating mechanisms that lead to the opening of the channels. Particularly, the latter consists of state-transition models (including closed, pre-open, open, and desensitized states) that statistically represent the stochastic distribution of the current traces recorded by electrophysiological experiments. We selected and integrated one kinetic model for both AMPA and NMDA receptors, proposed by Koike et al. and Amico-Ruvio and Popescu [36,37], respectively. Then, we translated both models into systems of first-order differential equations, implemented in a single larger kinetic model using the python PySB package (see methods for further details). Finally, a numerical integration was performed, allowing the simulation of the receptor's behavior with a high temporal resolution (integration step of 1 μs). We tested the reliability of these *ex-novo* implementations by comparing the behaviors predicted by our model, for both AMPA and NMDA receptors, with the behaviors reported in the works by Koike et al. and by Amico-Ruvio and Popescu [36,37] (Table S1). We observed a strong consistency between the kinetic features of both AMPA and NMDA receptors predicted by our PySB-based model and the respective original models, pointing to a high reliability of our implementation.

#### 2.1.3. EPSCs/EPSPs Calculation (CPC)

The third module of our framework consists of a system of equations used to explicitly calculate the EPSCs and the respective EPSPs generated during the simulation of the synaptic transmission. The EPSCs are estimated by calculating, over the simulation, the ion fluxes that permeate each open channel (predicted with the RGS module described in Section 2.1.2). This estimation is made according to the channel-specific conductance, the channel-specific reversal potential, and the depolarization level of the post-synaptic membrane. The EPSPs are then derived from the EPSCs (see methods Section 3.1.3 for further details). All the depolarization potentials, which include the EPSPs and, eventually, the bAPs arising from the post-synaptic stimulation, are summed together to assess the global changes in the membrane depolarization value. In this module, the equation for the explicit estimation of the NMDA-mediated Ca2+ current is used to assess the post-synaptic changes in the Ca2+ concentration according to a simple model proposed by Shouval et al. [35] (see methods, Section 3.1.3 for further details).

#### 2.1.4. CaMKII Activation Simulation (CAS)

The last module of our pipeline aims to simulate a kinetic description of the postsynaptic molecular interactions that controls the CaMKII kinase autophosphorylation events. For this purpose, as previously described for the RGS module (Section 2.1.2), we selected from the literature a detailed kinetic model based on its reproducibility, and we transcribed all its reactions into a second PySB model as a system of first-order differential equation. We chose to implement a model for the CaMKII activation proposed by Pepke et al. [38], and we integrated it into the simulation pipeline. This kinetic model includes a large number of reactions, mainly characterizing the interactions between free Ca2+ ions, calcium-binding messenger CaM, and the CaMKII enzyme. Particularly, the Ca2-CaM mediated autophosphorylation of CaMKII enzyme, which leads to its own activation, directly plays a pivotal role in inducing the early phase of synaptic plasticity [22–24]. Although the changes in the synaptic strength are currently not explicitly assessed in our model, the variations in the activated CaMKII accumulation allows one to assess the relative efficiency of the simulated synaptic transmission.

#### *2.2. Kinetic Behavior Analysis of AMPA and NMDA Receptors under Different Pre-Synaptic Stimulation Conditions*

We explored how AMPA and NMDA receptors kinetically behave under different stimulations patterns, exploiting the RGS module (Section 2.1.2). For this purpose, we simulated the model using different pre-synaptic stimulation patterns, consisting of either a single glutamate pulse or bursts of multiple glutamate pulses, delivered at different frequencies (ranging from 10 to 100 Hz). The amplitude of the glutamate pulses was set into a physiological range of 1–2 mM [39,40], while the time width was varied in a range between 1 ms and 1.5 s.

We first focused on the kinetic behavior of AMPA receptors under a single glutamate pulse of 1 mM, simulated with 1, 5, and 10 ms width. The desensitization kinetics of AMPA receptors predicted by the gating model shows a much slower time course (τ = ~25 ms, fitted with single exponential function) compared to the deactivation kinetics (τ = ~0.55 ms, fitted with single exponential function) after the end of a single glutamate pulse (Figure 3A). Moreover, both the exposure time of the glutamate (defined by the pulse width) and the stimulation frequency seem to strongly affect the number of desensitized receptors reached after a single pre-synaptic event [36] (Figure 3B). The faster deactivation, compared to the desensitization predicted by the model, points to the property of AMPA receptors to preferentially undergo a temporal accumulation of desensitized states instead of the open states.

We then analyzed how the variation of the glutamate pulses duration affects the summation of desensitized states under a single pre-synaptic burst stimulation. The latter was simulated by a single burst composed of 5 glutamate pulses of 1 mM amplitude and 1,5, and 10 ms width, with an intra-burst frequency of 100 Hz. We observed a significant increase in the temporal summation of desensitized AMPA receptors as the glutamate exposure values increased (Figure 3C–E, respectively).

Next, we analyzed the predicted kinetic behavior of NMDA receptors. By simulating a single glutamate pulse of 1 mM amplitude and 1 ms, 500 ms, and 1.5 s width, we observed a significatively slower deactivation and desensitization kinetics compared to AMPA receptors (Figure 4). Fitting the curves with a single exponential function, we found time constants of 163, 195, and 210 ms for the deactivation kinetics after 1 ms, 500 ms, and 1.5 s of glutamate exposure, respectively, and a time constant of 1.95 s for the desensitization kinetics (Figure 4A,B). From these results, we go<sup>t</sup> a ratio between the desensitization and the deactivation time constant (τdesens/τdeact) of ~12 for the NMDA receptors and ~45 for the AMPA receptors. The lower value found for the NMDA receptors leads to a more efficient temporal summation of its open states. In fact, when we simulated the model with a single pre-synaptic burst of 5 glutamate pulses of 1 mM amplitude and 1, 5, and 10 ms width, with intra-burst frequencies of 10, 50, and 100 Hz, we observed, effectively, summation of the open NMDAs (Figure 4C–E).

**Figure 3.** Opening and desensitization kinetics of the AMPA receptors. (**A**) Open fraction kinetics following a stimulation with a single glutamate pulse of 1 mM amplitude and width of 1 ms (black), 5 ms (dark grey) and 10 ms (light grey). Blue dotted trace shows the desensitization kinetics, while red dotted traces show the deactivation kinetics following glutamate removal from the synaptic cleft. (**B**) Desensitized fraction kinetics following a stimulation with a single glutamate pulse of 1 mM amplitude and width of 1 ms (black), 5 ms (dark grey), and 10 ms (light grey). (**C**–**E**) Kinetics of open and desensitized fractions following pre-synaptic stimulations with a burst composed of 5 glutamate pulses, with glutamate pulse amplitude of 1 mM, an intra-burst frequency of 100 Hz, and a pulses width of 1 ms (**C**), 5 ms (**D**), and 10 ms (**E**).

**Figure 4.** Opening and desensitization kinetics of the NMDA receptors. (**A**) Open fraction kinetics following a stimulation with a single glutamate pulse of 1mM amplitude and width of 1 ms (black), 500 ms (dark grey), and 1.5 s (light grey). Blue dotted trace shows the desensitization kinetics, while red dotted traces show the deactivation kinetics following glutamate removal from the virtual synaptic cleft. (**B**) Desensitized fraction kinetics following a stimulation with a single glutamate pulse of 1mM amplitude and width of 1 ms (light grey), 5 ms (dark grey), and 10 ms (black). (**C**–**E**) Kinetics of open and desensitized fractions following pre-synaptic stimulations with a burst composed of 5 glutamate pulses, with glutamate pulse amplitude of 1 mM, pulse width of 1 ms (black), 5 ms (dark grey), and 10 ms (light grey) and an intra-burst frequency of 10 Hz (**C**), 50 Hz (**D**), and 100 Hz (**E**).

To have a better insight into the difference between the kinetic behavior of AMPA and NMDA receptors, we simulated our model with a single pre-synaptic burst of 5 glutamate pulses of 1 mM amplitude and 1, 5, and 10 ms width, varying the intra-burst frequency between 10 and 100 Hz. For each intra-burst frequency, we calculated the ratio between the total number of desensitized and open receptors. According to our model, these simulations pointed out that the desensitized/open ratio of AMPA receptors depends more on the stimulation frequencies and on the glutamate pulses durations compared to the desensitized/open ratio of the NMDA receptors (Figure 5).

**Figure 5.** Desensitized/Open ratio expressed as a function of stimulation frequency. Simulations were performed using single pre-synaptic glutamate pulses of 1 mM amplitude and 1, 5, and 10 ms width. For each simulation, the ratio between the desensitized and the open fraction has been calculated for (**A**) AMPA (triangles, circles and diamonds symbols refer, respectively, to 1, 5 and 10 ms pulse width) and (**B**) NMDA receptors (crosses, circles and stars symbols refer, respectively, to 1, 5 and 10 ms pulse width).

#### *2.3. Temporal Relationship between Pre- and Post-Synaptic Stimuli Strongly Impacts Synaptic Transmission Efficiency*

During the stimulation of the synapse, the equations implemented in the CPC module (Section 2.1.3) allow us to explicitly assess the individual contribution of both AMPA and NMDA receptors to the global electrical transmission. Pre-synaptic-induced excitatory potentials and post-synaptic dendritic back-propagation events, programed during the stimulation pattern design, are integrated together to continuously estimate the variations in the NMDA permeability, as well as in the Ca2+ flux driving force (see methods Section 3.1.3 for further details). We explored, through several simulations, how the temporal relationship between pre- and post-synaptic stimuli can shape the efficiency of the electro-chemical transmission.

2.3.1. AMPA-Mediated EPSPs Are Not Sufficient to Efficiently Relieve the Mg2+ Block from NMDA Receptors

The pronounced voltage-dependent affinity of NMDA receptor for the extracellular Mg2+ ions causes the actual permeation of the channel to be strongly modulated by the depolarization level of the membrane [19]. We have previously observed that the kinetic equations implemented in the RGS module predict no effective temporal summations of open AMPA receptors because of their fast desensitization and deactivation kinetics, as observed in other studies [36,41]. Analyzing the output of the RGS module using the equations implemented in the CPC module (Sections 2.1.2 and 2.1.3), we then observed that, coherently, the AMPA-mediated responses also tend not to summate (Figure S1).

This observation prompted us to investigate if the amplitude of an AMPA-mediated EPSP evoked by a single pre-synaptic event was high enough to relieve the Mg2+ block from NMDA receptors. Because the EPSPs amplitudes of AMPA and NMDA receptors are

influenced by their levels of expression on the post-synaptic spine surface, we performed multiple simulations of a single glutamate pulse of 1 mM amplitude and 1, 5, and 10 ms width, varying the level of available AMPA receptors in a range between 20 and 200 [42]. Simulation results reported that the maximum AMPA-mediated EPSPs peaks elicited by single-pulse pre-synaptic stimulations reach −40 mV with 200 units of AMPA receptors (Figure 6A). According to the Mg2+ unblocking probability function that we have incorporated into the model (see methods Section 3.1.3 for further details), such a depolarization level can effectively release the Mg2+ ion from NMDA receptors only if the extracellular Mg2+ concentration is very low compared to the physiological concentration (Figure 6B), which is near to 1 mM [19].

**Figure 6.** (**A**) Simulated AMPA-mediated EPSPs evoked by different numbers of available AMPA receptors, ranging from 20 (lower trace) to 200 (upper trace). Solid, dashed, and dotted traces refer to single pulse stimulation performed with a glutamate pulse width of, respectively, 1, 5, and 10 ms. (**B**) Sigmoidal unblocking probability function for Mg2+ block, expressed as a function of membrane voltage. Each trace corresponds to a different value of extracellular Mg2+ concentration.

These results emphasize the fact that pre-synaptic events on their own may not be enough to ensure an effective Ca2+ permeation. As supported by the STDP paradigm, temporal coordination between pre- and post-synaptic events must occur in order to allow a significant Ca2+ influx that can effectively trigger plasticity [43].

2.3.2. Synchronization between Pre- and Post-Synaptic Stimulation Significantly Increases the NMDA Receptor Contribution to Synaptic Transmission

We further investigated how the synchrony between pre- and post-synaptic activity can affect the efficiency of synaptic transmission, particularly by increasing the amplitude of the NMDA receptors-mediated EPSCs and EPSPs. For this purpose, we compared the individual responses of the AMPA and NMDA receptors obtained from two different stimulation patterns, one including only pre-synaptic stimulation and one including coupled pre- and post-synaptic stimulations. In both stimulation pattern, the pre-synaptic stimulation consisted of a single theta burst composed of 5 glutamate pulses of 1 mM amplitude and 1, 5, and 10 ms width, with an intra-burst frequency of 100 Hz [44]; post-synaptic stimulation was designed as a single dendritic back-propagation event, which occurs in the post-synaptic spine 1 ms after the first pre-synaptic stimuli was delivered. Simulations were performed in the presence of 20 AMPA and 15 NMDA receptors [42,45], with extracellular Mg2+ concentration set to 1 mM. As expected, significant increases in the total NMDA receptor-mediated current peak (~2.5 fold), as well as in the Ca2+ that permeated the channel (~4.5 fold), were observed during the coupled pre- and post-synaptic stimulation compared to the pre-synaptic stimulation alone, showing the impact of bAP-mediated synaptic facilitation on the NMDA receptors conductance (Figure 7).

**Figure 7.** Simulation of synaptic transmission elicited by a single pre-synaptic burst of 5 glutamate pulses, in the absence of (**A**–**C**) or in the presence of (**D**–**F**) a single post-synaptic back-propagating action potential (bAP). (**A**,**D**) Time course of the individual AMPA-mediated EPSC. (**B**,**E**) Time course of the individual NMDA-mediated EPSC. (**C**,**F**) Time course of the Ca2+ molar flowrate that permeate NMDA receptors during the simulations. Pre-synaptic bursts were composed of 5 glutamate pulses of 1 mM amplitude and 1 ms (black pulses), 5 ms (dark grey pulses), and 10 ms (light grey pulses) width; in each plot, the responses elicited by 1 to 10 ms widths are represented by different colors, respectively, from the darkest to the brightest. Post-synaptic activity (red trace) was programmed as a single dendritic back-propagation event that occurs 1 ms after the first pulse of the pre-synaptic burst began. Both simulations were performed in the presence of 20 AMPA, 15 NMDA, and 1 mM of extracellular Mg2+.

Because we had observed that the presence of a bAP during stimulation significantly increases the NMDA receptor mediated EPSC, we analyzed how variations in temporal coordination level between pre- and post-synaptic stimuli impacts the amplitude of the elicited Ca2+ influx. For this purpose, we performed multiple simulations varying the time interval between pre- and post-synaptic stimuli (Δt=tpost − tpre). For each simulation, we then evaluated the effect of the bAP-induced synaptic facilitation by calculating the maximum Ca2+ concentration reached in the post-synaptic spine. Simulating a single pre-synaptic glutamate pulse of 1 mM amplitude and 1 ms width, together with a single post-synaptic bAP, we found that post-synaptic Ca2+ rises from a value of ~200 nM (the postsynaptic Ca2+ concentration elicited by a single pre-synaptic event alone) to a maximum of ~1.4 μM (Figure 8). This value is obtained when the pre-synaptic event precedes the post-synaptic event (positive Δt) of ~20 ms, in agreemen<sup>t</sup> with the Hebbian STDP paradigm for synaptic plasticity (see *Feldman 2012* [20] for a review).

**Figure 8.** Relationship between pre/post-synaptic stimulation timing and the Ca2+ concentration peaks reached in the post-synaptic spine. Simulations were performed in the presence of 20 AMPA, 15 NMDA, and 1 mM of extracellular Mg2+. Maximum post-synaptic Ca2+ concentration was reached with Δt ≈ 20 ms. Gray rectangle highlights negative Δt values in which post-synaptic stimuli precede pre-synaptic stimuli.

#### *2.4. Kinetic and Pharmacological Analysis of NMDA Variants: Multiscale Integration*

Deactivation time course defines the time required by the receptor-mediated current to decay after the removal of the agonist from the synaptic cleft. This kinetic feature, together with EC50 values of the agonist, constitute a prominent quantitative feature used to perform functional analysis of ion channels [46]. Many published studies on rare NMDA receptor variants have tried to assess the severity of a certain mutation, considering its impact on both glutamate potency and deactivation time constant [25,47–49].

We used our model to predict the glutamate affinity (Kd) and the weighted deactivation time constant (τw) in NMDA receptor variants, based on the EC50 values that have been reported in different experimental and computational studies [25,47,48,50]. In particular, we focused on two rare variants: Glu413Gly and Cys461Phe that fall inside the GluN2B binding pocket (Figure 9). These variants have been shown to decrease the glutamate potency, which may result from a decrease in the glutamate affinity [47,48,50,51].

**Figure 9.** Structure of human GluN1/GluN2A NMDA receptor (PDB accession code: 4TLM). The GluN2B subunit is colored in light blue. The insights show the glutamate binding domain of the wild-type (WT) receptor and the two structural variants Glu413Gly (E413G) and Cys461Phe (C461F). Each window focuses on the docked glutamate (white molecule) and the crucial residues that directly participate to the interaction. Red arrows point to the residue substitution of each of the two structural variants.

Therefore, we tuned the NMDA kinetic model to reproduce the same concentrationresponse behaviors experimentally observed for both the Glu413Gly and Cys461Phe variants.

Exploiting our kinetic model, we were able to computationally assess the NMDAglutamate concentration-response relationship by using the following approach: firstly, we sampled concentration values in a range between 0.01 and 1000 mM; next, for each value, we ran the RGS module, simulating a single glutamate pulse, with amplitude corresponding to the current glutamate concentration value and width of 1.5 s, as reported by experiments [47], setting the number of AMPA receptors to 0 (because we were interested in isolating the NMDA response). Finally, calculating from each simulation the peak of the evoked current, EC50 value was obtained by fitting the concentration-response data with the logistic function.

To predict the shifts in the NMDA receptor-glutamate affinity associated with the rare variants Glu413Gly and Cys461Phe, knowing their experimental EC50 values (75–79 μM for Glu413Gly [47,48] and 169 μM for Cys461Phe [47]), we progressively increase, during a sequence of multiple simulations, the ratio between the rate constants koff and kon (i.e., the Kd) of the equations describing the interaction between the NMDA receptor and the glutamate. For each simulation, we computed the EC50 value, and at the end of all the simulations we selected the Kd that rendered the EC50 values closest to the experimental ones.

As a result, we found that the NMDA receptor kinetic behavior generated by predicted Kd values shows a current deactivation time constantly very close to the experimental ones (Table 1).


**Table 1.** Predicted Kd and deactivation time constant for NMDA Wt and variants. Deactivation decay was fitted with a two-component exponential function, and the weighted Tau was then calculated (see methods Section 3.2.2).

The kinetic model of the NMDA receptor was tuned by only increasing the koff rate constant of the glutamate binding reactions. Therefore, we reasoned that the coherence between our results and the experimental data points to the fact that the analyzed variants are likely to affect the affinity of the receptor (thus causing an EC50 shifting) by negatively altering the glutamate residence time inside the binding pocket of the receptor.

For the wild-type NMDA receptor, we found a Kd value of 2.5 μM and a deactivation constant of 328 ms, whereas, for the Glu413Gly and Cys461Phe variants, we found Kd values of 190.5 and 446.5 μM and deactivation constants of 29 and 27 ms, respectively (Figure 10 and Table 1). As these results imply, the Glu413Gly and Cys461Phe variants increase the Kd of the glutamate ~75 and ~180-fold (Table 1).

**Figure 10.** Dose-response curves of the effect of glutamate on wild-type (WT) and variant-NMDA receptors. Simulation data was fitted with logistic regression. EC50 values of 2.7, 76, and 169 μM for WT, Glu413Gly, and Cys461Phe-NMDA receptors, respectively.

The next step in our multiscale analysis of NMDA Glu413Gly and Cys461Phe receptor variants consisted in further investigating if the calculated affinity alterations can impact the synaptic plasticity mechanism. To address this question, we simulated the effects of the structural variants on the amplitude of the rise in the post-synaptic Ca2+ concentration and on the amount of activated CaMKII, an enzyme that directly plays a pivotal role in triggering synaptic plasticity events in CA3-CA1 synapses. This latter estimation was done by exploiting the CAS module (Section 2.1.4). This module contains a mathematical description of the Ca2+-dependent CaM-CaMKII transduction pathway, which, starting from Ca2+ transients, leads to activation of CaMKII kinase (see methods Section 3.1.4 for further details). We stimulated our virtual synapse with a pair of single pre- and postsynaptic stimuli (glutamate pulse of 1 mM amplitude and 1ms width, time interval between pre- and post-synaptic stimuli of 20 ms). As expected, we found that the predicted decrease in the NMDA glutamate affinity significatively attenuates the amplitude of the elicited post-synaptic Ca2+ variation ~5 and ~8.5 fold for the Glu413Gly and Cys461Phe variants,

respectively (Figure 11A). Moreover, the kinetic model for the Ca2+-mediated activation of the CaMKII enzyme predicted much lower amounts of activated CaMKII for Glu413Gly (~13 fold) and Cys461Phe (~23 fold) variants compared to the wild type (Figure 11B). Considering the key role that the CaMKII enzyme plays in the molecular mechanism underlying the synaptic plasticity process, the predicted drastic decrease in the activation efficiency of such enzyme points out the severity of these rare structural variants. In fact, because CaMKII-driven neuroplasticity seems to be negatively affected in a significant way, severe neuropathological phenotypes, including learning and memory impairment, are likely to arise.

**Figure 11.** Variation of (**A**) Ca2+ concentration and (**B**) activated CaMKII over time for the WT and variant NMDA receptors. All simulations were performed under one pair of single pre- and post-synaptic pulses, with a pre-synaptic pulse of 1 ms of glutamate exposure, a delay between the pre- and the post-synaptic stimuli of 20 ms, and 1 mM of Mg2+.

In the last part of our in silico experiment, we were interested in reporting a more general representation of the relationship between NMDA-glutamate affinity and CaMKII enzyme activation efficiency. Here, our rationale was to search for an NMDA-glutamate affinity threshold that can be used for discriminating between high and low-severity variants, knowing their respective Kd. We proceeded, for this purpose, to simulate the whole synaptic model with the same basic stimulation pattern previously adopted for the analysis of the Glu413Gly and Cys461Phe variants but varying the Kd affinity value in a range between 1 and 1000 μM. For each simulation (i.e., for each Kd value), we selected the maximum amount of activated CaMKII observed. Data were first normalized to the maximum response observed across all the simulations and then fitted with the four-parameter logistic function (see methods Section 3.2.3 for further details) (Figure 12). Finally, the threshold was calculated by finding the bending point of the fitted curve, which corresponded to a Kd value of ~19 μM (Figure 12).

**Figure 12.** Variation of activated CaMKII as a function of the NMDA-glutamate Kd values. All simulations were performed under one pair of single pre- and post-synaptic pulses, with a presynaptic pulse of 1 ms of glutamate exposure, a delay between the pre- and the post-synaptic stimuli of 20 ms, and 1 mM of Mg2+.

The identification of this type of thresholds can be very useful for a rapid assessment of the downstream effects of variants and can be easily integrated into larger analytical pipelines. We are currently working on a further implementation of this synaptic model that also integrates a detailed kinetic description of the reactions controlling the phosphorylation of AMPA receptors by the CaMKII enzyme, an event which is known to directly control synaptic strength modulations (LTP and LTD) by altering the conductance and trafficking of these receptors. With this further extension, we aim to explicitly quantify synaptic plasticity events that can occur during the stimulations.
