Role of NMDAR plasticity in a computational model of synaptic memory

A largely unexplored question in neuronal plasticity is whether synapses are capable of encoding and learning the timing of synaptic inputs. We address this question in a computational model of synaptic input time difference learning (SITDL), where N‐methyl‐d‐aspartate receptor (NMDAR) isoform expression in silent synapses is affected by time differences between glutamate and voltage signals. We suggest that differences between NMDARs’ glutamate and voltage gate conductances induce modifications of the synapse’s NMDAR isoform population, consequently changing the timing of synaptic response. NMDAR expression at individual synapses can encode the precise time difference between signals. Thus, SITDL enables the learning and reconstruction of signals across multiple synapses of a single neuron. In addition to plausibly predicting the roles of NMDARs in synaptic plasticity, SITDL can be usefully applied in artificial neural network models.

. The SITDL hypothesis: An example of how activity-dependent changes in NMDAR population can affect synaptic current timing. (A) Left: A developing silent synapse starts with a majority of slow NMDARs with slower glutamate gate activation than fast NMDARs. The synapse receives a glutamate signal followed later by a dendritic voltage signal, causing a difference in activation of NMDARs' voltage and glutamate gates. Right: As the synapse develops and experiences the same timing between voltage and glutamate signals, its NMDAR population may change, caused by the difference in glutamate and voltage gate activation. In this case, fast NMDARs replace slow NMDARs until the gate conductance difference is minimized, thus aligning the overall NMDAR current peak with the peak of the voltage signal. This allows the developing synapse to learn the timing difference between the voltage and glutamate signals by encoding it in the NMDAR population's glutamate gate activation time. (B) Left: Each NMDAR has a glutamate gate and a voltage gate that can be activated independently, with corresponding gate conductances shown. Coincident activation of both gates allows influx of Ca 2+ (NMDAR Current). A slow NMDAR has a slower glutamate gate activation, therefore experiencing slower change in glutamate gate conductance than a fast NMDAR. Right: Circuit diagram of a post-synaptic spine containing only NMDARs. An NMDAR (R NMDAR ) is shown as variable serial resistances of its glutamate gate (R G ) and voltage gate (R V ) that depend on the glutamate signal (S Glu ) and dendritic spine's voltage (V), respectively. C M denotes the membrane capacitance, R M the membrane resistance, E Ca2+ the cell's Nernst potential for Ca 2+ , and I(t) the external input.

Scientific Reports
| (2021) 11:21182 | https://doi.org/10.1038/s41598-021-00516-y www.nature.com/scientificreports/ post-synaptic membrane and also couples them to intracellular signaling systems involving calmodulin 4 . The binding of the Ca 2+ -calmodulin complex to an NMDAR results in a Ca 2+ -dependent reduction of the NMDAR's channel opening frequency and channel open time. Furthermore, NMDAR activity regulates casein kinase 2 phosphorylation of GluN2B subunits, which disrupts the subunits' interactions with PSD-95 and decreases the GluN2B surface expression in neurons [9][10][11] . The wide diversity of ways in which NMDAR expression and dynamics can be modulated, through binding, post-translational modifications, and changing subunit composition, provides a wealth of opportunity for activity-dependent changes in synaptic current timing through NMDARs. While NMDARs are present even at "silent" synapses, which do not generate detectable excitatory postsynaptic potentials (EPSPs) in response to neurotransmitter release 12 , specific NMDAR subunit expression varies over time and across brain regions 7 . For instance, GluN2B expression tends to be highest in early development, while GluN2A subunits are more widely expressed in adulthood. Notably, there is evidence of rapid activitydependent bidirectional switching between GluN2A and GluN2B containing NMDARs, on the order of seconds in neonatal synapses 13 , and even in adult hippocampal synapses 14 . Distinct NMDAR expression is implicated not only in brain structures important for memory, but in other systems that rely on precise timing, such as in mammalian Calyx of Held and avian auditory brainstem [15][16][17] , important for inter-aural time and intensity difference calculations in sound localization, as well as in weakly electric fish relay cells, which are important for precise temporal regulation of the jamming avoidance response 18,19 . Notably, these delay line systems can resolve temporal disparities in the microsecond range, which may arise at the level of individual neurons. For instance, single neurons of the pre-pacemaker nucleus in a weakly electric fish are sensitive to temporal disparities as small as 1 µs, and their signaling and precision may involve NMDARs 20,21 . In one of the few studies on changes in NMDAR subunit composition during the development of a delay line system, it has been shown that GluN2A replace GluN2B subunits in chicken cochlear nuclei, after onset of hearing 17 .

