Modeling Working Memory in a Spiking Neuron Network Accompanied by Astrocytes

We propose a novel biologically plausible computational model of working memory (WM) implemented by a spiking neuron network (SNN) interacting with a network of astrocytes. The SNN is modeled by synaptically coupled Izhikevich neurons with a non-specific architecture connection topology. Astrocytes generating calcium signals are connected by local gap junction diffusive couplings and interact with neurons via chemicals diffused in the extracellular space. Calcium elevations occur in response to the increased concentration of the neurotransmitter released by spiking neurons when a group of them fire coherently. In turn, gliotransmitters are released by activated astrocytes modulating the strength of the synaptic connections in the corresponding neuronal group. Input information is encoded as two-dimensional patterns of short applied current pulses stimulating neurons. The output is taken from frequencies of transient discharges of corresponding neurons. We show how a set of information patterns with quite significant overlapping areas can be uploaded into the neuron-astrocyte network and stored for several seconds. Information retrieval is organized by the application of a cue pattern representing one from the memory set distorted by noise. We found that successful retrieval with the level of the correlation between the recalled pattern and ideal pattern exceeding 90% is possible for the multi-item WM task. Having analyzed the dynamical mechanism of WM formation, we discovered that astrocytes operating at a time scale of a dozen of seconds can successfully store traces of neuronal activations corresponding to information patterns. In the retrieval stage, the astrocytic network selectively modulates synaptic connections in the SNN leading to successful recall. Information and dynamical characteristics of the proposed WM model agrees with classical concepts and other WM models.


Introduction
In neuroscience, the understanding of the functional role of astrocytes in CNS is still open to debate (Savtchouk and Volterra, 2018), but now there is much evidence demonstrating the involvement of astrocytes in local synaptic plasticity and coordination of network activity (Durkee and Araque, 2019), and as a result in information processing and memory encoding (Santello et al., 2019).Astrocytes sense synaptic activity and respond to it with the transient elevation of the intracellular Ca 2+ concentration (lasting from hundreds of a millisecond to a dozen of seconds).Such Ca 2+ signals in astrocytes have been observed in different brain regions and also in the cortex, appearing there in response to the mechanic sensory stimulation (Stobart et al., 2018;Wang et al., 2006a;Takata et al., 2011) and visual sensory stimulation (Chen et al., 2012;Schummers et al., 2008;Perea et al., 2014).Ca 2+ activation can trigger the release of gliotransmitters from astrocyte, which in turn affect the dynamics of presynaptic and postsynaptic terminals resulting in modulations of synaptic transmission (Araque et al., 2014).The gliotransmittermediated synaptic modulation lasts from a dozen of seconds (Jourdain et al., 2007;Perea et al., 2014) to a dozen of minutes (Perea and Araque, 2007;Navarrete et al., 2012;Stellwagen and Malenka, 2006) contributing to both short-and long-term synaptic plasticity.Obviously, there is a qualitative coincidence of time scales of astrocyte-mediated synaptic modulation with the working memory (WM) timings during decision making.Based on this and the other following facts of astrocyte participation in neuronal signalling, we hypothesized that the astrocytes may be involved in the WM formation.In particular, recent in vivo studies have shown the participation of astrocytes in synchronization of certain cortical network activity (Takata et al., 2011;Chen et al., 2012;Perea et al., 2014;Paukert et al., 2014), cognitive functions and behaviour (Sardinha et al., 2017;Poskanzer and Yuste, 2016).Experimental evidence shows that astrocyte pathology in medium PFC impairs the WM and learning functions (Lima et al., 2014), increasing of astrocyte density enhances short-term memory performance (Luca et al., 2020), recognition memory performance and disruption of WM depend on the gliotransmitter release from astrocyte in the hippocampus (Robin et al., 2018;Han et al., 2012).Despite these numerous experimental insights of the contribution of astrocytes to synaptic modulations in neuronal signalling, the possible role of astrocyte in information processing and learning is still a subject of discussion (Kastanenka et al., 2019;Kanakov et al., 2019).
Considering the significance of WM processes, the challenge of finding alternative mechanisms and the experimental evidence of the astrocytic role in information processing in CNS, it is interesting to study astrocyte-induced modulation of synaptic transmission in the WM organisation.Specifically, we assume that the NMDAR-mediated potentiation of excitatory synapses induced by the D-serine released from astrocyte in PFC (Fossat et al., 2011;Takata et al., 2011;Chen et al., 2012) plays an essential role in WM.To test this hypothesis, we developed a novel neuron-astrocyte network model for visual WM to reflect experimental data on the structure, connectivity, and neurophysiology of the neuron-astrocytic interaction in underlying cortical tissue.We focused on the implementation of a multi-item WM task in DMS framework representing classical neuropsychological paradigm (Miller et al., 1996), which was previously studied using a recurrent neural network (Brunel and Wang, 2001;Amit, 2003;Amit et al., 2013;Fiebig and Lansner, 2016).In our model, memory is associated with item-specific patterns of astrocyteinduced enhancement of excitatory synaptic transmission.We show how the biologically relevant neuron-astrocyte network model implements loading, storage and cued retrieval of multiple items with significant overlapping.The memory items are encoded in neuronal populations in the form of discrete high-frequency bursts rather than persistent spiking. 2 In the following, we review some related works (Section 2), describe the proposed model and methods in detail (Section 3), present the results (Section 4), and finally, we conclude this work in Section 5.

