Astrocyte-mediated spike-timing-dependent long-term depression modulates synaptic properties in the developing cortex

Astrocytes have been shown to modulate synaptic transmission and plasticity in specific cortical synapses, but our understanding of the underlying molecular and cellular mechanisms remains limited. Here we present a new biophysicochemical model of a somatosensory cortical layer 4 to layer 2/3 synapse to study the role of astrocytes in spike-timing-dependent long-term depression (t-LTD) in vivo. By applying the synapse model and electrophysiological data recorded from rodent somatosensory cortex, we show that a signal from a postsynaptic neuron, orchestrated by endocannabinoids, astrocytic calcium signaling, and presynaptic N-methyl-D-aspartate receptors coupled with calcineurin signaling, induces t-LTD which is sensitive to the temporal difference between post- and presynaptic firing. We predict for the first time the dynamics of astrocyte-mediated molecular mechanisms underlying t-LTD and link complex biochemical networks at presynaptic, postsynaptic, and astrocytic sites to the time window of t-LTD induction. During t-LTD a single astrocyte acts as a delay factor for fast neuronal activity and integrates fast neuronal sensory processing with slow non-neuronal processing to modulate synaptic properties in the brain. Our results suggest that astrocytes play a critical role in synaptic computation during postnatal development and are of paramount importance in guiding the development of brain circuit functions, learning and memory.

Introduction Synaptic long-term plasticity, defined as the activity-dependent change in the strength or efficacy of the synaptic connection between a pre-and postsynaptic neuron, is expressed in the brain in diverse forms across multiple timescales [1]. Action potential (or spike) timing is one of the many factors governing synaptic plasticity induction [2,3]. In spike-timing-dependent plasticity (STDP), the order and precise temporal difference between pre-and postsynaptic action potentials determine the direction and magnitude of long-term plasticity. Depending on the form of synaptic plasticity and the brain area, a large number of cellular and molecular level mechanisms are involved [4][5][6][7]. In the developing mouse barrel area of the somatosensory cortex, spike-timing-dependent long-term depression (t-LTD) [8] at layer 4 (L4) to layer 2/3 (L2/3) synapses has been shown to require activation of presynaptic mechanisms [9][10][11][12][13][14] and involve astrocytic functions [15]. This t-LTD has been shown to emerge in the first postnatal week, be present during the second week, and disappear in the adult, whereas spike-timingdependent long-term potentiation (t-LTP) persisted into adulthood [10]. Long-term depression may provide an important mechanism for synapse pruning and subsequent neuron and circuit remodeling during postnatal development [16].
Astrocytes, a type of non-neuronal cells in the mammalian brain, are recognized as important homeostatic, metabolic, and neuromodulatory elements that are also coupled to the neurovascular system [17,18]. In the developing central nervous system, astrocytes promote the formation of excitatory synapses and the establishment of synaptic connectivity [19]. Astrocytes can also sense and modulate synaptic functions [20]. Astrocytes maintain glutamatergic synaptic transmission by glutamate uptake [21] and clear excess extracellular potassium ions (K + ) to spatially transfer K + from regions of elevated concentration to regions of lower concentration [22]. In addition, there is ample evidence to indicate that astrocytes actively contribute to the information processing capabilities of neural circuits and ultimately affect animal behavior [23,24]. Astrocytes have, for example, been shown to influence brain state transitions [25], promote the coordinated activation of neuronal networks [26], and modulate sensory-evoked neuronal network activity [27] and brain rhythms during sleep [28]. Recent research has the potential to revolutionize our current understanding of the role of astrocytes in the modulation of brain network activity [17,29].
Nevertheless, the central questions still remain: Do cortical astrocytes in vivo have subcellular mechanisms capable of synapse modification at fast enough timescales comparable to neuronal ones? Does this modulation depend on the brain area and circuitry in question? Is this modulation significant only in developing brain circuits or does it also happen in mature circuits? The answers to these questions will significantly increase our understanding of mammalian neocortical network functioning.
Computational modeling is an important complementary method for linking the dynamics of different biochemical and biophysical reactions and processes together and for unraveling the complexity of synaptic functions. Our goal here is to better understand through computational modeling the role of cortical astrocytes in sensory processing, particularly in synaptic plasticity, during postnatal development. To address this question we propose a new biophysicochemical model of a somatosensory cortical L4 to L2/3 synapse and study the role of astrocytes in t-LTD in vivo. We made several assumptions based on the experimental electrophysiological, Ca 2+ imaging, and other data (see Materials and methods and S1 Appendix). The computational model was built in component-by-component manner for the presynaptic L4 spiny stellate cell and postsynaptic L2/3 pyramidal cell as well as for the nearby fine astrocyte process. After careful validation of each model component, all the components were brought together to describe all the necessary elements of a somatosensory cortical synapse. The integrated model takes into account the well-established biophysical and biochemical mechanisms for this particular synapse, such as the voltage-gated ion channels, transmitter-gated receptors, Ca 2+ -mediated signaling pathways including the neuronal endocannabinoid and astrocytic inositol 1,4,5-trisphosphate (IP 3 ) receptor (IP 3 R) signaling, as well as other crucial subcellular mechanisms. These mechanisms are described using deterministic differential equations. The integrated model is carefully validated against experimental data on synaptic plasticity [11,12,15]. Here we show that cortical astrocytic Ca 2+ dynamics can be modified by presynaptic L4 spiny stellate cell and postsynaptic L2/3 pyramidal cell activity through the endocannabinoid signaling pathway. The subsequent downstream signaling pathways in astrocytes have an influence on synaptic longterm plasticity, particularly on the t-LTD in somatosensory cortex, through presynaptic Nmethyl-D-aspartate receptors (NMDARs) and calcineurin (CaN) signaling. Our study provides several predictions that can be tested in future electrophysiological, Ca 2+ imaging, and molecular biology experiments.

Results
We simulated a synapse model containing neuronal pre-and postsynaptic terminals and a fine astrocyte process. Specifically, our computational model includes the axonal compartment of a presynaptic L4 spiny stellate cell, the dendritic and somatic compartments of a postsynaptic L2/3 pyramidal cell, and the nearby fine astrocyte process. Several previous modeling studies [52][53][54][55] have had an influence on our synapse modeling project and the choices we made during the work. In our in silico experiments, we studied which mechanisms are important in the induction of t-LTD at L4-L2/3 synapses in somatosensory cortex, including key Ca 2+ -dependent intracellular processes. We used stimulation protocols equivalent to the protocols applied in electrophysiological experiments in vitro and in vivo to activate our in silico synapse model [12]. We showed that t-LTD at an L4-L2/3 synapse can be explained by the activation of Ca 2+ -dependent mechanisms in the fine astrocyte process and this further has an influence on the probability of neurotransmitter release in the presynaptic neuron through NMDARs and calcineurin signaling. In the absence of the Ca 2+ -dependent mechanism in the fine astrocyte process, the synapse did not show t-LTD similarly to experimental data.