SITDL hypothesis.
Much is still unknown about NMDARs for the intricacies of their dynamics and regulation, the roles of different subunits, and their specific roles in memory formation. Here, we propose a hypothetical model, Synaptic Input Time Difference Learning (SITDL), of how activity-dependent changes in fast and slow NMDAR expression could affect synaptic current timing. In particular, we assume that synaptic NMDAR expression depends on differences in voltage and glutamate gate conductances, and that matching of these conductances provides optimal Ca 2+ influx through NMDAR gates. This process is driven by changes in subunit composition of the NMDARs, which determine glutamate gate characteristic activation time.
For example, the left panel of Fig. 1A shows a silent synapse having a majority of slow NMDARs, with slower activation times in response to glutamate than fast NMDARs. The slow NMDARs can resemble NMDARs containing GluN2B subunits, while fast NMDARs resemble those with GluN2A subunits. The silent synapse receives a synaptic glutamate signal followed by a voltage signal carried by the dendrite from a non-silent synapse, whose peaks occur several milliseconds apart. This difference in timing accompanies a difference in NMDAR voltage and glutamate gate conductance. We propose that as the synapse develops and experiences the same timing difference between voltage and glutamate signals, the difference in glutamate and voltage gate conductances causes its NMDAR population to change.
In the right panel of Fig. 1A, fast NMDARs replace slow NMDARs until there is a minimization of the difference in glutamate and voltage gate conductance, thus aligning the NMDAR glutamate gate conductance peak with the peak of the voltage signal. In effect, the synapse learns the timing difference between the voltage and glutamate signals, encoding it in the NMDAR population's glutamate gate activation time, which is a function of the numbers of fast and slow NMDARs. Then, stabilization of a synaptic NMDAR population, such that there are no more removals and insertions of fast or slow NMDARs, possibly due to PSD-95 anchoring, would mean a stable memory of the timing difference. Note that silent synapses are particularly useful substrate for this SITDL mechanism, as glutamate signals alone do not cause depolarization of the synapse, thus providing a greater degree of independence between NMDAR glutamate gate and voltage gate activation.
Furthermore, maturation, which involves the insertion of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic receptors (AMPARs), may provide the synapse with a way of "recalling" the memorized timing difference. In this case, the NMDAR population no longer must rely on passing dendritic voltage signals for activation of their voltage gates, as AMPARs can provide depolarization in response to the glutamate signal. Thus, in theory, a glutamate signal alone could provide activation of NMDARs' glutamate gates, and also indirectly provide depolarization to activate NMDARs' voltage gates through the AMPARs. The resulting overall NMDAR current depends on the glutamate gate activation time of the stabilized NMDAR population, which encodes the previously memorized timing difference, in effect, recalling it.
In particular, SITDL and recall mechanisms may be useful for signal reconstruction when there are many silent redundant synapses and few non-silent synapses, as there often are in developing neurons 22 . Notably, networks of developing neurons have been shown to exhibit coordinated activity and rhythmic firing patterns 23,24 , essential to normal development of networks and synaptic connections. These developmental activity patterns are significantly different from those in mature circuits 25 . Thus, if a developing neuron receives rhythmic glutamate signals from presynaptic neurons, activation of a non-silent synapse would produce a rhythmic dendritic voltage signal in the neuron. Each of its silent synapses could therefore potentially receive both rhythmic voltage and glutamate signals, which would produce consistent overlap between signals. If these synapses could learn and encode the timing differences of these rhythmic signals during development, then once they mature, they would be able to recall and reproduce certain parts or patterns from the original developmental signals.
We explore the hypothesis of SITDL further and provide a computational single-compartment model (Methods). We also show that with many redundant silent synapses, as in a developing neuron, SITDL mechanisms, along with activity-dependent synaptic elimination and maturation, can enable reconstruction and recall of www.nature.com/scientificreports/ original developmental signals. Supplemental Slides S1 provide an overview and visualization of the SITDL mechanism and multi-synaptic signal reconstruction.