Related work
The concept of WM proposes the ability to temporarily store and process information in goaldirected behaviour.WM is crucial in the generation of higher cognitive functions for both humans and other animals (Baddeley, 1986;Conway et al., 2003;Baddeley, 2012).In primates, visual WM has been studied in delay tasks, such as delayed matching to sample (DMS), which require memory to be held during a brief delay period lasting for several seconds (Miller et al., 1996).Recordings in the monkeys' prefrontal cortices (PFCs) during the delay task showed that some neurons displayed persistent and stimulus-specific delay-period activity (Fuster and Alexander, 1971;Funahashi et al., 1989;Barak et al., 2010;Shafi et al., 2007;Funahashi, 2017).Delay persistent activity is considered the neural correlate of WM (Goldman-Rakic, 1995;Constantinidis et al., 2018).
The classical theoretical memory models suggest that an information item can be stored with sustained neural activity which emerges via activation of stable activity patterns in the network (e.g.attractors) (Hopfield, 1982;DJ, 1995;Wang, 2001;Wimmer et al., 2014) recently reviewed by Zylberberg and Strowbridge (2017); Chaudhuri and Fiete (2016).These WM models propose that the generation of the persistent activity can be the result of an intrinsic property of the neurons (including the generation of the bistability mediated by the voltage-gated inward currents (Kass and Mintz, 2005) and Ca 2+ -triggered long-term changes in neuronal excitability (Fransén et al., 2006)) and can be induced by the connectivity within the neural circuit with feed-forward (Ganguli and Latham, 2009;Goldman, 2009) or recurrent (Kilpatrick et al., 2013;Koulakov et al., 2002) architecture.In such models, memory recall is impossible from a silent inactive state.For many WM models of persistent activity based on recurrent connectivity, small deviations in the network structure destroy the persistence.Moreover, a spiking form of storage information is energetically unfavourable because of the high metabolic value of action potentials (Attwell and Laughlin, 2001).
In theoretical studies a concept of oscillatory sub-cycles storing 7 ± 2 in oscillatory neuronal networks were proposed by Lisman and Idiart (1995).Other models employ oscillatory activity of spiking neuron networks with afterdepolarization to memorize set of information patterns at different phases of rhythmic oscillations (Klinshov and Nekorkin, 2008;Borisyuk et al., 2013).
Recently, the persistent activity hypothesis has been undergoing critical reviews (Lundqvist et al., 2018) based on the experimental findings in the rodents, rats, and primates showing that the robust persistent activity does not last for the entire delay period, but rather sequential neuronal firing is observed during the delay period suggesting that the PFC neural network may support WM based on dynamically-changing neuronal activity (Fujisawa et al., 2008;Runyan et al., 2017;Park et al., 2019;Ozdemir et al., 2020;Lundqvist et al., 2016).Despite the considerable progress which has been made in identifying the neurophysiological mechanisms contributing to WM in mammals (D'Esposito and Postle, 2015;Zylberberg and Strowbridge, 2017), the ongoing debate focuses on the generation mechanisms of the delay period activity that appears to underlie the WM (Constantinidis et al., 2018;Sreenivasan and D'Esposito, 2019).Currently, one of the recognised experimentally based hypotheses of the WM mechanism underlining the delay activity (not necessarily persistent) is the synaptic plasticity in the PFC (Tsodyks and Markram, 1997;Wang et al., 2006b;Hempel et al., 2000;Erickson et al., 2010).Synaptic plasticity implies a rapid regulation of the strengths of individual synapses in response to specific patterns of correlated synaptic activity and contributes to the activity-dependent refinement of neural circuitry.Following these findings, alternative synaptic-based WM models have been proposed (Mongillo et al., 2008;Manohar et al., 2019;Barak and Tsodyks, 2014;Koutsikou et al., 2018).In these models, memory items are stored by stimulus-specific patterns of synaptic facilitation in neuronal circuit.Synaptic plasticity does not require neurons to show a persistent activity for the entire period of the memory task, which results in a robust and more metabolically efficient mechanism.Some synaptic WM models based on short-term non-associative synaptic facilitation (Mongillo et al., 2008;Lundqvist et al., 2011;Mi et al., 2017) allow reading out and refreshing existing representations maintained in the synaptic structure.Others have proposed fast Hebbian activitydependent synaptic plasticity (Fiebig and Lansner, 2016;Sandberg et al., 2003) for encoding and maintenance of novel associations.
There are a few attempts to investigate theoretically the role of astrocyte-induced modulation of synaptic transmission in memory formation.Shen and Wilde (2007) demonstrate one of the first results of simulating the coupling of Hopfield neural network, astrocyte, and cerebrovascular activity.Although there remained a far stretch to a biophysical model, the resulted suggested that a modification of the synapse strengths allows the neuronal firing and the cerebrovascular flow to be compatible on a meso-scale; with astrocyte signalling added, limit cycles exist in the coupled networks.Tewari and Parpura (2013) and Wade et al. (2011) study how bidirectional coupling between astrocytes and neurons in small neuron-astrocyte ensembles mediates learning and dynamic coordination in the brain.Recent interesting theoretical study proposes a selfrepairing spiking astrocyte-neural network combined with novel learning rule based on the spiketiming-dependent plasticity and Bienenstock, Cooper, and Munro learning rule (Liu et al., 2019).