Synapse model dynamics before, during, and after t-LTD induction: Fitting the model to experimental data
In our simulations, we closely followed experimental stimulation protocols [12]. During our stimulation protocol before t-LTD induction, we simulated our synapse model with five pulses of presynaptic stimulus at a frequency of 0.2 Hz (Fig 2A). Our t-LTD induction protocol consisted of 100 post-pre pairings at a frequency of 0.2 Hz where a postsynaptic stimulus occurred between 10 ms and 200 ms before a presynaptic stimulus, thus the temporal difference (ΔT) had values between −10 ms and −200 ms ( Fig 2H). The protocol after t-LTD induction included five pulses of presynaptic stimulus at a frequency of 0.2 Hz (Fig 2O), similarly as with the protocol before t-LTD induction. All these different stimuli triggered changes in the preand postsynaptic membrane potentials, similarly to experimental data [12,15] (Fig 2B, 2C, 2I, 2J, 2P and 2Q), that led, for example, to the opening of pre-and postsynaptic Ca 2+ channels and glutamate release from the presynaptic neuron (Fig 2D-2G, 2K-2N and 2R-2U). The simulated presynaptic Ca 2+ concentration values followed the experimental values [81][82][83][84] (Figs 2D, 2K and 2R and 3F and 3M). The glutamate concentration in the synaptic cleft increased to about 500 μM after stimuli, which is close to the measured experimental values [85] (Fig 2G,  2N and 2U). The release probability of presynaptic glutamate vesicles and the concentration of glutamate in the synaptic cleft were the lowest for the shortest ΔT due to ongoing astrocytemediated molecular dynamics during depression (Fig 2L, 2N, 2S and 2U). The effect of depression is clearly seen after t-LTD induction (Fig 2Q and 2S-2U).
During the t-LTD induction protocol, the released glutamate in the synaptic cleft activated AMPARs, NMDARs, and mGluRs in the dendritic membrane of the postsynaptic neuron, in addition to presynaptic NMDARs. The activation of these postsynaptic receptors together with Ca 2+ influx via Ca LHVA and Ca LLVA channels into the postsynaptic neuron induced a G-protein signaling cascade that activated phospholipase C (PLC) (Fig 3A and 3H). This led to the production of biologically realistic concentrations of diacylglycerol (DAG) and IP 3 in the postsynaptic neuron. Then IP 3 activated, together with Ca 2+ , Ca 2+ -induced Ca 2+ release via IP 3 Rs Presynaptic membrane potential depends on currents via Ca NHVA , Na + , and K + channels as well as via NMDARs. Presynaptic action potential and Ca NHVA -and NMDAR-mediated Ca 2+ concentrations together with the influence of CaN affect the vesicular release. (2) The released glutamate in the synaptic cleft activates postsynaptic mGluRs, NMDARs, and AMPARs in addition to presynaptic NMDARs. (3) Postsynaptic membrane potential in the soma depends on currents via Na + , Na P , and K DR channels, whereas postsynaptic membrane potential in the dendrite depends on currents via Ca LHVA , Ca LLVA , Na + , and K A channels as well as via NMDARs and AMPARs. The activation of postsynaptic mGluRs and NMDARs, together with the Ca LHVA -and Ca LLVA -mediated Ca 2+ influx, triggers a G-protein signaling cascade where GαGTP dissociates from mGluR-bound Gβγ and activates PLC and production of DAG and IP 3 . Increases in Ca 2+ and IP 3 concentrations activate Ca 2+ release via IP 3 Rs from the ER to the cytosol. On the other hand, PMCA and SERCA pumps transfer Ca 2+ away from the cytosol and leak fluxes transfer Ca 2+ back to the cytosol. The production of DAG leads to a production of endocannabinoid 2-AG. (4) Endocannabinoid 2-AG released from the postsynaptic neuron binds to the astrocytic CB 1 Rs and triggers Ca 2+ signaling in the astrocyte. We modeled this step by directly modifying astrocytic IP 3 concentration based on the postsynaptic 2-AG concentration. (5) Astrocytic IP 3 and Ca 2+ activate similar ER-related events as in the postsynaptic neuron. Astrocytic Ca 2+ increase then induces glutamate exocytosis to the extrasynaptic space. (6) Glutamate in the extrasynaptic space and the spillover of glutamate from the synaptic cleft activate presynaptic NMDARs. (7) Presynaptic NMDAR-mediated Ca 2+ concentration activates CaN, and CaN has an effect on vesicular release together with presynaptic action potential and Ca NHVA -mediated Ca 2+ concentration.   [12]. The simulation results are shown for six key model variables during the first two stimulus pulses of our protocol before t-LTD induction in (B-G), during a single post-pre from the ER to the cytosol in the postsynaptic neuron, which led to an increase in Ca 2+ concentration in the cytosol. The production of DAG, on the other hand, resulted in a production of 2-arachidonoylglycerol (2-AG) (Fig 3B and 3I), and ultimately in the release of endocannabinoid 2-AG from the postsynaptic neuron. All these are well-established signaling pathways known to exist in cortical neurons.
Endocannabinoid 2-AG can bind to the astrocytic CB 1 Rs and trigger Ca 2+ signaling in the astrocyte [15,42]. We modeled this step by directly modifying IP 3 concentration in the astrocyte based on the postsynaptic concentration of 2-AG ( Fig 3B, 3C, 3I and 3J), followed by an increase in astrocytic Ca 2+ concentration (Fig 3D and 3K) and ultimately glutamate exocytosis, thus inducing glutamate release from the astrocyte to the extrasynaptic space ( Fig 3E and 3L). We chose the astrocytic Ca 2+ threshold for glutamate release based on experimental data [86]. Similarly, we chose the maximum value of glutamate concentration in the extrasynaptic space based on experimental findings [87].
Astrocytes have been shown to have an effect on presynaptic glutamate release by modifying release probabilities [15,36,88]. In somatosensory cortex, astrocytes have exhibited reduction in the presynaptic release probabilities as a response to the t-LTD induction protocol [15]. In our synapse model, glutamate release from the presynaptic neuron depended, among other things, on the presynaptic Ca NHVA -and NMDAR-mediated Ca 2+ concentrations (Figs 2D, 2K and 2R and 3F and 3M), release probability of presynaptic glutamate vesicles (Fig 2E, 2L and 2S), presynaptic calcineurin concentrations [14] (Fig 3G and 3N), and fraction of presynaptic glutamate release inhibition (f pre , see Materials and methods and S1 Appendix). The presynaptic NMDARs were activated by the glutamate in the extrasynaptic space and the spillover of glutamate from the synaptic cleft. Our simulations showed that the glutamate in the extrasynaptic space substantially increased the presynaptic NMDAR-mediated Ca 2+ concentration (Fig 3F and 3M).