Methods
Input signals. To construct dendritic and glutamate signals to test the plasticity model, periodic input signals, S N and S Glu ( Figure S2), were generated using a gaussian-like shape for the spikes in voltage, and a rightskewed shape for synaptic glutamate concentration 26-28 : where S N = S N (t) determines the shape of the voltage input spike with parameters t N0 and σ N. Glutamate signal, S Glu = S Glu (t) has characteristic right-skewed shape, with constant α L that determines the width and decay of the spike, and t L0 is simply shifted relative to t N0 , with t L0 = t N0 − 1/α L , to ensure that the spike peaks of both S N (t) and S Glu (t) occur at the same exact time ( Figure S2). The input signals, S N and S Glu were normalized by their respective maximums, S NMax and S GluMax such that S N = S N /S NMax and S Glu = S Glu /S GluMax . The dendritic signal, I D = I D (t), depends on I 0 , a constant current, and the signal S N . Values for all constants are given in Table 1.
NMDAR voltage and glutamate gate dynamics of a single synapse. NMDARs play important role in learning, memory, and development by regulating synaptic properties at the post-synaptic membrane. Each NMDAR is a heteromeric cation channel with a glutamate-activated gate and voltage-activated gate. Activations of the glutamate and voltage gates are largely independent, and either alone does not open the NMDAR channel completely. However, the activation of either gate will cause specific changes in conformation or channel structure and may potentially affect the NMDAR's binding and phosphorylation sites. With near-coincident activation of both gates, the NMDAR channel opens, permitting Ca 2+ influx. We consider the glutamate and voltage gates as separate variable resistors, with each NMDAR composed of the glutamate gate resistor in series with the voltage gate resistor (Fig. 1B, Right). Therefore, in a single synapse, the overall post-synaptic NMDAR conductance, g = g(t), can be expressed as follows: www.nature.com/scientificreports/ where g Glu = g Glu (t) represents the conductance of the NMDAR population's glutamate gates, and g V = g V (t) represents the conductance of its voltage gates. g determines the level of permeability and Ca 2+ influx through channels of the NMDAR population ( Fig. 1B, Left). The product of g Glu and g V in Eq. 4 also implies necessity of coincidence of glutamate and voltage gate activation. The voltage dependence of an NMDAR has been previously modeled with experimental data as a logistic function of synaptic Mg 2+ concentration and the post-synaptic voltage 29,30 . We can similarly describe the conductance of a population of NMDARs' voltage gates, g V , as a simplified logistic function of the normalized voltage along the dendrite, V D = V D (t): with constants a v and b v provided in Table 1. In Eq. 5, we assume that Mg 2+ concentration stays relatively constant in the synapse, as expected for physiological conditions. The evolution of g Glu can be written as: with approximate solution in the form of a single-time-step mapping: where τ Glu is glutamate gate conductance characteristic time, Δt is a single time step, and a L and b L are constants between 0 and 1, the values for which are provided in Table 1. Here, we assume that glutamate gate conductance limit, g L , acts as an eligibility trace of glutamate signal S Glu (Eq. 8), since NMDAR conductance depends on synaptic glutamate concentration and prior activations. On a sub-second time scale, local changes in τ Glu represent modification of receptors that transiently affect their dynamics and membrane stability, such as phosphorylation of NMDARs 9 . On a longer scale of seconds or more, overall changes in τ Glu define the changes in numbers of slow and fast NMDARs, which significantly affect the dynamics of NMDAR glutamate gate conductance, g Glu . A key point and assumption of the SITDL hypothesis is that NMDAR subunit composition, and consequent value of τ Glu , depend on the difference between the NMDAR gate conductances, g Glu and g V . So, the optimal Ca 2+ influx occurs when the value of g Glu matches that of g V . Following the assumption, temporal evolution of τ Glu can be written as: where Δτ is time step for τ Glu , γ is a scaling constant. The primary goal of the SITDL model is to find value of τ Glu that minimizes the gate conductance mismatch, (g Glu − g V ), which would effectively represent the estimated time difference between the synaptic glutamate signal and the dendritic voltage signal. Equation 9 results from gradient descent minimization of the function F(τ Glu ) = (g Glu − g V ) 2 , with dF/ dτ Glu =− 2t(g Glu − g V )(g L − g Glu )/ τ Glu 2 . Equation 9 has two fixed points, at g Glu = g V and g L = g Glu . Figure S4 shows evolution of (g L − g Glu ) and (g Glu − g V ) toward stable solution, where the peaks of g Glu and g V are aligned. Optimal solution is reached when g Glu → g V and g L ≈ g V .
There are few studies modeling the dynamics of the NMDAR glutamate gate and how they change with different NMDAR subunit compositions. For instance, they primarily involve kinetic modeling of NMDARs with different subunit compositions [31][32][33] , but provide little suggestion for modeling the synaptic population of NMDARs, its overall dynamics, and how the NMDAR population may change over time in an activity-dependent manner. While τ Glu describes the dynamics of NMDAR glutamate gate conductance, it depends on the glutamate signal and is difficult to compare to the experimentally known characteristics of slow and fast NMDARs. Therefore, it is useful to look at the NMDAR dynamics in response to a fixed signal: a single glutamate spike. For each synapse, we are particularly interested in the g Glu rise-to-peak time, τ Syn , in response to a single glutamate spike, which is defined through the numbers of slow and fast NMDARs (n Slow and n Fast ) 34 : where n Total is the total number of NMDAR receptors, τ Fast and τ Slow are known time constants of fast and slow NMDARs, respectively.
Using Eq. 14, and the assumption that n Total stays constant, numbers of slow and fast NMDARs can be calculated from τ Syn as follows: www.nature.com/scientificreports/ τ Syn is calculated as the rise-to-peak time of NMDAR glutamate gate conductance, g Glu , for a single glutamate spike as signal S Glu , similarly to previous studies 34 . τ Syn dependence of τ Glu is shown in Figure S3. τ Syn values were calculated from g Glu response to a single glutamate spike, at fixed τ Glu vaues (ranging from τ Glu = 1 ms to 2000 ms).
Synaptic current and dendritic voltage. We can calculate the synaptic current, I Syn = I Syn (t), that results from NMDAR activation and consequent Ca 2+ influx, using overall NMDAR conductance, g, and dendritic voltage, V D : Changes in dendritic voltage are described via a single-compartment conductance-based model: where V Rest is resting membrane potential, τ R = CR is a time constant, k D and k S are current contribution constants, τ D is the dendritic delay time constant, indicating the time difference between synaptic and dendritic signal, I D (t − τ D ) is delayed dendritic signal, and I Syn is synaptic current resulting from NMDAR activation.
Stabilization of glutamate gate conductance characteristic time τ Glu . The change in τ Glu described by Eq. 9 is always dominated by the gate conductance mismatch. Therefore, unless there is a consistent and perfect match of the gate conductances, g Glu and g V , τ Glu will constantly change. We propose the existence of a stabilization mechanism that depends on the overall NMDAR conductance, g, such that when g is consistently large, τ Glu will eventually stabilize. We modify Eq. 9 to include a stabilization variable, P, with the following set of equations: where P and σ define the stabilization mechanism: when g is consistently greater than zero (Eqs. 17,18), and concurrently the change in τ Glu , Δτ Glu , is close to zero, τ Glu stops changing. Δτ Max is the maximum possible value for Δτ Glu . Its value and the values for constants a P , b P , and τ P are provided in Table 1.