Neuron-astrocyte network model
Even though the balance of inhibition and excitation play was shown to play an important role in WM stabilization and to influence the WM capacity (Barak and Tsodyks, 2014), we focus on the properties of astrocyte-induced modulation of excitatory synaptic transmission in the PFC.We take spiking neuronal network with dimension W × H consisting of synaptically coupled excitatory neurons formed by the Izhikevich model (Izhikevich, 2003).Neurons in the network are connected randomly with the connection length determined by the exponential distribution.It has been experimentally estimated that there is a little overlap in the spatial territories occupied by individual astrocytes in the cortex (Halassa et al., 2007).An individual cortical astrocyte contacts on average 4-8 neuronal somata and 300-600 neuronal dendrites (Halassa et al., 2007).A cortical astrocyte has a "bushy" ramified structure of the fine perisynaptic processes, which cover most of neuronal membranes within their reach (Allen and Eroglu, 2017).This allows the astrocyte to integrate and coordinate a unique volume of synaptic activity.Following the experimental data, astrocytic network compartment of our model is organised as a two-dimensional square lattice with only nearest-neighbour connectivity.Each astrocyte interacts with the neuronal ensemble of N a neurons with small overlapping.We consider the bidirectional communication between neuronal and astrocytic networks.The scheme of the network topology is shown in Fig. 1.
Model equations are integrated using the Runge-Kutta fourth-order method with a fixed time step, ∆t = 0.1 ms.A detailed listing of model parameters and values can be found in Tables 1 (neural network model), 2 (astrocytic network parameters), 3 (neuron-astrocytic inter-  Blue lines show connections between elements in each layer.We consider the bidirectional interaction between the neuronal and astrocytic layers.Each astrocyte is interconnected with a neuronal ensemble of N a = 16 neurons with dimensions 4 × 4 with overlapping in one row (red lines).The input signal is fed to the neural network.

Neuronal network
Dynamics of the membrane potential of a single neuron is described by the Izhikevich model (Izhikevich, 2003): dU (i, j) dt = a(bV (i, j) − U (i, j) ); (1) with the auxiliary after-spike resetting if V (i, j) ≥ 30 mV, then where the superscripts (i = 1, . . ., 79, j = 1, . . ., 79) correspond to a neuronal index, the transmembrane potential V is given in mV and time t in ms.The applied currents I (i, j) app simulate the input signal.Neurons receive a number of synaptic currents from other presynaptic neurons in the network, N (i, j)  in , which are summed at the membrane according to the following equation Parameter g (i, j) syn describes the synaptic weight, The astrocyte (m, n) modulates the synaptic currents of the neuron (i, j).The variable ν Ca introduces the astrocyte-induced modulation of synaptic strength and will be discussed below.Here, purely for illustration of astrocyte effects, we did not include intrinsic short-term synaptic plasticity in the model.The synaptic reversal potential for excitatory synapses is taken with E syn = 0. V pre denotes the membrane potential of the presynaptic neuron.For simplicity, we neglect the axonal and synaptic delays.
The architecture of synaptic connections between neurons is non-specific (e.g.random) with the following parameters.The number of output connections per each neuron is fixed at N out = 40.Each neuron innervates N out local postsynaptic targets, the distance R to which is determined according to the exponential distribution: (4)

Astrocytic network
In the model we try to implement the biologically plausible organisation of astrocytic network and neuron-astrocyte interaction.The astrocytic network is configured in the form of a two-dimensional square lattice with dimension M × N. Cortical astrocytes are coupled via CX43 gap junctions mostly permeable to inositol 1,4,5-trisphosphate (IP 3 ).Hence, in the model, we consider local diffusive coupling.Besides, each astrocyte is interconnected with neuronal ensemble of N a neurons.It was experimentally shown that the sensory stimulation evokes fast intracellular Ca 2+ signals in fine processes of cortical astrocyte in response to local synaptic activity in neuronal circuit (Stobart et al., 2018;Takata et al., 2011;Wang et al., 2006a).Multiple rapid, spatially restricted Ca 2+ events in the astrocytic process are induced by intense neuronal firing.Local events are spatially and temporally integrated by the astrocytic cell, which results in a global, long lasting Ca 2+ event.In turn, this event induces the release of gliotransmitters affecting synaptic transmission in the local territory of individual astrocytes (Henneberger et al., 2010;Bekar et al., 2008;Araque et al., 2014).For simplicity, we did not model detailed process of spatial-temporal integration of the rapid Ca 2+ signals in the morphological structure of astrocyte modelled earlier by Gordleeva et al. (2019Gordleeva et al. ( , 2018) ) and Wu et al. (2018).Here we employ a mean-field approach to describe the emergence of a global Ca 2+ signal and its impact on synchronisation of neuronal ensemble controlled by a certain astrocyte.
As pyramidal neurons generate the spike, glutamate is released from the presynaptic terminal into the synaptic cleft (Fig. 2b,d).Amount of glutamate, G, that diffused from synaptic cleft and reached the astrocytic process can be described by the following equation (Gordleeva et al., 2012;Pankratova et al., 2019): here α glu is the glutamate clearance constant, kglu is the efficacy of the release, Θ denotes the Heaviside step function and V (i, j) is the membrane potential of the corresponding presynaptic neuron (i, j).Binding of glutamate to metabotropic glutamate receptors (mGluR) on the astrocytic membrane, which is located close to the synapse, triggers the production of IP 3 in the astrocyte (Fig. 2e).We use the approaches from earlier studies to describe the dynamics of the intracellular concentration of IP 3 in the astrocyte (Nadkarni and Jung, 2003;Ullah et al., 2006): with m = 1, . . ., 26, n = 1, . . ., 26.Parameter IP * 3 denotes the steady state concentration of the IP 3 and J PLCδ describes the IP 3 production by phospholipase Cδ (PLCδ) (Ullah et al. (2006)): The variable J glu describes the glutamate-induced production of the IP 3 in response to neuronal activity and is modelled as a rectangular-shaped pulse with amplitude A glu µM and duration t glu ms: here t 0 denotes the periods when the total level of glutamate in all synapses associated with this astrocyte reaches a threshold: here we use the parameter G thr = 0.7.[] denotes the Iverson bracket.F act is the fraction of synchronously spiking neurons of the neuronal ensemble corresponding to the astrocyte.For the emergence of the calcium elevation F act = 0.5 is required.In other words, the production term, J glu , is activated when correlated activity in the neuronal ensemble reaches a certain level of coherence.
Increase of IP 3 concentration in the astrocyte induces the release of Ca 2+ from internal stores, mostly from the endoplasmic reticulum (ER), to cytosol.For simplified description of the biophysical mechanism underlying the calcium dynamics in astrocytes, we use the Ullah model (Ullah et al., 2006).Changes of the intracellular Ca 2+ concentration, Ca, are described by the following equations: where h is the fraction of the activated IP 3 receptors (IP 3 Rs) on the ER surface.Flux J ER is Ca 2+ flux from the ER to the cytosol through IP 3 Rs, J pump is the flux pumped Ca 2+ back into ER via the sarco/ER Ca 2+ -ATPase (SERCA), J leak is the leakage flux from the ER to the cytosol.Fluxes J in and J out describe calcium exchange with extracellular space.The fluxes are expressed as follows: Biophysical meaning of all parameters in Eqs. ( 6), ( 7),( 10),( 11) and their values determined experimentally can be found in (Li and Rinzel, 1994;Ullah et al., 2006) and in Table 2. Cortical astrocytes are coupled by CX43 gap junctions (Nimmerjahn et al., 2004).Thus, diffusion of active chemicals becomes possible between the neighbouring astrocytes.Currents diff Ca and diff IP 3 describe the diffusion of Ca 2+ ions and IP 3 molecules via gap junctions between the astrocytes in the network and can be expressed as follows: where parameters d Ca and d IP 3 describe the Ca 2+ and IP 3 diffusion rates, respectively.Following experimental data, we assume that CX43 is less permeable to Ca 2+ than to d IP 3 .(∆Ca) (m,n) and (∆IP 3 ) (m,n) are the discrete Laplace operators: Equations ( 6)-( 8),( 10)-( 12) predict that the synchronized activity in the neuronal ensemble triggers the astrocytic Ca 2+ signals, and in the absence of neuronal stimulus in the astrocytic network steady state Ca 2+ concentration is maintained.Next, we account for the effect of enhancement of excitatory postsynaptic currents (ePS C) generation through modulation of postsynaptic NMDARs by D-serine released from the astrocyte IP 3 diffusion rate 0.1 s −1 (Bergersen et al., 2011;Henneberger et al., 2010).In the model, we propose that global events of Ca 2+ elevation in astrocyte result in D-serine release, which can modulate the synaptic strength of all synapses corresponding to the morphological territory of this astrocyte.For simplicity, the relationship between the astrocyte Ca 2+ concentration and synaptic weight of the affected synapses g syn , is described as follows: where the parameter ν * Ca denotes the strength of the astrocyte-induced modulation of the synaptic weight, Θ(x) is the Heaviside step-function.The feedback from the astrocytes to the neurons is activated when the astrocytic Ca 2+ concentration is larger than Ca thr and the fraction of synchronously spiking neurons of neuronal ensemble corresponding to the astrocyte F astro during the time period of τ syn = 10 ms.The duration of the feedback is fixed and is equal to τ astro = 250 ms.

Stimulation protocol
The term I app in Eq. ( 1) represents specific and non-specific external inputs.A non-specific noisy input simulates the input signals from networks of other brain areas and is applied continuously to all neurons in the form of independent Poisson pulse trains of a certain rate, f bg , with Average size of sample is 1078 neurons (18% of network) stimulated by the specific input, with average 35.2% overlapping in the population.For visual representation of samples, we take binary images of numerals (0,1,2,3,4,..) with size W × H pixels, where each pixel corresponds to a neuron in the neuronal layer.Neurons corresponding to the shape of the numerals receive rectangular excitatory pulse with length t stim and amplitude A stim .The shape of each sample was spatially distorted by 5% random noise such as "salt and pepper noise".Then a transient inputs were applied to simulate the nonmatching test items and the cue (length t test , and amplitude A test ).In the cued recall for simulating the loss in saliency, we applied shorter input with lower amplitude and higher level (20%) of random noise.10

Memory performance metrics
To measure memory performance of the system we count the correlation of recalled pattern with ideal item in the following way: here w = 10 frames = 1 ms, P -a set of pixels belonging ideal pattern, W, H -network dimension, thr -spike threshold, I -indicator function, T P -a set of frames in the tracking range of pattern P. In a sense, this correlation metric can be associated with 1 − d averaged between pattern and background, where d is the Hamming distance.

Results
Let us show how the neuron-astrocyte network model exhibits memory formation.First, we will consider a simple single-item memory task illustrating information loading, storage, and retrieval.Next, we will demonstrate how the network can be successfully trained to memorize and recall several patterns with significant overlaps.Finally, we will analyse model performance metrics, capacity, and characteristic of pattern remembering on different parameters.

Single-item WM
First, we will test the neuron-astrocyte network in the most common experimental paradigm of WM studies -the DMS task.This task requires a single item to be held in memory during a brief delay period.Before specific stimulation, neural network demonstrates irregular, low-rate background activity (see activity beginning in Fig 3).At the 500 ms mark, we load an item by applying transient external input to the corresponding neuronal population for 200 ms (Figs.2a,  3).During training, each astrocyte tracks the activity of neuronal subnetwork associated with it.As soon as the extracellular concentration of glutamate (Figs.2b,d) and correlated firing in neurons achieve a certain level, which satisfies the condition Eq. ( 9), Ca 2+ concentration in matching astrocytes elevates (Fig. 2e).In accordance with the experimental data (Bindocci et al., 2017), we tuned the model parameters in such way that the onset of the calcium elevation in astrocytes induced by synchronous neuronal discharge has a delay of ≤ 2 sec.Following the increased firing in the stimulus-specific part of the neuronal network upon reaching a threshold of 0.15 µM the astrocyte releases gliotransmitters modulating the synaptic strengths in corresponding locations (Fig. 2).The calcium pulse in astrocyte lasts for several seconds.Its duration determines the length of the delay interval in the DMS task, during which the item is maintained in the memory.After the training stimulus ends, we test maintenance of the single-item memory by applying two nonmatching items and cue item with t test durations and 250 ms inter-item interval (Figs. 2, 3).Because the astrocytic feedback also depends on the activity of neuronal subnetwork, the model responds differently to the applied items.A short presentation of the cue to the neural network evokes the astrocytic-induced increase in the synaptic strength between stimulus-specific neurons and results in a local spatial synchronization in the whole stimulus-specific neuronal population (Fig. 2a with comparison to Fig. 2c).Similar to experimental data (Miller et al., 1996), delay activity in our model is sample-selective.We observe that the pattern-specific firing rate in the neuronal network increases and is equal to 270 Hz in comparison with response to non-specific stimulus (80 Hz) (Fig. 3,b).Such a high frequency is determined by the choice of fast-spiking neuron model (Izhikevich, 2003).The firing rates in simulations with regular spiking neuron model (Izhikevich, 2003) are almost 10 times lower: 30 Hz for stimulus-specific and 4.5 Hz for non-specific stimulus.The elevation of the frequency in the stimulus-specific neuronal population can continue after the end of the cue, which is determined by duration of the astrocyte-induced enhancement of the synaptic weight.
For a visual representation of memory formation, we follow space-time distribution of sampleselective delay activity.Figure (4) illustrates the spatial distribution of activity in neuronal and astrocytic layers at the different moments of training and cued recall for the same single-item memory task as presented in the Figure (3).Training induces the emergence of synchronized calcium activity of spatially clustered astrocytes (Fig. 4c).Note that locally synchronised astrocytes have been found in the neocortex and hippocampus in situ and in vivo (Sasaki et al., 2011;Takata and Hirase, 2008).Such calcium activity correlated in time and space can lead to the spatial-temporal synchronization in the neuronal network (Araque et al., 2014).This mechanism of neuron-astrocyte network interaction underlies the sample-selectivity and pattern retrieval in the model.Note that 20 % noisy cue item (Fig. 4d) can be identified and cleared from noise by the neuronal network due to the astrocyte-induced feedback (Fig. 4e).In other words, the spiking neuronal network accompanied by astrocytes can filter a cue pattern distorted by noise.Video of the single-item memory encoding and cued recall in the neuron-astrocyte network can be found in the supplementary material.

Multi-item WM
Next, we consider the multi-item WM formation.In this case, we loaded four items, images of numerals 0,1,2,3 in the following way.The images were loaded consequently by applying external inputs of t stim durations with 100 ms inter-item intervals (see Figs. 5,6a).Due to the coincidence of different stimulus-specific neuronal populations in space, the spatial calcium patterns in astrocytic layers for different items overlap significantly (Fig. 7h).After a 700 ms training stimulus was applied, we tested the maintenance of the memory by applying matching and nonmatching items of t test durations and 250 ms inter-item intervals (Figs. 5, 6a).The images were distorted by 20% noise.The astrocyte-mediated feedback modulating coherent neuronal activity provided the selectivity of the model response.The system remembered the correct image.Thus, we observed that all items were successfully filtered only in the cued recall.Needless to say that firing rate increases significantly in the cued recall due to the selective increase of synaptic strengths (Figs.6b,c).To evaluate the performance of the neuron-astrocyte WM, we used as metric the correlation between recalled item and the ideal item during the multiitem WM task (see section Memory performance metrics) (Fig. 6d).It is important to note, that during multi-item remembering spurious correlations never dominate in that sense accuracy of our system is always equal to 100%.There was an increase in correlation with the target image and no attraction to the wrong image or chimeras.Maximal correlation reached 95% in training and 93% in testing sets on average for 4 samples.Video of the multi-item WM in the neuronastrocyte network can be found in the supplementary material 2.
To characterise the quality of memory formation in the model, we examined the dependencies of correlation of retrieval pattern in cued recall on variable parameters of input patterns, astrocytic, synaptic, and network structure (Figs.8,9).First, we investigated its dependence on noise parameters.The dependence of correlation of recalled pattern on the noise level in training and test experiments is shown in Figure (8a).Specifically, the correlation difference between the recalled pattern and noisy input is presented.In other words, the model can improve test images  The morpho-functional structure of connections between neurons and astrocytes can affect pattern retrieval in the model (Figs.8b,c).The key parameters determining this structure are the fraction of synchronously spiking neurons of neuronal ensemble corresponding to the astrocyte required for the emergence of the calcium elevation in the astrocyte, F act , and the fraction of synchronously spiking neurons of neuronal ensemble corresponding to the astrocyte required for the emergence of astrocyte-induced enhancement of synaptic transmission, F astro .Here we do not account for various spatiotemporal properties of gliotransmitter release and astrocyte Ca 2+ signals evoked by different levels of neuronal activity (Araque et al., 2014).We assume that simultaneous activation of synapses induce multiple Ca 2+ events at different processes of astrocyte, which are spatially and temporally integrated and result in the generation of a global, long lasting Ca 2+ elevation (Bindocci et al., 2017) that can affect synaptic transmission in the territory of individual astrocytes.For this purpose, parameters F act and F astro estimate the correlation level of synapses activity in the model.The optimal range of F astro for correlation of recalled pattern is [40-60]% (Fig. 8b).Smaller values of the parameter, F astro , lead to the effect of astrocyteinduced synchronization initiated even by non-stimulus specific uncorrelated noise activity in a small ensemble of neurons.On the contrary, the use of larger F astro values implies that a highly correlated activity of almost all neurons located in the territory of a given astrocyte is required for existence of astrocytic modulation of synapses.Hence, the neuron-astrocytic network can not perform a correct recall of noisy cue.Another point is that the Figure 8b was obtained for training set with a low 5% noise level and did not reveal the dependence of the correlation of recalled pattern on the parameter F act .We studied the influence of the parameter F act on the correlation in the simulations with different noise in training samples (Fig. 8c).For lower noise level in training samples, the network memorises items regardless of the value of the parameter, F act .
Increasing the level of noise in training samples for small values of the parameter, F act , leads to the Ca 2+ elevation in randomly distributed astrocytes first, and then in the whole astrocytic layer.Nevertheless, such nonstimulus-specific astrocytic activation can result in high correlation of recalled pattern because of the moderate noise level in the cue and optimally chosen value of the parameter F astro .On the contrary, for larger value of the F act > 85% Ca 2+ signal in astrocyte can be evoked only by relatively "clean" sample with a small percentage of noise < 5%, therefore increasing the level of noise in training samples results in poor correlation (Fig. 8c).We found that range of F act = 80 − 85% was optimal for performing the WM tasks by neuron-astrocyte network.In this range astrocyte activations were stimulus-specific and the astrocyte layer could memorise training samples with a low noise level.Note that recent experiments in the investigation of Ca 2+ activity in astrocytes with precise spatial-temporal resolutions (Bindocci et al., 2017) revealed that global calcium event originated from the multiple foci on the most structural part of astrocyte periphery depends on the synchronous neuronal discharges.
Next, we studied the influence of synaptic connectivity architecture in the neural network, specifically the number, weight, and distribution of synaptic connections, on the correlation in the multi-item WM task (Fig. 9).The minimal number of synaptic connections, N out , required for existence of cued recall is 20 (Fig. 9a).A smaller number of connections is not enough to activate all the neurons from stimulus-specific population.Simultaneous increase of weights and number of connections induces the generation of large synaptic currents resulting in selfsustained overactivation of neuron-astrocyte network.Therefore there exists the optimal range of synaptic weight values to ensure high correlation.We found that for our model this range is η ∈ [0.005 − 0.05].Figure 9b illustrates the dependence of correlation on the number of output synaptic connections and their distribution.The smaller the parameter λ from Eq.( 4), the lower the probability of long-distance connections.The highest correlation was observed for local connections, λ < 7, due to the fact that short-range connections do not lead to blurring of the pattern retrieval boundaries in the neural network.Figure 9b also determined the optimal range for the number of synaptic connections: N out ∈ [25,55].
The key parameters, which determine the WM capacity in the proposed neuron-astrocyte network, are the duration of calcium signals in astrocyte and duration of the astrocyte-induced modulation of synaptic transmission.Duration of astrocytic calcium elevations is determined by the intrinsic mechanisms of the IP 3 -evoked Ca 2+ -induced Ca 2+ release from the astrocyte endoplasmic reticulum stores, which is described by the biophysical model (Li and Rinzel, 1994) used in this study.Fragmentary experimental data on duration of the gliotransmitter-induced modulation of synaptic transmission shows that the short-lived version of this modulation lasts from fractions of a second to few minutes, while long-term plasticity can last for tens of minutes (for review see (Pittà et al., 2016)).
To characterise memory capacity, we subjected it to longer trains of samples.Samples were applied to neuronal ensembles with average 35.2% overlapping in population.To check the memorisation, we presented cue items in the reverse order compared to learning mode (e.g.learning: 0, 1, 2, ..., 7, 8; test: 8, 7, ..., 2, 1, 0).The number of items with a correlation of recalled pattern higher than 90% indicated the capacity of the system.Figure 10 shows the capacity as a function of the sample number in the training sequence.For a chosen set of parameters, the capacity of WM ranges from five to six.Note, that such a limited capacity coincides with psychological studies indicating that human ability to keep information in readily accessible WM is limited, ranging between three and five items for most healthy people (Cowan, 2010).We investigated the influence of the duration of astrocyte-induced modulation of synaptic transmission, τ astro , on the capacity.The number of items stored in the system memory is maximum for parameter val-ues in the range: τ astro ∈ [60, 360] ms.For small values of τ astro excitation does not have time to spread to the stimulus-specific neural population; for long astrocytic modulation the different items in cued recalls interfered with each other.
In the case of non-overlapping neuronal ensembles, the WM capacity is unequivocally determined by the duration of the astrocytic calcium signal and can be obtained analytically.We estimated the Ca 2+ signal duration in τ Ca = 3.8 sec.We consider the case of training in K samples of a fixed order: 1, 2, ..., K.During the test, we look at a permutation p consisting of K patterns.For the permutation, we estimate the number of correctly recalled patterns, K p .Pattern is considered correctly recalled if no more than τ astro has passed since its presentation.The average capacity, C, for this case is defined as the average number of correctly recalled patterns over all possible permutations p.According to this description the average capacity can be calculated by the following equation: where τ 11 -train sample duration, τ 12 -duration between train samples, τ 21 -test sample duration, τ 22 -duration between test samples, τ shi f t -delay between train and test, τ Ca -calcium event duration, t i train -time of i-th train sample finishes, t j,i test .Figure 11 shows the capacity as a function of the sample number, K.After it reaches a maximum of 6.6 for 8 samples, the capacity begins to decrease monotonically, while the number of samples increases.

Discussion
We proposed a biologycally motivated spiking neuron network model accompanied by astrocytes that demonstrates the working memory formation.The model acts at multiple timescales: at a millisecond scale of firing neurons and the second scale of calcium dynamics in astrocyte.Neuronal network consists of randomly sparsely connected excitatory spiking neurons with nonplastic synapses.Astrocyte induced activity-dependent short-term synaptic plasticity results in local spatial synchronization in neuronal ensembles.WM realised by such astrocytic modulation is characterised by one-shot learning and is maintained during seconds.The astrocyte influence on the synaptic connections during the elevation of calcium concentration implements Hebbianlike synaptic plasticity differentiating between specific and non-specific activations.Note that the proposed model is crucially different from the attractor-based network memory models (Hopfield, 1982;DJ, 1995;Wang, 2001;Wimmer et al., 2014) and works similarly to WM models based on synaptic plasticity (Mongillo et al., 2008;Manohar et al., 2019;Lundqvist et al., 2011;Mi et al., 2017).In particular, in its functionality, the model is quite close to short-term associative (Hebbian) synaptic facilitation (Fiebig and Lansner, 2016;Sandberg et al., 2003).The concept of our WM model operation is schematically summarised in Fig. 12. Composed of two building blocks, e.g.fast-spiking neurons and slow astrocytes, the proposed memory architecture eventually demonstrated synergetic functionality in loading information and its readout by the neuronal block and storage implemented by the astrocytes.In contrast with solely neuronal circuit models where memory is encoded in synaptic connections and their plasticity, which inevitably leads to the problem of overlapping, our model splits functionality using astrocytes as a pool for stored patterns.Even with significant overlaps, they can be successfully retrieved due to coherent synaptic modulations by the astrocytes and synchronous neuron firing, which provide the selectivity.When the memory is maintained at the time scale of calcium elevation in astrocytes, the synapses are not specifically modulated and theoretically can be employed for other tasks.Note also that the memory is transient with no long-lasting changes in structure and parameters.Thus, new information can be successfully uploaded and stored without interference with traces of previously memorised information.Complete memory overwrite interval is estimated at several seconds and is also defined by the duration of astrocyte activations.
The proposed model clearly confirms the theoretical hypothesis that astrocytic modulation of synaptic transmission can be involved in the formation of functional cortical WM.We show that the NMDAR-mediated potentiation of excitatory synapses induced by the D-serine released from astrocyte in PFC could serve as a possible molecular mechanism for WM.Stages of the multiitem WM, which include loading, storage, and cued recall, manifest in brief oscillatory bursts, which are functionally similar to WM activity in nonhuman primate PFC (Lundqvist et al., 2016) rather than sustained neuron spiking.Interestingly, the activity of the neuron-astrocyte network corresponding to memorised pattern exhibits a sufficient degree of stability, which ensures memory retention despite the presence of significant overlaps in the stimulus-specific subnetworks.
Needless, to say that astrocyte-induced modulation of synaptic transmission proposed in this study as a mechanism for WM organisation does not exclude but rather complements other synaptic and neural plasticity mechanisms (fast Hebbian synaptic plasticity/short-term synaptic plasticity, facilitation, augmentation, dendritic voltage bistability, etc.) and may well act in parallel to them.
On the one hand, there has been much experimental evidence that astrocytes contribute to synaptic plasticity, coordination of neural network oscillatory activity, and memory function (Santello et al., 2019).It was shown recently that astrocytic impact is circuit-specific (Martin et al., 2015) and stimulus-specific (Mariotti et al., 2018).Improved Ca 2+ imaging approaches have identified a spatiotemporal diversity of astrocytic signals that may underlie the capacity of astrocytes to encode and process different patterns of activation (Bindocci et al., 2017;Stobart et al., 2018).Besides, the temporal scale of the astrocytic calcium dynamics and dynamics of the neuron-astrocyte bidirectional communication including the effects of astrocytic influence on synaptic plasticity fit very well in timing required in WM processes.
On the other hand, the ongoing intense debate about principles of WM organisation challenges the canonical theory of persistent delay activity in network attractors with recurrent excitation (Bouchacourt and Buschman, 2019) and offer alternative models incorporating different biophysical network mechanisms of WM (Lundqvist et al., 2018;Barak and Tsodyks, 2014).The principal reasons for such a debate is the reexamination of experimental data, which shows a large heterogeneity in the delay neuronal activity during WM tasks (Stokes et al., 2013).A nonclassical WM model includes the short-term synaptic plasticity (Mongillo et al., 2008;Hansel and Mato, 2013), the balance of inhibition and excitation (Boerlin et al., 2013), the NMDA currents affecting on the neuronal excitability (Durstewitz, 2009), and other parameters.These models, however, have a number of shortcomings: inability to describe encoding of novel asso-ciations in synaptic facilitation-based model; unclear mechanisms for achieving precise tuning of recurrent excitation and inhibition; the time constant of the NMDA receptor is appropriate to maintain memories for 1-5 s, but not longer.The investigation of the synaptic mechanisms underlying WM is an ongoing process.Therefore, incorporation of astrocytes as spatiotemporal integrator and modulator of synaptic transmission in neural network models may help advance the theoretical framework of WM encoding and maintenance mechanisms.
Future research direction in the framework of the proposed WM model in the neuron-astrocyte network will be focused on the interplay of excitation and inhibition that can stabilise WM [100]; the effects of synaptic plasticity (Boerlin et al., 2013); the effects of synaptic plasticity (Hansel and Mato, 2013;Mongillo et al., 2008) namely associative short-term potentiation (a fast-expressing form of Hebbian synaptic plasticity) that can provide an encoding of novel associations (Fiebig and Lansner, 2016); and on the structure of the cortical microcircuit reflecting columnar organisation of the neocortex.
This research was supported by the Ministry of Science and Higher Education of the Russian Federation (project No. 075-15-2020-808, No. 0729-2020-0061).SG thanks the RFBR .
In p u t s ig n a l N e u ro n a l n e tw o rk A s tr o c y ti c n e tw o rk

Figure 1 :
Figure1: Neuron-astrocyte network topology.Neuron-astrocyte network consists of two interacting layers: neural network layer and astrocytic layer.Neuronal network with dimension W × H (79 × 79) consists of synaptically coupled excitatory neurons modelled by the Izhikevich neuron.Neurons in network are connected randomly.The astrocytic network consists of diffusely connected astrocytes with dimension M × N (26 × 26).Blue lines show connections between elements in each layer.We consider the bidirectional interaction between the neuronal and astrocytic layers.Each astrocyte is interconnected with a neuronal ensemble of N a = 16 neurons with dimensions 4 × 4 with overlapping in one row (red lines).The input signal is fed to the neural network.

Figure 2 :
Figure 2: Model of neuron-astrocyte interaction.(a) Spike train and (b) concentration of neurotransmitter, G(t), of stimulus-specific neuron.(c,d) Same as in (a,b) but for an unspecific neuron.(e) Intracellular concentration of Ca 2+ and IP 3 in stimulus-specific astrocyte.Black bars at the top indicate periods when each of the stimuli (training stimulus -sample, nonmatching test items -nonmatch, test cue -match) was presented.In response to presynaptic spike train (a,c), the neurotransmitter, glutamate, G releases (b,d) into extracellular space and the concentration of IP 3 increases in astrocyte (e, blue line) inducing the elevation of intracellular Ca 2+ (e, red line).

Figure 3 :
Figure 3: Delayed matching to sample WM task.(a) One trial of the task simulated in the network.Spike raster of neuronal network showing sample-selective delay activity.Neurons belonging to stimulus-specific population are indicated by red color.Black bars indicate periods when each of the stimuli was presented.(b, c) The averaged firing rate of the stimulus-specific and unspecific neurons over time, respectively.(20 ms bins)

Figure 5 :
Figure 5: Neuron-astrocyte network simulation with four loaded memory items.(a,c,e) Spike train and (b,d,f) concentration of neurotransmitter, G(t), of three neurons belonging to different stimulus-specific populations.(a,b) Stimulusspecific neuron to sample 0. (c,d) Stimulus-specific neuron to sample 2. (e,f) Neuron unspecific to all samples.Black bars at the top indicate periods when each of the stimuli (training stimulus -sample, nonmatching test items -nonmatch, test cue -match) was presented.

Figure 6 :
Figure 6: Multi-item WM in neuron-astrocyte network.(a) Spike raster of neuronal network with 4 training patterns.Neurons are coloured according to their pattern selectivity.Pattern overlapping in neuronal populations is 35.2% on average for 4 patterns.Black bars indicate periods when each of the stimuli was presented.(b, c) The averaged firing rate of the stimulus-specific and unspecific neurons over time, respectively.(20 ms bins) (d) Correlation of filtered items.The different colours correspond to the correlations with different ideal samples.The dotted line shows the correlation of the testing item (for 20% noise level in test).

Figure 7 :
Figure 7: Snapshots of the training (a, b, c, f, g, h) and testing (d, e, i, j) neuron-astrocyte network in the multi-item working memory task.We used the following training set consisting of 4 samples: 0,1,2,3.(a,f) The example of first and last training samples, respectively.(b,g) The response of the neuronal network to samples.The values of the membrane potentials are shown.(c,h) The intracellular Ca 2+ concentrations in astrocytic layer.(d,i) The testing items with 20% of salt and pepper noise.(e,j) The cued recalls in the neuronal network.The firing rate averaged on the test time interval for each neuron is shown.

Figure 8 :
Figure 8: Correlation between recalled pattern (see section 3.5) and ideal item in a multi-item WM task performed by the neuron-astrocyte network.The correlation averaged over 4 patterns is shown.(a) Noise-resistance ability of the model.Dependence of correlation on noise level in training and test.The correlation difference between cued recall pattern and noisy input is shown.(b),(c) The influence of the neuron-astrocytic interaction structure.(b) Dependence of correlation on the number of spiking neurons required for the calcium elevation in the astrocyte, F act , and on the number of spiking neurons required for the emergence of astrocyte-induced enhancement synaptic transmission, F astro .(c) Dependence of correlation on noise level in training samples and on the parameter, F act .F astro = 0.5.For (b,c) noise level in cue is 20%.

Figure 9 :
Figure 9: Influence of synaptic connectivities architecture in the neural network on the correlation of recalled pattern in multi-item WM task performed by the neuron-astrocyte network.The correlation averaged over 4 patterns is shown.(a) Dependence of correlation on the number of output synaptic connections for each neuron, N out , and synaptic weight, η.(b) Dependence of correlation on the number of output synaptic connections for each neuron, N out , and synaptic connection distribution parameter, λ Eq.(4).

Figure 10 :
Figure 10: Capacity of the multi-item WM in neuron-astrocyte network.Capacity as a function of the sample number in training.The number of images with correlation of recalled pattern higher than 90% is shown.

Figure 11 :
Figure 11: Capacity of the multi-item WM in neuron-astrocyte network in the case of non-overlapping neuronal populations obtained analytically.Capacity as a function of the sample number in training sequence.

Figure 12 :
Figure 12: Concept of WM operation in spiking neuron network model accompanied by astrocytes.

Table 4 :
Stimulation protocol and recall testing parameters Specific input contains training samples in the form binary spatial patterns.The patterns represent different spatial distributions relative to background state with non-specific input only.
amplitudes randomly and uniformly distributed in the interval[−10, 10]µA.This input evokes background network state with low-rate, spontaneous spiking.