t-LTD amplitude depends on the temporal difference between post-and presynaptic activity: Confirming the broad t-LTD time window
In our in silico experiments, we followed the experimental t-LTD stimulation protocols [12]. First, we estimated the amplitude of the excitatory postsynaptic potential (EPSP) before t-LTD induction when the stimulation protocol consisted of only a presynaptic stimulus repeated five times at a frequency of 0.2 Hz (Fig 2A). The EPSP before t-LTD induction is presented in Figs  2C and 4A (right). Then t-LTD was induced by the post-pre pairing protocol consisting of a postsynaptic stimulus followed by a presynaptic stimulus with a temporal difference ΔT from −10 ms to −200 ms and the pairing was repeated 100 times at a frequency of 0.2 Hz (Fig 2H). In this case, the postsynaptic action potential in the soma was followed by an EPSP which is shown during one presynaptic stimulus in Fig 4A (left) for ΔT from −10 ms to −200 ms (see also Fig 2J). Presynaptic activity and thus EPSPs were delayed by ΔT in respect to the postsynaptic action potential. After t-LTD induction, we estimated the changes in the EPSPs by stimulating the synapse model with the same protocol as before t-LTD induction, including only a pairing with five different ΔT occurring about in the middle of the 100 post-pre pairings of the t-LTD induction protocol in (I-N), and during a single stimulus pulse of our protocol after t-LTD induction in (P-U, note the different x-axis in U). The presynaptic membrane potential (V pre ) in (B, I, P), postsynaptic membrane potential in the soma (V soma,post ) in (C, J, Q, note the different y-axis in Q), presynaptic Ca NHVA -mediated Ca 2+ concentration ([Ca 2+ ] CaNHVA,pre ) in (D, K, R), release probability of presynaptic glutamate vesicles (P rel,pre ) in (E, L, S), fraction of releasable presynaptic vesicles (R rel,pre ) in (F, M, T), and glutamate concentration in synaptic cleft ([Glu] syncleft ) in (G, N, U) responded to the stimuli shown in (A, H, O), respectively. In (J), the postsynaptic action potential in the soma was followed by an EPSP with a delay corresponding to ΔT. The lowest release probability of presynaptic glutamate vesicles in (L, S), the lowest glutamate concentration in synaptic cleft in (N, U), and the highest fraction of releasable presynaptic vesicles in (M, T) were obtained with the shortest ΔT due to the astrocyte-mediated signaling during depression.
https://doi.org/10.1371/journal.pcbi.1008360.g002 The post-pre pairings induced presynaptically expressed t-LTD, sensitive to the temporal difference between the pre-and postsynaptic activity (ΔT). The strongest LTD change was observed for the shortest ΔT: the EPSP decreased from 4.9 mV (before t-LTD induction) to 3.1 mV (after t-LTD induction) (Fig 4A (right)). The time courses of the presynaptic glutamate release inhibition fraction (f pre ) for ΔT from −10 ms to −200 ms are shown in Fig 4B. A shorter astrocytic IP 3 concentration ([IP 3 ] astro ) in (C, J), and astrocytic Ca 2+ concentration ([Ca 2+ ] astro ) in (D, K). The postsynaptic 2-AG concentration shown in (B, I) triggered an increase in the astrocytic IP 3 concentration shown in (C, J) which was followed by an increase in the astrocytic Ca 2+ concentration shown in (D, K). (E, L) After the astrocytic Ca 2+ concentration reached the threshold, astrocyte released a fixed amount of glutamate to the extrasynaptic space and this can be seen as an increase in the glutamate concentration in the extrasynaptic space ([Glu] extsyn ). (F, M) Both glutamate in the extrasynaptic space and the spillover of glutamate from the synaptic cleft were able to activate presynaptic NMDARs. The smaller changes in the presynaptic NMDAR-mediated Ca 2+ concentration ([Ca 2+ ] NMDAR,pre ) occurred due to the spillover of glutamate from the synaptic cleft, and the larger changes occurred due to the glutamate in the extrasynaptic space. (G, N) The presynaptic NMDAR-mediated Ca 2+ influx increased the presynaptic calcineurin concentration ([CaN] pre ).

Fig 4. Shorter temporal difference between pre-and postsynaptic activity leads to stronger t-LTD through astrocyte-mediated cellular and subcellular mechanisms.
The t-LTD stimulation protocols were obtained from experimental literature [12] and their use in in silico modeling is shown in Fig 2A, 2H and 2O. (A) In the left, the postsynaptic membrane potential in the soma is shown during a single post-pre pairing of the t-LTD induction protocol with a temporal difference ΔT between −10 ms and −200 ms at every 10 ms. The postsynaptic stimulus evoked a somatic action potential followed by an EPSP generated by the presynaptic stimulus. The longer the temporal difference ΔT, the longer the delay for the EPSP. In the right, the postsynaptic membrane potential in the soma, in other words in this case the postsynaptic EPSP generated by the presynaptic stimulus, is shown during a single presynaptic stimulus occurring before (black) and after (color bar) t-LTD induction. The shorter the temporal difference ΔT, the smaller the amplitude of the EPSP. (B) The fraction of presynaptic glutamate release inhibition (f pre ) had the highest values with the shortest ΔT, i.e. the strongest t-LTD, during the 100 post-pre pairings in the t-LTD induction protocol. Color bar is given in (A). (C) The final value of −f pre is shown as a function of ΔT. (D) The ΔEPSP percentage is shown as a function of ΔT. We calculated the ΔEPSP percentage for every ΔT as the percentage change between the somatic EPSP amplitude evoked by the presynaptic stimulus occurring before t-LTD induction (shown in (A, right) as black) and the somatic EPSP amplitude evoked by the presynaptic stimulus occurring after t-LTD induction (shown in (A, right) with different colors given in the color bar). Our synapse model confirmed the experimental data [12]. The shorter the temporal difference ΔT, the stronger the t-LTD. ΔT led to a larger increase in f pre (Fig 4B), and thus stronger t-LTD ( Fig 4D). The dependence of final f pre values on ΔT is shown in Fig 4C. The fraction f pre had different resulting values depending on ΔT used (Fig 4C), for example f pre = 0.5 for ΔT = −10 ms, f pre = 0.34 for ΔT = −100 ms, and f pre = 0.03 for ΔT = −200 ms. The changes in the EPSP, estimated in Fig (Fig 4D). Thus the stimulation protocol induced t-LTD for ΔT values shorter than −200 ms and t-LTD was the strongest for the shortest ΔT, which was consistent with the experimental results [12] ( Fig  4D). The broad time window for t-LTD in somatosensory cortex has been reported in several experimental studies [8,12].