SITDL multi-synaptic memory and recall.
To test whether SITDL mechanisms could be used to memorize peak times of a glutamate signal and reconstruct it, several alterations are made. Multiple synapses are used, with each synapse represented as a single-compartment SITDL model with stabilization mechanisms (Eqs. 1-18), with its own dendritic delay time constant, τ D , and initial glutamate gate characteristic time, τ Glu . A set of "Learning Phase" simulations starts with a full set of synapses (i = 1, …,N), with 49 τ D values uniformly distributed between 4 and 100 ms, and initial τ Glu = 20 ms. In these simulations, each synapse receives a periodic glutamate signal as described before, and a sparse voltage signal with one spike repeating over the same period as the glutamate signal (Figs. 5, S2 bottom). For each synapse i, g Avg , the overall NMDA conductance averaged over the last 0.1% of the simulation time (last 150,000 points), is calculated. If a synapse's g Avg is lower than a threshold factor, δ, times the value of g Avg averaged across the entire set of synapses, such that g i Avg < δ N N j=1 g j Avg , then it is eliminated. This is similar to the performance-based synaptic elimination rule used in the Clusteron model 35 , and to minimal-value deletion algorithms 36 . All synapses that have not been eliminated are considered stabilized, with τ Glu permanently fixed. We assume that stabilized synapses express AMPARs, which depolarize the synapse upon binding glutamate. Therefore, for "Recall Phase" simulations, each stabilized synapse receives a single glutamate spike coincident with a single voltage spike (Fig. 5C). Note that k D and k S , the constants for dendritic and synaptic current contribution to voltage, respectively, are set to 2.0 and 3.0 for these simulations. The resulting synaptic voltage traces are shifted by their corresponding τ D and summed over all stabilized synapses to give a "recall signal", which is compared against the original glutamate signal used during Learning Phase simulations. A smaller set of 12 synapses was also used for SITDL learning and recall simulations, with τ D values uniformly distributed between 8 and 96 ms, and initial τ Glu = 20 ms. "No-SITDL" simulations use the same Learning Phase and Recall Phase simulations, except there are no SITDL mechanisms during the Learning Phase, such that there are no changes in τ Glu for any synapse.

Mutual information (MI) analysis.
A mutual information (MI) estimator (Adaptive partition using Interspike intervals MI Estimator (AIMIE)), was used to compare recall signal and original glutamate signal 37 . Short signals were repeated with proper periodicity, such that each signal had at least 4000 spikes for more accurate (13)  www.nature.com/scientificreports/ MI estimation. A higher MI estimate for the two signals suggests that they have a greater dependency and higher degree of similarity. AIMIE has been shown to work well with spike time series with disparate firing rates and may provide more accurate estimates of MI than several other commonly used MI estimators. Details on AIM-IE's calculations and its use in information flow analysis of a spiking network are provided in previous work 37 .
All simulations of the SITDL model were programmed and run in Spyder 3.0.0, a Python 3.5 environment. SITDL data were analyzed using MATLAB R2016b and Excel 2016.

