**3. Methods**

In this section, we provide a full and detailed description of all the individual modules that compose the proposed mathematical model, each of which implements a different fragment of the whole synaptic transmission process.

This modular rationale at the base of the framework implementation guarantees an easy customization of the simulation pipeline, as well as the further extensibility of the system.

The current build of the framework includes:


The kinetic equations used to describe the reactions contained in both the RGS and CCS modules are implemented, exploiting the PySB python package [53] as systems of first-order differential equations. Numerical integration is performed using the SciPy ODE integrator [54]. All of the data analysis and fittings were performed using SciPy and Numpy packages [54,55]. Finally, all the plots were generated using the Matplotlib library [56].

All the code is stored in a publicly available github repository (https://github.com/ pietromicheli/CA3-CA1\_SynapticModel), where a jupyter notebook file for running simulations and performing basic analysis can also be found.

#### *3.1. Mathematical Model Implementation*

## 3.1.1. SPD Module

In this module, the stimulation pattern of the virtual synapse can be designed. Bidirectionality is a crucial feature of neuronal communication. The functional and topological properties of the brain neural network can be significantly shaped by the temporal relationship between forward and backward signals, as the STDP paradigm for the synaptic plasticity claims [20,21,57]. Therefore, integration of pre- and post-synaptic stimuli constitute a logic core of our implementation. For this purpose, patterns of pre- and post-synaptic stimuli can be programmed and simulated independently in order to analyze how the system behaves under different levels of synchronization between pre- and post-synaptic activities. Each pattern is modeled as a train of bursts. Numbers of stimuli per burst, intra-burst, and inter-bursts frequencies can be specified to design custom stimulation patterns (Figure 1).

In our model, pre-synaptic stimuli have been idealized as glutamate pulses, representing the instantaneous rise and fall in the free neurotransmitter concentration available inside the cleft compartment following pre-synaptic action potentials. Amplitude (i.e., the amount of available free glutamate) and width (i.e., the exposure time of the free glutamate) of the pre-synaptic glutamate pulses can be defined during the stimulation design.

On the other hand, post-synaptic stimuli have been modeled as transient depolarizations of the post-synaptic spine generated by dendritic back-propagating action potentials (bAP). Each bAP is shaped using a two-component exponential function, taken from the work by Shouval et al. [35]:

$$\text{bAP}(t) = V\_{\text{max}} \ast \left( \left( I\_{fast} \ast \exp\left(\frac{-t}{\tau\_{fast}}\right) \right) + \left( I\_{slow} \ast \exp\left(\frac{-t}{\tau\_{fast}}\right) \right) \right) \tag{1}$$

where *Vmax* is the maximum depolarization value for bAP value set to +67 mV [58], *Ifast* and *Islow* are the relative magnitudes of the fast and slow components of the bAP that sum to one, and <sup>τ</sup>*fast* and τ*slow* are the relative time constants that describe the exponential decays of the two components.

## 3.1.2. RGS Module

This module contains a system of kinetic equations describing the interactions between glutamate and AMPA/NMDA receptors, which take place inside the cleft compartment. The aim of this module is to accurately simulate both the receptors-neurotransmitter binding reactions and the gating mechanism that lead to opening or desensitization of the receptors.

Individual models describing the kinetic behavior of both AMPA and NMDA receptors have been selected from the literature based on their reproducibility and subsequently implemented as systems of first-order differential equations inside the PySB framework. To reproduce the kinetic behavior of AMPA receptors, we chose a model proposed by Koike et al. [36] for homomeric GluR2(flip) receptors. The model assumes two glutamate binding steps, one pre-open transient state, three desensitized states, and one open state of the receptor (Figure S2B). For the kinetic description of the gating mechanism of NMDA receptors, we used the model for GluN1/GluN2B NMDA receptor proposed by Amico-Ruvio and Pospescu [37]. This kinetic scheme includes two sequential glutamate binding steps, three pre-open transient states, two desensitization states, and one open state of the receptor (Figure S2A). Because we assume a saturating concentration of glycine inside the clef compartment, the binding steps with this molecule are not included in the kinetic model. Thus, all the resting NMDA receptors are considered glycine bound.

## 3.1.3. CPC Module

In this module, we implemented a set of equations with the aim of assessing the EPSCs and the respective EPSPs generated by the open fractions of both AMPA and NMDA receptors. EPSPs are then integrated with the back-propagating action potentials programmed during the stimulation design. Finally, the sum of all the depolarizing contributions is used to assess the variations of the post-synaptic membrane potential.

Many synaptic models that have been proposed in the past estimated the EPSCs and/or the EPSPs simply by using two-component exponential functions fitted on electrophysiological recordings [35,58–60]. On the contrary, in our model, the open probabilities of the receptors vary according to a system of deterministic rate equations that represent massaction kinetics of receptors-neurotransmitter interactions [53]. For this reason, the rising and decay phases of both receptor-mediated EPSCs and EPSPs responses are shaped by the complex receptor-specific interaction kinetics with the neurotransmitter. This confers more flexibility to our model, allowing us, for example, to explore the responses generated by mutant forms of the receptors by tuning the rate constants of some of the kinetic equations.

We defined the EPSCs of AMPA and NMDA receptors as follows:

$$EPSC\_{AMPA}(t) = O\_{AMPA}(t) \* G\_{AMPA} \* \left(V\_m(t - \Delta t) - V\_{E\_{AMPA}}\right) \tag{2}$$

$$EPSC\_{NMDA}(t) = O\_{NMDA}(t) \* G\_{NMDA} \* \left(V\_m(t - \Delta t) - V\_{E\_{NMDA}}\right) \* B(V\_m(t - \Delta t)) \tag{3}$$

where *OAMPA*(*t*) and *ONMDA*(*t*) are the number of open NMDA and AMPA receptors at each time step; *GNMDA* and *GAMPA* are the single channel conductance set to 40 pS and 15 pS, respectively [61–63]; *Vm*(*t* − Δ*t*) is the membrane potential at time (*t* − Δ*t*); *VE* is the channel-specific equilibrium reversal potential and defines the value of the membrane potential for which the electrochemical equilibrium is reached and, thus, the net flux through the channel is 0 (we assume that *VEAMPA* = *VENMDA* = 0) [64]; and *B*(*Vm*) describes the voltage dependence of the NMDA current given by the Mg2+ blocks defined by [35]:

$$B(V\_m) = \frac{1}{1 - \exp\left(-K\_M V\_m\right) \* \left(\frac{\left[M\_K^{2+}\right]}{3.27}\right)}\tag{4}$$

Once the EPSCs have been calculated, the relative EPSPs are determined simply by applying the law of Ohm:

$$EPSP\_{AMPA}(t) = EPSC\_{AMPA}(t) \* R\_{\\$} \tag{5}$$

$$EPS\_{NMDA}(t) = EPSC\_{NMDA}(t) \* R\_s \tag{6}$$

where *Rs* is the spine's resistance set to 500 M Ω [65].

Finally, the total membrane potential, defined as the sum of the partial depolarization contributions, is calculated according to the equation:

$$V\_{\rm mf}(t) = V\_{\rm r} + EPSP\_{\rm AMPA}(t) + EPSP\_{\rm NMDA}(t) + bAP(t) \tag{7}$$

where *Vr* is the resting membrane potential of the spine ( −65 mV).

In CA3 Schaffer collateral-CA1 synapses, the key mediator of the post-synaptic response is the elicited intracellular Ca2+ variation. Because NMDA receptors are the major

source of Ca2+ during spine stimulation [66], we explicitly calculate the NMDA receptormediated Ca2+ molar flowrate as follows:

$$I\_{\mathbb{C}a^{2+}}(t) = O\_{\text{NMD}A}(t) \* G\_{\mathbb{C}a^{2+}} \* \left(V\_M(t - \Delta t) - V\_{\mathbb{E}\_{\mathbb{C}a^{2+}}}\right) \* B\left(V\_m(t - \Delta t)\right) \tag{8}$$

where *GCa*2+ expresses the permeability of the NMDA receptor to Ca2+ ions, set to 2 nM·ms<sup>−</sup><sup>1</sup> ·mV−<sup>1</sup> [58] and *VECa*2+ is the reversal equilibrium potential for Ca2+ set to +130 mV [58].

Finally, the calcium dynamics in the post-synaptic cell is integrated by a simple first-order differential equation [35,58]:

$$\frac{d\left[\mathbb{C}a^{2+}(t)\right]}{dt} = I\_{\mathbb{C}a^{2+}}(t) - \frac{\left[\mathbb{C}a^{2+}(t)\right]}{\mathbb{T}\_{\mathbb{C}a^{2+}}} \tag{9}$$

where <sup>τ</sup>*Ca*2+ is the passive decay time constant of post-synaptic Ca2+ concentration, set to 20 ms [35].

A full list of all the parameters used in the equations described above is provided in the supplementary data (Table S2).

## 3.1.4. CAS Module

The last module of the pipeline contains a compartmentalized kinetic description of a reaction network that takes place inside the post-synaptic spine. Here, our rationale was to assess the variability in the amount of activated CaMKII enzyme upon different stimulation conditions. Because CaMKII plays a crucial role in the positive regulation of the early phase of LTP in CA3 Schaffer collateral-CA1 synapses [22–24], this estimation allows us to qualitatively infer the strength and the efficiency of the synaptic transmission. As previously described for the RGS module, we selected from the literature a kinetic model based on its reproducibility; we translated it inside the PySB framework and, finally, we appended the new block to the pipeline. For this purpose, we selected from the BioModels database [67] a model describing a set of interactions that, starting from post-synaptic rise in Ca2+ concentration, leads to the autophosphorylation (i.e., the activation) of monomeric CaMKII [38]. Particularly, the set of reactions implemented includes:


## *3.2. Data Fitting*

3.2.1. Concentration-Response Curves

We computed the glutamate concentration-response curves for NMDA receptors by stimulating the system with a glutamate pulse of 1.5 s in the absence of Mg2+ [47]. We run multiple simulations varying the amplitude of the glutamate pulse, with a concentration range between 0.01 and 1500 μM, and calculating the NMDA receptor-mediated current peak values. The EC50 values were then calculated by fitting the concentration-response data with the following equation:

$$Respons\,\%=\frac{100}{\left(1+\frac{EC\_{50}}{\left[\%Intamate/\text{s}\right]}\right)^{n}}\tag{10}$$

where *n* is the Hill slope.

#### 3.2.2. Two-Component Exponential Function Fitting

The deactivation time constant for the NMDA wild-type receptor and Glu413Gly and Cys461Phe variants were estimated as weighted time constants of the double exponential fit of the NMDA receptor current decay after the exposure of 1 mM glutamate for 1.5 s. The two-component exponential function used for the fitting takes the form:

$$I(t) = I\_{fast} \* \exp\left(\frac{-t}{\tau\_{slow}}\right) + I\_{slow} \* \exp\left(\frac{-t}{\tau\_{slow}}\right) \tag{11}$$

where *I* is the current, *Ifast* and *Islow* are the amplitudes of the fast and slow components, respectively, and <sup>τ</sup>*fast* and τ*slow* are the respective decay time constants. The weighted time constant of decay (<sup>τ</sup>*w*) was calculated using the following equation:

$$\pi\_{\rm IF} = \frac{I\_{fast}}{I\_{fast} + I\_{slow}} \ast \pi\_{fast} + \frac{I\_{slow}}{I\_{slow} + I\_{fast}} \ast \pi\_{slow} \tag{12}$$

#### 3.2.3. Four-Parameter Logistic Function and Bending Points

The data generated by the simulation of the relationship between different glutamate-NMDA Kd values and the concentration peaks of activated CaMKII enzyme (see results Section 2.4) were fitted with the four-parameter logistic function:

$$Y = \frac{a - d}{1 + \left(\frac{\chi}{c}\right)^b} + d \tag{13}$$

where *Y* represents the activated CaMKII response, *X* represents the affinity value Kd (expressed in μM), *a* is the lower asymptote, *d* is the upper asymptote, *c* represents the Kd that generates a mid-way response between the estimated *a* and *d*, and *b* is a slope factor. The bending point of the curve was then computed as follow:

$$X\_{bend} = \frac{a-d}{1+k} + d \tag{14}$$

$$\gamma\_{bend} = c \ast \left(\frac{a - \mathcal{Y}\_{bend}}{\mathcal{Y}\_{bend} - d}\right)^{\frac{1}{b}} \tag{15}$$

where *k* is a constant value, set to 4.6805 [68].