t-LTD requires astrocytic signaling and presynaptic NMDARs
Previous experimental studies have reported that presynaptic GluN2C/D-containing NMDARs are required for t-LTD, whereas postsynaptic GluN2B-containing NMDARs are necessary for t-LTP at the vertical L4 input onto L2/3 neuron [9,12,15,[89][90][91]. Our model simulations showed that blocking postsynaptic GluN2B-containing NMDARs, thus changing their conductance to zero, did not prevent the increase in fraction of presynaptic glutamate release inhibition (f pre ) (Fig 5B and 5C), and therefore did not abolish t-LTD (  The EPSP percentage is given for all six models described in (A). We calculated the EPSP percentage for every ΔT by normalizing the EPSP amplitude occurring after t-LTD induction by the EPSP amplitude occurring before t-LTD induction, and multiplied them by 100%.
Previous experimental studies have also shown that t-LTD requires astrocytic CB 1 R activation by neuronal endocannabinoid release followed by an increase in astrocytic Ca 2+ signaling and the exocytosis of glutamate from astrocytes [15]. The released glutamate then activates presynaptic NMDARs and leads to t-LTD [15]. We therefore tested whether interfering with the astrocytic activity leads to the inhibition of t-LTD by blocking the astrocyte, thus keeping the astrocyte model in a steady state by setting all the astrocytic differential equations to zero in our synapse model. The simulation results showed that f pre stayed at low levels (Fig 5B and  5C) and did not lead to t-LTD (Fig 5A (bottom left) and Fig 5D). Thus, our synapse model confirmed the experimental findings that blocking the fine astrocyte process activity entirely prevented t-LTD.
In addition, we tested the model by applying either a presynaptic (Fig 5A (bottom middle)) or a postsynaptic (Fig 5A (bottom right)) stimulus at a frequency of 0.2 Hz for 500 s, thus repeating both stimulation protocols 100 times. In both cases, f pre did not increase substantially, failing to induce t-LTD ( Fig 5D). Our synapse model confirmed the experimental data that an unpaired synaptic pathway remains unmodified [12].

Astrocytes sense the temporal difference of t-LTD and modify their Ca 2+ signaling
Finally, we studied in more detail how Ca 2+ concentration behaves in the fine astrocyte process during the t-LTD induction protocol depicted in Fig 2H. The delay in the astrocytic Ca 2+ peak responses to the post-pre pairing onset varied with the temporal difference of post-pre pairings ( Fig 6A). The delay increased with the lengthening of the post-pre pairing temporal difference (Fig 6B-6D).
We then addressed the occurrence of peaks in the fine astrocyte process. Our synapse model showed that astrocytic Ca 2+ peaks for different ΔT of the t-LTD induction protocol occurred within 2 s of each other in the beginning of the protocol (Fig 6B), whereas the peaks occurred within 10 s in the middle of the protocol ( Fig 6C) and within 14 s towards the end of the protocol (Fig 6D). The normalized mean peak value of the astrocytic Ca 2+ concentration increased with the shortening of the temporal difference between post-pre pairings, having the highest values during the first 50 post-pre pairings and the lowest during the last 50 post-pre pairings in the t-LTD induction protocol ( Fig 6F). It is of interest to explore these in silico results further in future wet-lab experiments to make it possible to build more sophisticated and biologically relevant models for astrocytes.
Previous experimental results have shown that Ca 2+ transients do not occur at a certain fixed time point after each individual post-pre pairing, but are rather evenly distributed in the 5 s long period between each post-pre pairing [15]. We confirmed this experimental finding by calculating the probability of astrocytic Ca 2+ peaks occurring in the 5 s long period between each pairing with different ΔT (Fig 6G). One of the reasons behind this distribution is that astrocytes are activated slower than the individual post-pre pairings because of the slow endocannabinoid signaling [15]. Our synapse model predicted that the astrocytic Ca 2+ concentration oscillated every 13 s during the t-LTD induction protocol (13.34 s for ΔT = −10 ms and 13.67 s for ΔT = −200 ms), which is close to the reported experimental oscillation rate of every 15 s for ΔT = −25 ms [15]. Note that both the experimental and computational values are about 3 times longer than the length between each post-pre pairings. The number of times the astrocyte released glutamate during the whole 100 post-pre pairings in the t-LTD induction protocol increased with the shortening of ΔT (Fig 6E and 6H). In our synapse model, this is due to the fact that astrocytic Ca 2+ peaks were slightly higher with shorter ΔT. There is experimental evidence showing that Ca 2+ peaks are not higher with shorter ΔT but instead Ca 2+ transients are more frequent and an individual Ca 2+ transient lasts longer [15]. Our astrocyte model is based on the same mechanisms as the models published so far [78,79]. This issue with more frequent and longer Ca 2+ transients clearly requires further experimental clarification, so that future computational models may be extended to incorporate more realistic Ca 2+ transients using available simulation tools, for example [92].
In summary, the model simulations confirmed several experimentally obtained results, such as t-LTD sensitivity to ΔT [12,15] and the role of astrocytic signaling in t-LTD [15]. Moreover, the model simulations predicted the time courses of astrocytic Ca 2+ signals and the putative roles and time courses of presynaptic mechanisms in t-LTD. These predictions will be useful in planning the future studies of astrocytes and synapses in somatosensory cortex in vivo.

Discussion
Astrocytes have been shown to dynamically modulate synaptic transmission and plasticity in some cortical synapses, but how this occurs in time and space has remained unclear [15,17,27,93]. We demonstrate with a new somatosensory cortical synapse model that a well-established feedback signal from a postsynaptic neuron to a presynaptic neuron via a fine astrocyte process can induce, maintain, and modulate spike-timing-dependency of long-term depression during postnatal development at cortical layer 4 to layer 2/3 synapses. This modulation occurs through astrocyte-mediated molecular mechanisms to the presynaptic axonal terminal. We predict for the first time the dynamics of these molecular mechanisms underlying spiketiming-dependent LTD and link complex biochemical networks at the pre-and postsynaptic as well as astrocytic sites to the electrophysiology and time window of spike-timing-dependent plasticity induction at vertical L4-L2/3 synapses [12,15]. The removal of any of the key mechanisms, including the astrocytic mechanisms, impaired synaptic t-LTD. Our results indicate that multiple biophysical and biochemical plasticity mechanisms at the L4-L2/3 neuronal synapse and nearby fine astrocyte process contribute to enabling synaptic LTD in a developing somatosensory cortex.
Our study highlights several important advancements in neuroscience. First, we link together the dynamics of known cellular and molecular players of t-LTD during postnatal development and describe each model component by mathematical equations and data from a multitude of experimental and modeling studies. Second, we combine unique experimental results on the time-dependency of t-LTD in a developing somatosensory cortex, obtained by two independent research groups [11,12,14,15], to validate our model. Third, our analysis using the biophysically and biochemically detailed synapse model confirms the experimental findings on astrocytes' ability in setting the temporal difference of t-LTD at L4-L2/3 synapses [15]. In summary, we confirm with our in silico synapse model the following experimental findings and predictions (1-4).
1. The fine sensitivity of t-LTD to the temporal difference in a developing somatosensory cortex is achieved through complex molecular signaling, similarly to experimental data and predictions [11,12,15]. PLOS COMPUTATIONAL BIOLOGY 2. At the L4 spiny stellate cell-L2/3 pyramidal cell synapse, t-LTD is orchestrated through the postsynaptic release of endocannabinoid molecules (agonist 2-AG) detected by CB 1 Rs on the fine astrocyte process [15].
3. Astrocytic Ca 2+ transients induced by endocannabinoids and subsequent exocytosis of glutamate from the fine astrocyte process are appropriate to induce and maintain t-LTD (comparable to experimental validation data [15]).
4. Glutamate release from the fine astrocyte process can be detected by presynaptic NMDARs at the time courses appropriate for the modulation of synaptic release through calcineurinrelated signaling [14].
We modeled all the above-mentioned mechanisms using biologically realistic time constants validated against published experimental data (see Materials and methods and S1 Appendix). The predictions made by our synapse model are readily available for further experimental wet-lab testing. In addition, we provide all mathematical equations and their relationships, all parameter values, all references used in the model construction, and commented code upon publication in order to enable the reproduction of our results and facilitate reproducible science [76][77][78][79]94].
In a developing somatosensory cortex, t-LTD has been shown to require activation of CB 1 Rs by postsynaptically released endocannabinoids, and increased astrocytic Ca 2+ concentration [15]. However, the spatial location and distribution of the CB 1 Rs is under debate. Earlier it was assumed that t-LTD requires CB 1 Rs located on the presynaptic neuron [95], but more recent evidence from several brain areas and spinal cord shows that CB 1 Rs are also located on astrocytes [15,[96][97][98]. Agonists of CB 1 Rs have been found to evoke Ca 2+ transients in astrocytes [97] and in the micro-domains of astrocyte processes [15,98]. Furthermore, it has been shown that a prerequisite for t-LTD in the somatosensory cortex [15], and also in the hippocampus [96], are astrocytic CB 1 Rs, not the presynaptic CB 1 Rs. Based on these most recent findings we modeled the postsynaptically released endocannabinoid activation only on the astrocytic CB 1 Rs [15]. In addition, we made an assumption that an increase in astrocytic Ca 2+ levels, due to endocannabinoids in our model, is mediated by IP 3 Rs on the ER membrane [40] and subsequent Ca 2+ -dependent glutamate exocytosis [42,99,100]. Recently, studies have found multiple types of Ca 2+ signals in astrocyte processes [44], also in somatosensory cortex in vivo [15,101,102]. These multiple types of Ca 2+ signals may be explained by the activation of different subtypes 1, 2, and 3 of IP 3 Rs [40]. Although evidence against [103,104] and for [47] IP 3 R-mediated Ca 2+ -dependent glutamate exocytosis in plasticity exist, we decided to test with our model whether the kinetics of Ca 2+ -dependent glutamate exocytosis in astrocyte processes can take part in mediating t-LTD observed in somatosensory cortex during postnatal development. The controversial results between different studies may be explained by various factors, including differences in the postnatal developmental stage of the rodent, the brain area, the type of a synapse and brain circuitry, the motility of astrocyte processes, and the experimental conditions as well as the different measurement techniques, including the use of transgenic animals. Furthermore, other Ca 2+ -related mechanisms may coexist in astrocyte processes [40,[42][43][44][45][46][47][48][49] which can also contribute to the modulation of plasticity.
In addition to astrocytic CB 1 Rs, the activation of presynaptic NMDARs is required for t-LTD in the developing somatosensory cortex [9][10][11][12][13][14]. These presynaptic NMDARs have been shown to be tightly linked with presynaptic Ca 2+ , proteins and associated signaling cascades to control the release of neurotransmitters from the vesicles, the size of the vesicle pool, and/or the replenishment of synaptic vesicle pools [14]. The exact signaling between astrocytic CB 1 Rs and presynaptic NMDARs cooperatively leading to synaptic depression is, however, not fully understood. Moreover, the presynaptic NMDAR-dependent LTD (in the vertical pathway) seems to be developmentally regulated and disappears by 3-4 weeks of age in the mouse barrel cortex [10] and visual cortex [105] as well as in the mouse hippocampus [41]. Regardless of the few missing components and debates on the exocytosis of astrocytic glutamate and the presynaptic NMDARs [14,17,29,[106][107][108], we conclude that there is a growing body of evidence suggesting the involvement of astrocytes in t-LTD during postnatal development. Based on recent reconstruction studies of astrocytic morphologies [32] and imaging of IP 3 R-mediated events in fine astrocyte processes [40], astrocytes indeed seem to make an important contribution to synapses. Using computational modeling we present the links between different molecular pathways contributing to the temporal difference of t-LTD and the required time courses of the molecular players. The full synapse model couples the following key signaling cascades: (1) the signaling cascade from the postsynaptic terminal, thus the release of endocannabinoids, to the astrocyte, (2) the signaling cascade from the fine astrocyte process, including Ca 2+ -dependent glutamate exocytosis, to the presynaptic terminal, and (3) the signaling cascade from the presynaptic NMDARs and calcineurin (a protein phosphatase) to the vesicular release of synaptic glutamate. Based on our simulation results endocannabinoid-induced, Ca 2+ -mediated glutamate release from fine astrocyte processes in vivo can thus have a pivotal impact on synaptic properties and thereby on neuronal activity, most profoundly in the developing somatosensory system.
There is plenty of experimental evidence that NMDAR-dependent synaptic plasticity can be induced by several different mechanisms [14]. Studies with neocortical and hippocampal synapses show that presynaptic NMDARs typically induce LTD and postsynaptic NMDARs LTP. This indicates that presynaptic NMDARs control synaptic release and plasticity, particularly in glutamatergic synapses. The expression of presynaptic NMDARs is, however, highly heterogeneous and synapse specific [109]. For example, it has been shown that presynaptic NMDARs can selectively modulate L4-L2/3 synapses in the somatosensory cortex, but not L4-L4 or L2/3-L2/3 synapses [109]. Moreover, presynaptic NMDARs have been shown to operate in unconventional ways in some synapses [110]. At the L4-L2/3 synapse, NMDARs may therefore support a special form of plasticity, also confirmed by our modeling. Taken together, these previous results on the heterogenous expression of presynaptic NMDARs may explain the lack of presynaptic NMDAR-mediated plasticity in some studies. Furthermore, our results suggest that the astrocytic modulation of NMDAR-dependent t-LTD is highly synapse specific, and synapses that do not contain any presynaptic NMDARs cannot implement astrocytic modulation of t-LTD during postnatal development. We are not aware of any study showing that astrocytes are not modulating t-LTD at L4-L2/3 synapses. All these findings highlight the acute need for detailed mechanistic modeling such as our present study where we show astrocytic CB 1 R-and presynaptic NMDAR-dependent t-LTD in a developing somatosensory cortex. It is also likely that some additional plasticity mechanisms could be added to the model or that their role could be fulfilled by multiple redundant parallel plasticity pathways.
For more than ten years, STDP has been suggested to underlie the development of sensory representations and synapse maturation in the somatosensory cortex [111]. In particular, t-LTD at L4-L2/3 synapses in rodents has been shown to be vital for plasticity during postnatal development [89]. A growing body of evidence also suggests that astrocytes have a fundamental role in cortical postnatal development and map plasticity [15,51,112]. It has been suggested that the functional role of astrocytes in t-LTD at developing somatosensory L4-L2/3 synapses might be to act as a time buffer (or, delay factor) for neuronal activity and sensory processing that occurs on a fast millisecond timescale [15]. During these events, fast and correlated neuronal activity is integrated into slower astrocytic Ca 2+ dynamics. It can therefore be speculated that astrocytes monitor, integrate, and modulate the activity of synapses, on longer timescales, to enhance the capacity of information processing in the brain to build a complex cognitive, conscious experience of the acquired sensory information in higher animals and humans. We have here demonstrated how this monitoring, integration, and modulation of activity is orchestrated through biophysicochemical processes in a synapse to induce t-LTD. We conclude that modeling the dynamics of neuron-astrocyte signaling in a synapse can offer profound mechanistic insights into the development of synaptic computation and information processing in sensory systems.
Developing sensory circuits undergo synapse elimination, a process of pruning synapses during development. Synapse elimination is essential for the formation of mature neuronal circuits and proper brain functions in the cerebral cortex. Although less is known about the cortical pruning compared to other areas, a disruption of this process is likely involved in neurodevelopmental diseases such as schizophrenia, autism spectrum disorder, and epilepsy [16]. The specific molecular mechanisms that drive synapse elimination remain mostly unknown. Interestingly, hippocampal astrocytes have been found to contribute to synapse elimination in a subtype 2 IP 3 R-dependent manner through the activation of purinergic signaling [113]. On the other hand, dendritic spines that have contacts with astrocytes have been found to survive longer and be morphologically more mature than those without such contacts [114]. We argue that astrocytes, potentially together with microglia, might contribute to the elimination of synapses at L4-L2/3 using t-LTD. Overall, the astrocytic modulation of STDP may be one important phase in the development of synapses and functional circuits for mature cortical sensory processing.
Different forms of plasticity, including the Hebbian type of plasticity, have been studied both in experimental and computational settings for a long time. There is accumulating evidence that Hebbian framework and plasticity rules may depend on the 3rd and 4th factors, such as neuromodulatory agents or neuroglial cells [115,116]. The 3rd factor is usually included in phenomenological models of synaptic plasticity [115][116][117]. To the best of our knowledge, we present here the first computational study that provides strong supportive evidence on the role of astrocytes and their processes as a putative 3rd factor in t-LTD in the somatosensory cortex during postnatal development. Overall, our results highlight the importance of neuroglial mechanisms in STDP that may complement and stabilize developing somatosensory L4 to L2/3 synapses. The synapses in other cortical layers and brain areas as well as the inhibitory synapses deserve further study, both experimentally and computationally.
We argue that to understand how the brain functions, we need to understand both the structure and function of all the different spatial scales, from genes to the whole brain. Although a great deal of experimental work has been undertaken to study all these different scales, we still have not solved many of the puzzles the brain holds [118]. Computational modeling tightly integrated with experimental data is one of the tools that is used more and more to study brain functions on different scales [52][53][54][55]119]. Modeling approaches bridging different organizational levels and dynamical scales have been increasingly introduced to describe complex neuronal systems [54]. We have here shown how computational modeling can provide important additional insights into the newly developed experimental tools and protocols to study astrocytes and their genetic, molecular, morphological, and physiological profiling in in vivo [120,121]. With computational modeling, we can test different hypotheses, ease the planning of experimental studies, and, especially, explore the role of new mechanisms and their dynamics (temporal behavior) in different experimental settings and brain phenomena.
Our biophysically and biochemically detailed model provides several predictions that could be tested in future wet-lab experiments. The experiments should address the influence of molecular mechanisms, electrophysiological properties, and patterns of neuronal activity on the t-LTD time window (Fig 4). An additional testable key prediction of our work is the astrocytic Ca 2+ signals, shown in Fig 6, by using the same t-LTD induction protocol with different temporal differences. The testing of these predictions requires a combination of electrophysiological, Ca 2+ imaging, and molecular biology techniques. New experimental data could also be helpful in refining some of the model components, particularly the subcellular ones. There are new emerging techniques for single cells developed in the intersection of engineering and biology [122]. These techniques could be used to refine the description of signal transduction pathways, especially the calcineurin-related pathway in the presynaptic terminal. The concentration levels of key molecular species and the rates of molecular reactions could be measured during plasticity induction both in a single cortical neuron and cortical astrocyte using novel imaging techniques [121]. The NMDAR functioning should as well be further studied in wetlab experiments, particularly addressing the type, time courses, and density of presynaptic NMDARs. There is a great demand for new targets for treating neurodevelopmental disorders and diseases. Systematic collection of experimental data on the role of how astrocytic signaling pathways impair synaptic plasticity in developmental brain disorders is crucial. Taken together, all these future experiments will enable deeper insights into the players of long-term plasticity in developing circuits in health and disease by providing data for construction and validation of models.
It is extremely complex to model synaptic plasticity and the underlying biochemical networks in a biologically meaningful way. Despite the challenges, we were able to bring about a combination of experimentally verified neuronal and astrocytic mechanisms and show how they lead to the emergence of spike-timing-dependent long-term plasticity. Our analysis confirms the experimental findings on astrocytes' ability in setting the temporal difference of t-LTD at developing somatosensory L4-L2/3 synapses [15]. Furthermore, we predicted with our in silico synapse model (1) which are the key molecules related to t-LTD, (2) how the molecular reactions depend on the temporal difference of t-LTD, and (3) what are the time courses of molecular interactions. The synapse model can be used to design future wet-lab experiments and, ultimately, to clarify the controversies present in the field. Our study provides both neuronal and neuroglial elements to build sophisticated and biologically relevant large-scale neuronastrocyte network models. With such models bridging different scales, we will expect to link the molecular, synaptic, cellular, and network level dynamics to cognitive phenomena and to assess the roles of astrocytes in higher brain functions, such as learning, memory, decisionmaking, sleep, and, ultimately, consciousness.

Materials and methods
To study the role of astrocytes in modulation of t-LTD, we simulated an L4-L2/3 synapse in somatosensory cortex. We described major biophysical and biochemical mechanisms for the one-compartmental presynaptic L4 spiny stellate cell, two-compartmental (soma and dendrite) postsynaptic L2/3 pyramidal cell, and one-compartmental fine astrocyte process (Fig 1). We employed the following key assumptions to build our initial hypotheses about the testable cellular and subcellular mechanisms: (1) Endocannabinoid 2-AG activates astrocytic CB 1 Rs and triggers Ca 2+ signaling in astrocytes in somatosensory cortex [15,42,97], (2) astrocytic Ca 2+ -dependent glutamate exocytosis, together with a spillover of glutamate from the synaptic cleft, has an effect on presynaptic glutamate release by modifying the release probabilities [36,42,88,99], and (3) the link between the glutamate exocytosis from the astrocyte and the presynaptically released glutamate is the protein phosphatase calcineurin which is activated by the influx of Ca 2+ through the presynaptic NMDARs [13,14]. The model components are described using differential equations and validated against experimental data. We stimulated the synapse model using t-LTD stimulation protocols with a varying temporal difference between pre-and postsynaptic activity [12]. For clarity, only those differential equations that we developed or modified from previously published models are given next. A complete description of the model is given in S1 Appendix.

Presynaptic neuron model
The differential equation for the presynaptic membrane potential can be given as [56] C m;pre dV pre dt ¼ À I CaNHVA;pre À I K;pre À I Na;pre À I L;pre À I Ca;NMDAR;pre À I Na;NMDAR;pre þ I ext;pre ; where C m,pre is the presynaptic membrane capacitance per unit area, I CaNHVA,pre is the current density via Ca NHVA channels, I K,pre is the K + current density, I Na,pre is the Na + current density, I L,pre is the leak current density, I Ca,NMDAR,pre and I Na,NMDAR,pre are the Ca 2+ and Na + current densities via NMDARs, and I ext,pre is the stimulus current injected into the presynaptic neuron per unit area. The presynaptic channels are described by the Hodgkin-Huxley and Goldman-Hodgkin-Katz formalisms [56,57] as explained in S1 Appendix. The differential equations for the gating variables of different currents are given in S1 Appendix [56,57]. Presynaptic Ca 2+ concentrations were elevated by Ca 2+ influxes through presynaptic NMDARs and Ca NHVA channels. The differential equations for the presynaptic Ca 2+ concentration mediated by Ca NHVA channels and by NMDARs are based on a previously published study [123]. The concentration of Ca 2+ mediated by Ca NHVA channels ([Ca 2+ ] CaNHVA,pre ) activates vesicle exocytosis and glutamate release from the presynaptic neuron. The concentration of Ca 2+ mediated by NMDARs ([Ca 2+ ] NMDAR,pre ) activates presynaptic calcineurin [14], and the differential equation for the presynaptic calcineurin concentration ([CaN] pre ) is given in S1 Appendix [60].
Calcineurin has been shown to regulate a specific phase of synaptic vesicle cycling, thus influencing the vesicle release [11,[124][125][126]. We modeled this effect via a signaling pathway linking calcineurin to vesicle release and recycling in the presynaptic terminal with the following differential equation d½X� ac;pre dt ¼ p 1;pre ½CaN� n 2;pre pre K n 2;pre A;pre þ ½CaN� n 2;pre pre X total;pre À ½X� ac;pre where [X] ac,pre is the active concentration and X total,pre is the total concentration of the unspecified protein that affects the vesicle release, p 1,pre is the rate constant, K A,pre is the calcineurin concentration producing half occupation, and n 2,pre is the Hill coefficient. The differential equation for the fraction of releasable presynaptic vesicles (R rel,pre ) was taken from previously published models [62][63][64][65], and the differential equation for the release probability of presynaptic glutamate vesicles was combined and modified from previously published equations [62][63][64][65] and is given as where the fraction (f pre ), which is the active concentration ([X] ac,pre ) divided by the total concentration (X total,pre ) of the protein, affects the probability of presynaptic glutamate release. Through f pre , we modeled the inhibiting role of calcineurin in vesicle exocytosis and glutamate release [13,14] in the presynaptic terminal. Parameters k f,pre , K rel,pre , and n 1,pre describe the facilitation rate constant, Ca 2+ concentration producing half occupation used in calculation of glutamate release, and Hill coefficient, respectively. The presynaptic glutamate release occurs at the first time point t = t j such that [Ca 2+ ] CaNHVA,pre � C thr,pre and less than 10 ms has passed from the previous presynaptic membrane potential crossing 0 mV from negative to positive voltages (V pre � 0, dV pre dt > 0) at that time point t j . The δ function has units of 1 ms . The differential equation for the glutamate concentration in the synaptic cleft was combined and modified from previously published glutamate equations [63][64][65] and glutamateactivated postsynaptic equations related to mGluRs [53] and is given as where k Glu,f,post , k mGluR,f,post , and k mGluR,b,post are the rate constants for the postsynaptic mGluR glutamate uptake, and postsynaptic mGluR glutamate binding and unbinding, respectively. Parameter f Glu,pre represents the spillover of glutamate from the synaptic cleft, and thus the amount 1 − f Glu,pre denotes the part of glutamate in synaptic cleft that activates the postsynaptic receptors.
[mGluR] post and [Glu_mGluR] post denote the concentrations of postsynaptic mGluRs and glutamate-mGluR complex, respectively. Parameters G pre , N pre , k Glu,pre , N A , and V syncleft denote the number of glutamate per presynaptic vesicle, number of readily releasable presynaptic vesicles, scaling factor to convert from units M to μM, Avogadro's constant, and volume of synaptic cleft, respectively. We modeled presynaptic GluN2C/D-containing NMDARs [58,59], because experimental studies have reported that presynaptic GluN2C/D-containing NMDARs are required for t-LTD at L4 to L2/3 synapses [9,12,15,[89][90][91]. On the other hand, postsynaptic GluN2B-containing NMDARs are necessary for t-LTP at L4-L2/3 and L2/3-L2/3 synapses, and postsynaptic GluN2A-containing NMDARs are required in t-LTD at L2/3-L2/3 synapses [9,12,15,[89][90][91]127]. In our synapse model, presynaptic NMDARs are activated by the glutamate in the extrasynaptic space and the spillover of glutamate from the synaptic cleft. The differential equations for the presynaptic variables can be seen in S1 Appendix. Intermediate variables, parameter values, and initial values needed to solve the presynaptic neuron model are also given in S1 Appendix.

Postsynaptic neuron model
The postsynaptic neuron model had two compartments, a soma and a dendrite, modified from a previously published study [66]. The differential equations for the membrane potentials of these two compartments are where C m,post is the membrane capacitance per unit area, I KDR,soma,post is the somatic K DR current density, I Na,soma,post and I Na,dend,post are the somatic and dendritic Na + current densities, I NaP,soma,post is the somatic Na P current density, I L,soma,post and I L,dend,post are the somatic and dendritic leak current densities, I coupl,soma,post and I coupl,dend,post are the somatic and dendritic coupling terms, I ext,post is the current injected into the soma per unit area, I KA,dend,post is the dendritic K A current density, I CaLHVA,dend,post and I CaLLVA,dend,post are the dendritic Ca LHVA and Ca LLVA current densities, and I AMPAR,post and I Ca,NMDAR,post are the synaptic current densities via AMPARs and NMDARs in the dendrite. Hodgkin-Huxley formalism was used to describe the behavior of ionic currents [54,[67][68][69] as explained in S1 Appendix. The differential equation for the fraction of postsynaptic AMPARs in open state can be modified from previously published study [70] and is given as where α AMPAR,post and β AMPAR,post describe the rate constants of opening and closing postsynaptic AMPARs, respectively, and we can similarly write the equation for NMDARs. Other differential equations for the gating variables of different currents related to the membrane potential in the dendrite and in the soma as well as for the IP 3 R inactivation gating variable are given in S1 Appendix [54, 67-69, 72, 73]. The biochemical mechanisms related to mGluRs, the activation of the G-protein signaling cascade, and the production of endocannabinoids are included in the synapse model to study the effect of endocannabinoids on the adjacent astrocyte. The differential equations starting from the mGluR activation to the endocannabinoid 2-AG release were based on previously published models [53,71]. Glutamate in the synaptic cleft binds to postsynaptic mGluRs and induces dissociation of the G protein α subunit bound with guanosine-5'-triphosphate (GαGTP) from the mGluR-bound G protein with β and γ subunits (Gβγ). Calcium can bind to PLC, and, in addition, GαGTP can enhance its activity. The postsynaptic Ca 2+ equation was based on previously published models [53,[72][73][74]. Active PLC produces IP 3 and DAG from phosphatidylinositol 4,5-bisphosphate (PIP 2 ). After Ca 2+ binds to DAG lipase, the complex binds to DAG and catalyzes 2-AG synthesis. The differential equations, intermediate variables, parameter values, and initial values needed to solve the postsynaptic neuron model are given in S1 Appendix.

Astrocyte model
More and more evidence about the complexity of astrocyte processes is accumulating in vivo. Several studies have revealed Ca 2+ activity [27,44,101,102,[128][129][130][131] and complex molecular and biochemical mechanisms [40,[42][43][44][45][46][47][48][49] in the main and fine processes of astrocytes. Currently it is not clear how these signals and the underlying subcellular mechanisms are linked. These recent findings encouraged us to test the dynamical capacity and the role of astrocytic Ca 2+ -dependent glutamate exocytosis in the endocannabinoid-mediated signaling from the postsynaptic terminal to the presynaptic terminal to modulate synaptic functions. The fine astrocyte process in our model is assumed to contain IP 3 R-mediated Ca 2+ -dependent glutamate exocytosis [42,99,100], similar to earlier published models on astrocyte-neuron interactions (for a summary of models, see [78,79]). We modeled IP 3 R-mediated Ca 2+ -dependent glutamate exocytosis as a generic glutamate exocytosis, not as a biochemically detailed vesicular release due to lacking molecular details [132,133]. Contradictory evidence has been presented on the involvement of astrocytic IP 3 Rs on synaptic plasticity in hippocampal slices, using transgenic mice that lack the commonly expressed subtype 2 of the IP 3 Rs in astrocytes [103,134]. It is possible that other IP 3 R subtypes exist in astrocyte processes [40], indeed knocking out the subtype 2 of IP 3 Rs has been shown to abolish all Ca 2+ signals in astrocyte soma but only about half in the astrocyte processes [44]. These recent results about the diversity of IP 3 R subtypes in astrocyte processes [40] provide additional justification for our assumption to further examine the significance of the kinetics of IP 3 R-mediated Ca 2+ -dependent glutamate exocytosis in t-LTD in somatosensory cortex [15,42]. Other mechanisms that couple the endocannabinoids to astrocyte Ca 2+ signaling may coexist in somatosensory cortex in different phases of development and should be studied in the future, but our study focused on the most modeled and tested Ca 2+ -dependent mechanism in astrocytes.
We modeled Ca 2+ and IP 3 concentrations, and the gating variable for IP 3 R inactivation in the astrocyte based on previously published models [72,73,75,80]. The differential equation for the astrocytic IP 3 concentration was modified to be where IP ? 3;astro , τ IP3,astro , r IP3,astro , [2-AG] post , and AG ? post denote the resting concentration of IP 3 , time constant for IP 3 degradation, rate constant of IP 3 production, concentration of the endocannabinoid 2-AG released from the postsynaptic neuron, and resting concentration of 2-AG, respectively. The differential equation for the fraction of releasable glutamate resources in the astrocyte was taken from previously published models [62][63][64][65] and the glutamate concentration in the extrasynaptic space was also taken from the previously published model [64,65]. The differential equations, intermediate variables, parameter values, and initial values needed to solve the astrocyte model are given in S1 Appendix.

Stimulation protocols
The following protocols were used [12]: 1. The stimulation protocol before t-LTD induction consisted of five 10 ms long presynaptic stimuli at a frequency of 0.2 Hz and with an amplitude of 10 mA cm 2 keeping the fraction of presynaptic glutamate release inhibition (f pre ) as constant zero (Fig 2A).
2. The t-LTD induction protocol consisted of a 10 ms long postsynaptic stimulus with an amplitude of 25 mA cm 2 occurring between 10 ms and 200 ms before a 10 ms long presynaptic stimulus with an amplitude of 10 mA cm 2 and the post-pre pairing was repeated 100 times at a frequency of 0.2 Hz (Fig 2H). Thus, the temporal difference (ΔT) between the pre-and postsynaptic stimulus in this study had negative values meaning that the postsynaptic stimulus occurred before the presynaptic stimulus (ΔT had values between −10 ms and −200 ms). The initial value of f pre in these simulations was zero, but it increased during the stimulation protocol to above zero depending on ΔT used in the protocol.
3. The stimulation protocol after t-LTD induction consisted of five 10 ms long presynaptic stimuli at a frequency of 0.2 Hz and with an amplitude of 10 mA cm 2 (Fig 2O) keeping f pre as constant value that is the final simulation value obtained from the simulation with the t-LTD induction protocol. Thus, the fraction f pre had different constant values depending on ΔT used in the t-LTD induction protocol (Fig 2H).

Data analysis
Data was analyzed in MATLAB 1 . We calculated the amplitude of the postsynaptic EPSP in the soma during the first presynaptic stimulus in the stimulation protocol before t-LTD induction and also for every ΔT in the stimulation protocol after t-LTD induction. To obtain the EPSP percentage, we normalized for every ΔT the EPSP amplitude occurring after t-LTD induction by the EPSP amplitude occurring before t-LTD induction, and multiplied them by 100%. Furthermore, we calculated the ΔEPSP percentage for every ΔT as the percentage change between the somatic EPSP amplitude occurring before t-LTD induction and the somatic EPSP amplitude occurring after t-LTD induction. Thus, the ΔEPSP percentage was obtained by subtracting the somatic EPSP amplitude occurring before t-LTD induction from the somatic EPSP amplitude occurring after t-LTD induction, dividing this change with the somatic EPSP amplitude occurring before t-LTD induction, and finally multiplying by 100%. The shorter way to calculate the ΔEPSP percentage was to subtract 100% from the EPSP percentage. The experimental ΔEPSP data for comparison was calculated from the normalized EPSP slopes and obtained from the literature [12].
We calculated the Ca 2+ peak values, peak times, and oscillation frequencies during the t-LTD induction protocol for every ΔT. We divided the 500 s long post-pre pairing simulation data into 5 s long sweeps, where each sweep represents one post-pre pairing (0.2 Hz stimulus). We reorganized the peak times of Ca 2+ oscillations into the 5 s long sweeps. Then we reorganized the number of peaks occurring in ten 0.5 s long bins. We obtained the probability of astrocytic Ca 2+ peaks by normalizing the reorganized data by the total number of peaks [15].

Simulation details and code
Simulation code was written in Python 3.7. The code is available in the ModelDB [135] (http:// modeldb.yale.edu/266819) and in the author's GitHub page (https://github.com/ TiinaManninen/synapsemodel). State variables were updated using the forward Euler method with 0.05 ms step size.