Results
The SITDL mechanism increases overlap of NMDAR gate conductances. The core mechanism of the SITDL model is learning the timing difference between synaptic glutamate and dendritic voltage signals. This involves changes in numbers of fast and slow NMDARs, and consequent changes in NMDAR glutamate gate characteristic time, τ Glu , and τ Syn , the corresponding rise-to-peak time of the glutamate gate conductance in response to a single glutamate spike. Using periodic glutamate and voltage signals and computational singlecompartment model described in Methods (Eqs. 1-14), we ran multiple simulations for a range of initial values of dendritic delay, τ D , and τ Glu . It may be noted that the relative timing delay between signals, τ D , can result from travelling along post-synaptic dendrites, as well as from travelling along the pre-synaptic axon. Conduction delays in certain systems, such as corticothalamic projections and some unmyelinated axons, can range from less than 2 ms to as long as 100 ms [38][39][40] , and in development the regulation of myelination can also significantly affect conduction delays 41 . Figure 2 shows the evolution of τ Glu and NMDAR conductances for a single synapse receiving the glutamate signal, S Glu , followed by the dendritic voltage signal, V D . Both are copies of a similar periodic signal, with dendritic signal delayed by time τ D . The timing difference between signals produces a gate conductance mismatch, (g Glu − g V ) (Eq. 9), changing τ Glu and shifting g Glu peaks. This leads to greater overlap of the gate conductances.
With initial values τ Glu = 5 ms (corresponds to τ Syn ≈ 7 ms) and τ D = 15 ms ( Fig. 2A), glutamate gate conductance, g Glu , initially peaks before voltage gate conductance, g V . This results in overall increase of τ Glu over time. Coincidence appears to have been achieved at the end of the 4000 ms simulation. In Fig. 2B, g V initially peaks before g Glu , and τ Glu decreases over time, resulting in greater overlap of gate conductances. Figure 2C shows that SITDL mechanisms can still work to achieve greater gate conductance coincidence even when dendritic delay, τ D , is significant (τ D = 95 ms) and there is much less overlap in signals. Figure 3A-C shows simulations of the SITDL model with a greater range of initial delays and longer simulation time of 600,000 ms to explore the evolution of τ Glu . For initial τ Glu = 5 ms, and τ D = 1 ms to 17 ms, τ Glu seems to quickly attain relatively stable values, with little change towards the end of the simulation. The smaller the τ D , the more quickly τ Glu stabilizes. This can be explained by smaller gate conductance mismatch and greater overlap between glutamate and voltage signals. For τ D = 18 to 45 ms (Fig. 3B), τ Glu grows but does not stabilize as much as for smaller τ D . For some τ D , τ Glu keeps growing past 1000 ms without achieving greater stability, suggesting an insufficient overlap of signals. Likewise, for τ D = 46 to 100 ms (Fig. 3C), only certain τ D values, such as τ D = 49 to 60 ms, 70 to 75 ms, and 90 to 100 ms, seem to provide sufficient overlap for τ Glu to approach small stable values at the end of the 600,000 ms simulation time. Different initial τ Glu and longer simulation times may allow τ Glu to stabilize, paticularly for τ D values larger than 17 ms, and to potentially achieve coincidence of glutamate and voltage gate conductances.
It may be noted that while τ Glu approaches a more stable range of values in most cases, τ Glu does not fully stabilize in Fig. 3A-C. Thus, we made key alterations to the SITDL model, using Eqs. 15-18 for implementing τ Glu stabilization mechanisms, which are important for demonstrating stable memory formation and recall.  Fig. 3G-I shows that τ Glu can decrease from a higher initial value (τ Glu = 50 ms), and fully stabilize at values similar to those in Fig. 3D-F, due to SITDL stabilization mechanism. In some cases, the higher initial τ Glu allows for stabilization at additional τ D values (Fig. 3H). While we showed that τ Glu is able to stabilize at values lower than 5 ms (Fig. 3A,D,G) to demonstrate the range of capabilities of the SITDL algorithm, for all subsequent simulations we limited τ Glu to values between 5 and 1410 ms to correspond to the more biologically relevant NMDAR conductance rise-to-peak time limits of around 7 ms and 50 ms, respectively. NMDAR τ Glu corresponds to fractions of slow and fast NMDARs. We assumed that τ Glu is determined by the numbers of fast and slow NMDARs. Fast NMDARs can be considered those with GluN1/GluN2A subunit composition, with estimated rise-to-peak time, τ Syn , of 7 ms, and slow NMDARs as those with GluN1/ GluN2B subunit composition, with estimated rise-to-peak time τ Syn of 50 ms 8 . Over the course of development, slow NMDARs are typically replaced with fast NMDARs 7 . As a test, we ran a SITDL model simulation with initial values τ Glu = 150 ms and τ D = 10 ms, and calculated numbers of fast and slow NMDARs using estimated τ Syn and Eqs. 10-12 with n Total = 50, τ Fast = 7 ms, and τ Slow = 50 ms. τ Syn was estimated from τ Glu values, using linear interpolation of data points from the τ Syn vs. τ Glu plot of Figure S3.
Because τ Glu decreases, the g Glu maximum shifts significantly closer to g V maximum (Fig. 4A). Initially, τ Glu = 150 ms (τ Syn ≈ 30 ms), with 26 slow NMDARs and 24 fast NMDARs. As τ Glu decreases and stabilizes at its final value of ~ 12.7 ms (τ Syn ≈ 11 ms), there are 5 slow NMDARs and 45 fast NMDARs (Fig. 4B,C). Reduction of τ Glu can therefore represent a replacement of slow NMDARs with fast NMDARs, and τ Glu convergence can be considered as stabilization in expression of fast and slow NMDARs.

SITDL mechanisms across multiple synapses enable reconstruction of synaptic signal.
The potential for SITDL mechanisms to achieve greater coincidence between NMDAR glutamate and voltage gate conductances, even with very limited overlap of signals, suggests that they could be used to memorize peak times of a glutamate signal and reconstruct it. To explore this, large numbers of SITDL synapses, with stabilization mechanisms and τ D values uniformly distributed between 4 and 100 ms, were run through Learning Phase simulations like the previous simulations in Fig. 3, with ranges of different initial conditions. During the Learning Phase simulation, each synapse receives a periodic glutamate signal, and a sparse rhythmic voltage signal of the same period, delayed relative to the glutamate signal by time constant, τ D (Fig. 5A,C). We can consider this initial set of synapses to be located along a neuron's dendrite, receiving the same glutamate signal from the branched axon of a single pre-synaptic neuron (Fig. 5A). These types of redundant multi-synaptic connections are seen in neocortex, striatum and hippocampus [42][43][44][45] , and they may also arise as a result of long-term potentiation (LTP) 46 . Notably, this synaptic multiplicity appears to increase during development in hippocampal CA3-CA1 synapses 47,48 . As before, τ Glu and conductance values change over the course of the Learning Phase simulations (Fig. 5B,C). After simulations end, synapses with insufficient overall NMDAR conductance are eliminated, similar to NMDAR-dependent and Ca 2+ -based synaptic elimination in biological neurons 49,50 , and to synaptic elimination In each case, the synapse receives a synaptic glutamate signal (light green) followed by a dendritic voltage signal (light blue) delayed by time τ D . τ Glu changes over time due to the NMDAR gate conductance mismatch. (A) Simulation with initial values τ Glu = 5 ms and τ D = 15 ms. Initially, glutamate gate conductance, g Glu , peaks before voltage gate conductance, g V . Over time, τ Glu grows, delaying the g Glu peak and achieving greater coincidence. (B) Simulation with initial values τ Glu = 50 ms and τ D = 10 ms. In this case, g V initially peaks before g Glu . τ Glu decreases over time, resulting in greater coincidence of the gate conductances. (C) Simulation with initial values τ Glu = 50 ms and τ D = 95 ms. Though there is much less overlap between glutamate and voltage signals, due to periodicity of the signal, the SITDL mechanism still appears to achieve greater coincidence over time.  36 . Figure 5B shows the changes in τ Glu for synapses that were stabilized, with an initial set of synapses having initial conditions τ Glu = 20 ms, and τ D uniformly distributed from 4 to 100 ms, with a step size of 2 ms. The τ D step can represent how far apart the initial synapses are on the dendrite (Fig. 5A). The leftover stabilized synapses are assumed to express AMPARs. These were run through Recall Phase simulations, where each stabilized synapse received a single glutamate spike and a single voltage spike that coincided due to AMPAR expression (Fig. 5D). The summation of resulting Recall Phase voltage signals across the stabilized synapses (Fig. 6, traces in color), shifted by the corresponding τ D value, provides a "recall signal" that very closely resembles the original glutamate signal from the Learning Phase (Fig. 6A). For recall signals constructed from SITDL simulations with much larger τ D step, there are less spikes and peak times are slightly less similar to those of the glutamate signal (Fig. 6B). Furthermore, recall signals constructed from simulation sets (Learning Phase, elimination, Recall Phase) without SITDL mechanisms such that τ Glu cannot change, appear to be even less similar to the original glutamate signal, as they are missing the first, second, and fifth glutamate spikes, and have additional spikes in the vicinity of the third and fourth glutamate spikes (Fig. 6C,D). Using the MI estimator AIMIE 37 , we show that there is a greater degree of similarity between the original signal and recall signal reconstructed from SITDL simulations (Fig. 6E). SITDL mechanisms across a set of synapses, particularly when they are distributed closely enough, can therefore be used to effectively memorize and recall synaptic signals. In most cases, τ Glu increases, but there does not appear to be sufficient simulation time for it to reach more stable values as in (A). (C) Changes in τ Glu with initial conditions τ Glu = 5 ms and τ D value ranging from 46 to 100 ms, which provide much less overlap of glutamate and voltage signals. Certain τ D values, such as 49 to 60 ms, 70 to 75 ms, and 90 to 100 ms, seem to provide sufficient overlap for τ Glu to change significantly and potentially achieve coincidence of gate conductances with enough simulation time. (D-F) Same simulations and initial conditions as in Fig. 4A-C, but with stabilization mechanisms, which cause τ Glu to stop changing when there is sufficient overall NMDAR conductance. Note, τ Glu stabilization generally occurs more quickly for smaller τ D , and even at high τ D (18 to 100 ms), τ Glu still stabilizes for certain values at the end of the simulation, despite limited overlap of the periodic signals. (G-I) Same simulations with stability mechanisms as in Figure (D-F), except with initial τ Glu = 50 ms.

Discussion
SITDL explores the potential for a single synapse to learn the timing of input signals through changes of its receptor dynamics. In the computationally simple SITDL model we show that single synapse can "learn" the timing difference between glutamate and dendritic voltage signals. The functioning principle of single-compartment SITDL model is to reach optimal combination of fast and slow NMDAR subunits and corresponding glutamate gate conductance characteristic time, τ Glu , through minimization of NMDAR gate conductance mismatch, (g Glu − g V ). τ Glu stabilizes when optimal Ca 2+ influx is achieved. Decreases in τ Glu correspond to a replacement of slow NMDARs by fast NMDARs (Fig. 4), as is typical during neural development. However, depending on initial conditions, such as initial τ Glu value and dendritic delay τ D , τ Glu can also increase over time (Fig. 3), relating to the observed bidirectional switching of fast and slow NMDARs in hippocampus 13,51 . NMDARs are particularly suited to mediate relations of the SITDL hypothesis because of their largely independent glutamate and voltage gates. SITDL mechanisms would likewise function under the assumption of NMDAR's slow Mg 2+ unblock 52 . Furthermore, it is possible that other receptors with coincidence detector properties, such as the inositol 1,4,5-trisphosphate receptor (IP3R) 53,54 , could be involved in SITDL-like mechanisms. Additionally, a variety of ligand-gated ion channels also have active voltage-sensitive conductances that could dispose them toward various ways of regulation in sensitive ranges of voltage. While the assumptions of the SITDL model have not yet been tested in physiological experiments, the functional principles may significantly expand synaptic capabilities in both biological and artificial systems.
Much is still unknown about NMDAR properties, and their activity-dependent modifications. It is possible that specific changes in NMDAR conformation, occurring independently for glutamate and voltage gate activation, can convey the mismatch between gate conductances through interfacing with cellular signaling pathways. For instance, if τ Glu is too large, with too many slow NMDARs at a synapse, then glutamate gate activation peaks www.nature.com/scientificreports/ after the voltage gate activation, producing a specific combined conformational state in NMDARs, particularly slow ones, that indicates the NMDAR gate conductance mismatch. This specific conformational state could affect the NMDAR's anchoring with post-synaptic density proteins, such as PSD-95. Therefore, NMDARs with the highest mismatch, as indicated by their specific combined conformational state, can be removed from the membrane, and replaced with other NMDARs. This NMDAR cycling process can continue until the balance of slow and fast receptors is found, where glutamate and voltage gate conductances peak at the same time, for an optimal τ Glu . The different nanodomain organization of slow and fast NMDARs in the postsynaptic density could also provide more dynamic cycling and replacement of NMDARs 55,56 , such that unstable slow NMDARs with high mismatch would be more likely replaced with fast NMDARs, and vice versa. Moreover, SITDL provides potential mechanisms for new learning rules in ANNs, which currently employ primarily synaptic weight changes. The SITDL model can increase overlap of g Glu and g V in a synapse that receives very similar periodic signals with small relative time shifts (Fig. 2). This may be particularly useful for developing synapses in delay line systems like those of mammalian and avian auditory brainstem and sensory circuits of weakly electric fish 57,58 , all of which rely on very precise timing of signals. Some delay line systems resolve temporal disparities in the microsecond range even with just a single neuron 21 , and some depend on NMDARs in development 17 . Notably, there have been few computational models for achieving precise timing in delay line systems. One such model explores temporal precision in the barn owl auditory system 59 , using unsupervised Hebbian learning rules and a broad random distribution of transmission delays. While this may show how delay line systems with large numbers of neurons are able to achieve such temporal precision, it does not explain how it can be achieved with much more limited numbers of neurons and transmission delays, such as in pre-pacemaker nucleus neurons of weakly electric fish 21 . Furthermore, the vast diversity of neuronal morphology and circuit organization, as found in cerebral cortex and cerebellum 60,61 , may require fine-tuning of synaptic input timing. Even without τ Glu stabilization mechanisms, the SITDL model can still increase overlap of NMDAR gate conductances with larger dendritic delays, such as for some τ D greater than 45 ms, which provide very little overlap of glutamate and voltage signals (Figs. 2C and 3). Due to periodicity of the signals, less frequent but repeated overlap can cause significant changes in τ Glu . With τ Glu stabilization mechanisms depending on overall changes in τ Glu , and overall NMDAR conductance, the SITDL model stabilizes τ Glu even with quite small overlap between glutamate and voltage signals. Changes in τ Glu are bidirectional (Fig. 3), and represent replacement of fast and slow NMDARs (Fig. 4). This relates to mechanisms of rapid bidirectional switching in NMDAR subunit compositions of developing hippocampus 13 .
Achieving greater gate conductance overlap with very limited overlap of signals may be quite useful, especially in memory formation. For instance, in a neuron that receives highly delayed copies of the same or similarly periodic signal, any memorized timing differences may represent the commonly occurring inter-spike intervals of the signal. We show in SITDL simulations that it is possible to memorize and reconstruct the original synaptic glutamate signal (Figs. 5, 6). Notably, synapses which stabilized at the smallest τ Glu values provided the largest NMDAR conductances, saving them from synaptic elimination.
Potential tests of SITDL in biological systems. SITDL mechanisms would be particularly valuable in systems that rely on precise timing and sequence memorization, such as the neural circuits involved in avian song and human speech learning 40,62 , as well as spatial map formation 63 . If we consider the auto-associative networks of the CA3 region of the hippocampus, where synapses can exhibit bidirectional changes in NMDAR subunit compositions 51 , then the SITDL mechanism could potentially be used for sequence memorization in single neurons. Figure S5 shows one potential mechanism of sequence learning in CA3, with each synapse learning the timing between two distinct inputs. Notably, in hippocampal neurons, the distribution and nanoscale organization of GluN2B subunits vary significantly between proximal and distal synapses 64 . Likewise, spike timing dependent plasticity (STDP) rules of synapses also vary in a location-dependent manner 65 . SITDL mecha- Figure 6. Results of Recall Phase simulations for synapses with and without SITDL mechanisms. There are 4 different sets of synaptic simulations, each with their own starting conditions: 2 sets of SITDL simulations with τ D step = 2 ms and 8 ms, and 2 sets of simulations without SITDL mechanisms with τ D step = 2 ms and 8 ms. For all simulation sets, τ D values are initially distributed uniformly with the corresponding τ D step, and all initial τ Glu = 20 ms. For each simulation set, following Learning Phase simulations, synaptic elimination, and Recall Phase simulations, a recall signal is obtained by summating the resulting Recall Phase voltage signals (Fig. 5C), shifted by corresponding τ D (A-D traces in color), over all stabilized synapses. Note, to provide better visual comparison of original and recall signals, the original signal was slightly shifted relative to the recall signal, with a delay of 5.5 ms for (A,D), and 5 ms for (B,C). (A) The recall signal for SITDL simulation set with τ D step = 2 ms very closely resembles the original glutamate signal that was used during Learning Phase. (B) Recall signal for SITDL simulation set with τ D step = 8 ms is less similar to the original glutamate signal, because of the much larger τ D step used during Learning Phase. (C,D) Recall signals for no-SITDL simulation sets with τ D step = 8 ms and 2 ms, respectively. There are no changes in τ Glu in these simulations, thus synapses are unable to achieve greater overlap of glutamate and voltage gate conductances. Even with synapses that were stabilized, the peaks of the recall signal are significantly shifted from the peaks of the original glutamate signal. (E) Estimates of MI between original glutamate signal and recall signal for each simulation set, using MI estimator AIMIE 37 . SITDL simulation sets, particularly the one with smaller τ D step (2 ms), provide greater MI estimates than no-SITDL simulation sets, suggesting greater degree of similarity between their reconstructed recall signal and the original glutamate signal. www.nature.com/scientificreports/ nisms could potentially be used to establish specific NMDAR expression, and consequently the synaptic timing required for the distinct location-dependent STDP rules, along the dendritic arbor.
To test for the existence of SITDL mechanisms in real neurons, in vitro studies of neuronal cultures and brain slices, such as those of hippocampus and developing neocortex could be useful. In particular, it would be interesting to observe axonal and synaptic activity across multiple points in developing multi-synaptically connected neurons, through multi-electrode recordings or with calcium and voltage indicators. Furthermore, neuronal activity might be manipulated, pharmacologically or with stimulating electrodes, alongside single particle tracking of different NMDAR subunits in post-synaptic membranes 66 , to provide more insight into possible existence of neuronal SITDL mechanisms. Molecular analyses of NMDARs in these cases, including conformational changes and phosphorylation site alterations under different stimulation protocols, might provide a deeper understanding of the potential mechanisms involved. For instance, it is possible that NMDAR gate conductance mismatch could be due to the different conformational changes caused by glutamate binding and depolarization of an NMDAR. These specific conformational changes could alter the NMDAR's phosphorylation sites, modifying its dynamics and anchoring stability in the post-synaptic membrane, and leading to potential subunit switching.