To study the mechanism of stress on working memory, we used a spiking neural network for modeling and analysis. Given the diverse and complex nature of actual stress behavioral experiments, we approached the analysis from the perspective of neural activity, using a simple short-term working memory recall paradigm, as shown in Figure 1A. This task paradigm begins with a 1-second fixation period, followed by a 400-ms long-term input stimulus (memory period), a 1,000-ms delay (intermediate delay period), and a 5-ms cue stimulus (cue prompt period). If a neuronal group can maintain discharge when the stimulus input to the network is applied, it indicates successful memory of the stimulus information. The input to the network involves simultaneously stimulating all excitatory neurons, with inputs following a Poisson distribution, as shown in Figure 1C. The intensity of the input stimulus is determined by frequency and changes according to a specified gradient.

The theoretical maximum memory capacity of our spiking neuron cortical network model is 4, corresponding to 4 distinct groups of neurons. As shown in Figure 1B, each group consists of excitatory neurons (red) and inhibitory neurons (blue), with the central blue region serving as the boundary between groups. In the figure, in addition to excitatory neurons, inhibitory neurons also exhibit continuous discharge in the groups that are capable of sustained activity. The fact that only two groups of neurons fired continuously suggests that the inhibitory neurons suppress the activity of the other groups of neurons. The figure shows that a maximum of two neuron groups can discharge continuously, indicating that the network's current memory capacity is 2. The network connection structure sets up four different neuron groups. When one of the groups discharges continuously, it means that the memory is successful, and the network memory capacity increases by one. The network's current memory capacity is 2 because, in its steady state, only two neuron groups can discharge continuously. The remaining two groups are suppressed by GABA currents, preventing them from firing.

Based on the existing working memory experimental paradigm and network structure, we model acute and chronic stress separately. As shown in Figure 1D, from a computational modeling perspective, the characteristics of acute and chronic stress are distinguished by the structure of dendritic spines and synaptic currents. In the neuronal model, the number of dendritic spines and the three synaptic currents AMPA, NMDA, and GABA are modeled, with changes in these characteristics based on real stress experiments used to simulate the effects of stress.

Animals exhibit synaptic organization principles that group neurons in specific ways (Perin et al., 2011; Pals et al., 2020). These components combine in various ways to form unique neural circuits (Bouchacourt and Buschman, 2019). Our network model represents “capacity” in the form of neuronal groups. All excitatory neurons receive input simultaneously, but due to differences in the organizational structure of synaptic connections, they exhibit varying capabilities for maintaining discharge. The specific neuronal group connection relationship is shown in Figure 1E. Each group has excitatory neurons and inhibitory neurons, which are connected to each other within the group and to neurons outside the group. Individual neuronal groups follow specific connection probabilities, are not fully interconnected, and exhibit variable connection strength. As shown in Figure 1G, the connection strength of neurons within a group is higher than that of neurons between groups. When all excitatory neurons receive external Poisson stimulation, this synaptic organization results in some neuronal groups maintaining discharge while others do not. The number of groups that can maintain discharge defines the network's current memory capacity (the maximum number of items that can be remembered). During the stimulation window, excitatory neurons within the same group facilitate each other's firing, while inhibitory neurons within the group actively suppress the firing of neurons in other groups. The network consists of four groups of selective neurons, but not all groups can maintain information simultaneously. Furthermore, they can be inhibited by inhibitory neurons from other groups. Thus, the ability to maintain information concurrently is influenced by the network structure and synaptic strength.

To model acute and chronic stress, the neuron model uses a two-compartment structure. The specific connection relationships between the two compartments are shown in Figure 1G. From a computational modeling perspective, both the dendritic and somatic compartments receive external information. In our two-compartment model, the dendrite receives excitatory and inhibitory synaptic currents from the external environment, which are transmitted through the axon and attenuated before reaching the soma. The soma integrates these currents from the dendrites and the external input to support its discharge activity.

To model the impact of acute and chronic stress on working memory tasks, we developed four distinct network models, as shown in Figure 1F. These models are generally divided into two categories: one based on a heterogeneous network using neural data from the 2/3 layer of the rodent prefrontal cortex (Hass et al., 2016), and the other a homogeneous network built with random parameters. The biological data for layer 2/3 of the mouse prefrontal cortex are comprehensive and closely related to the primary function (Shrestha et al., 2015) of working memory in mice. In addition, chronic stress has been shown to significantly affect this layer, leading to significant changes in dendritic structure. Based on the heterogeneous network structure, we created models for acute stress, chronic stress, and a control network. The heterogeneous network incorporates biological data, such as the connection probability between neurons, the distribution of pyramidal neurons and interneurons, neuronal electrophysiological parameters, and the average number of dendritic spines in the prefrontal cortex, as detailed in the Methods Section. Some parameters of the homogeneous network are derived from existing working memory models (Brunel and Wang, 2001; Rolls et al., 2013). Specifically, for acute stress, we used both the acute stress heterogeneous and homogeneous networks, while for chronic stress, we used the chronic stress heterogeneous network, the control heterogeneous network, and the homogeneous network. Since we only simulated the effect of acute stress on excitatory synaptic currents, the control network with a scaling factor set to 1 can be directly incorporated into the acute stress heterogeneous network without needing to re-establish the control network.

Previous research has demonstrated that activation of the cortical glutamate system can enhance working memory under acute stress (Kim, 2021; Popoli et al., 2012). However, how can the primary mechanisms by which glutamatergic activity affects working memory be quantified? What changes in network dynamics induced by acute stress promote working memory? As shown in Figure 2A, in vivo stressors act on pyramidal neurons via glucocorticoid receptors (GR), enhancing AMPA and NMDA synaptic currents in the prefrontal cortex (Kim, 2021; Yuen et al., 2011). Various types of acute stress, such as swimming for 20 min (Roche et al., 2003), confinement in a cubicle for 2 h (Mitra et al., 2005), placement on an elevated platform for 20 min (Xu et al., 1997), or a single cortical injection of ketamine (Yuen et al., 2009), have been shown to significantly enhance NMDAR-EPSC and AMPAR-EPSC activity in prefrontal cortical pyramidal neurons, thereby improving working memory capacity. Using forced swimming stress as an example (descriptions of the remaining three conditions are detailed in the Supplementary material), acute stress significantly enhanced NMDAR-EPSC (control: 197 ≤ 15 pA, n = 14; swim stress: 425 ≤ 20.5 pA, n = 15, P < 0.001, ANOVA) and AMPAR-EPSC (control: 58.6 ≤ 4.4 pA, n = 12; swim stress: 98.8 ≤ 3.7 pA, n = 12, P < 0.001, ANOVA) magnitudes in PFC pyramidal neurons (Yuen et al., 2009). To simulate this physiological process, we adjusted the ratio of AMPA to NMDA receptor-mediated excitatory postsynaptic currents in our model based on the synaptic current ratio observed in acute stress biological experiments (see Methods for details). Correspondingly, we increased the intensity of AMPA and NMDA channels and recorded network dynamics to observe working memory changes over different time scales (Geißler et al., 2023). Additionally, we conducted experiments on a homogeneous network to investigate the factors influencing biological data from rodent prefrontal layer 2/3 during stress-related working memory tasks. As shown in Figure 2E, we recorded the AMPA and NMDA current values of excitatory neurons. While these values do not perfectly match actual biological data, we observed that increasing the intensity of the NMDA channel fraction with a longer decay time constant resulted in an overall enhancement of current values.

The simulation results are shown in Figures 2B–D. Task performance is measured using two indicators: working memory capacity and time scale. By varying the channel score strength of AMPA and NMDA, we can observe the impact of these parameters on working memory through the result matrix. Specifically, the horizontal axis of each sub-plot represents increasing NMDA channel score strength, while the vertical axis represents increasing AMPA channel score strength. Figure 2B shows that both homogeneous and heterogeneous networks exhibit a significant increase in working memory capacity as the AMPA and NMDA channel score strength increases. Overall, the heterogeneous network demonstrates a higher working memory capacity than the homogeneous network, emphasizing the role of rodent biological data in enhancing working memory performance. The ability of neurons to maintain information is quantified by the time scale: the larger the time scale, the stronger the ability to maintain information. The calculation formula is provided in the Methods Section. Figures 2C, D together illustrate the time scale distribution of excitatory and inhibitory neurons during the memory and cue prompt stages as the AMPA and NMDA channel scores are varied. An increase in AMPA and NMDA channel score strength leads to a greater time scale, improving the neurons' ability to maintain information. The heterogeneous network shows a larger time scale overall compared to the homogeneous network, with a particularly notable increase in the time scale of inhibitory neurons. The enhanced discharge capacity of inhibitory neurons contributes to improved overall network performance in working memory tasks.

As shown in Figure 3A, there are various ways to induce chronic stress in animal models. Taking restraint stress as an example, rodents are subjected to daily 6-h restraint with a wire mesh for 21 consecutive days (Radley et al., 2004). This procedure leads to a significant reduction (16%) in the density of apical dendritic spines on pyramidal neurons in the medial prefrontal cortex (Radley et al., 2006).

Considering these existing experimental results, we simulated the the disruptive effects of chronic stress on GABAergic networks. In our model, we explored the impact of this disruption on working memory tasks in both homogeneous and heterogeneous networks by varying the ratio of GABA channels. We then examined strategies to compensate for this disruption. As shown in Figure 3B, we tested four approaches to modify the connection strength between neurons within a group. In this section, we begin by discussing how changes in GABAergic signaling, in the form of synaptic currents, can alter network dynamics. Biological experiments in rodents subjected to chronic stress also show changes in average GABA currents. For instance, in the control group, the average GABA BR-GIRK current was 152 ≤ 6 pA, while in the stress-adapted group, it was 141 ≤ 12 pA, and in the anhedonia group, it was 100 ≤ 5 pA. The stress-adapted group and the anhedonic group were distinguished by the sucrose consumption test. Under chronic stress, some animals exhibited resilience, with increased sucrose intake, and were classified as “stress-resilient.” In contrast, animals showing reduced sucrose intake were classified as the “anhedonic” subgroup. Figure 3C shows the schematic of average GABA current in inhibitory neurons as the GABA channel score strength changes. In the homogeneous network, increasing the GABA channel strength led to a monotonic decrease in the GABA current received by inhibitory neurons. Interestingly, in the heterogeneous network, the response followed an inverted U-shaped curve. At low GABA channel scores, increasing the channel strength enhanced the GABA current. However, as the inhibitory current reached a certain threshold, further increases in channel strength amplified inhibition of excitatory neurons, ultimately reducing the intensity of the GABA current in the inhibitory neurons.

To assess the performance of the network's working memory task under chronic stress, we used the pulse count autocorrelation decay time constant as a measure of neuronal time scale (Cavanagh et al., 2018; Wasmuht et al., 2018). A larger neuronal time scale indicates longer sustained neuronal discharge, which corresponds to better memory performance. We calculated self-decay curves under different intensity factors and plotted the time constants (time scales) of all curves separately as line graphs. The self-decay curves were calculated at two initial time points: one at the start of the memory phase and another at the start of the cue phase. If the time scale calculated from the cue phase is higher than that from the memory phase, it suggests that neurons are better able to maintain discharge during the cue phase, leading to improved memory ability. Chronic stress impaired the regulatory function of GABA currents, weakening their association with reduced working memory capacity. As shown in Figure 3D, as the GABA current increases, the firing capacity of all inhibitory neurons is significantly enhanced, and the excitatory neurons are inhibited when the intensity increases to a certain extent. In heterogeneous networks, the timescale of excitatory neurons under chronic stress conditions is generally smaller than that under normal conditions. As shown in Figure 3E, as the GABA current increased, the overall working memory capacity of the network decreased significantly. Not only was the time scale of the homogeneous network smaller than that of the heterogeneous network, but its working memory capacity was also lower. After adding rodent data to the network, both the time scale and working memory capacity increased. As shown in the table in Figure 3F, chronic stress reduced the correlation between the strength of the GABA channel score and the changes in the working memory capacity of the network.s shown in Figure 3D, the time scale of excitatory neurons in heterogeneous networks under chronic stress is smaller than under control conditions. Additionally, the time scale of homogeneous networks is smaller than that of heterogeneous networks, and their working memory capacity is also lower, as shown in Figure 3E. The inclusion of rodent data in the network increased both the time scale and working memory capacity. Chronic stress reduced the correlation between the fractional strength of the GABA channel and changes in network working memory capacity, as shown in the table in Figure 3F. Under control conditions, changes in GABA current strongly correlate with reductions in working memory capacity. However, under chronic stress, this correlation weakens, consistent with biological findings that chronic stress impairs the regulatory function of the GABAergic network.

There are various treatments for chronic stress, one of which is pharmacological intervention. Ketamine, administered at an antidepressant dose, has been shown to partially restore the baseline structure of dendritic spines in prefrontal cortex .Chronic stress induces a reduction in dendritic and synaptic density in rodent models, which is consistent with the changes in brain structure observed in patients with depression (such as MDD and PTSD) in clinical studies, especially in the hippocampus (HPC) and prefrontal cortex (PFC) regions (Koolschijn et al., 2009). This suggests that chronic stress not only changes the function of neural circuits, but also affects working memory function by changing neural structure.In network modeling, we adjusted the connection strength between neurons.Previous studies have demonstrated that stronger connections between excitatory neurons can slow down the dynamic changes of neural networks (Chaudhuri et al., 2015; Wang, 2008). Therefore, we proportionally altered the connection strength between excitatory neurons (E-E) within the same group in the network (Lam et al., 2022). We also increased the NMDA current intensity to simulate the process of ketamine treatment of chronic stress. The experimental conclusions were similar to those of increasing the intensity of excitatory neurons. The results are shown in Supplementary Figure 6. To explore these effects, we applied two input stimuli with different durations but the same intensity, representing the memory phase and the cue prompt phase, to both homogeneous and heterogeneous networks. We then analyzed the resulting changes in network time scales and working memory capacity. Furthermore, to investigate the impact of various network connections on task performance, we conducted experiments altering I-E, I-I, and E-I connections, as illustrated in Figure 3B. The detailed results of these simulations are provided in the Supplementary material.

Figure 4A shows that in the homogeneous network model, excitatory neurons with enhanced EE connectivity exhibit a rapid increase in time scale. In heterogeneous networks, the time scale increase of excitatory neurons is more obvious. In addition, as the E-E connection strength increases, the time scale of inhibitory neurons in the heterogeneous network model changes in an inverted U shape: initially, the time scale of inhibitory neurons increases due to receiving more inputs, and the E-E connection strength increases; however, when the E-E connection strength is further enhanced, inhibitory neurons peak as a result of receiving stimulation time scales from more excitatory neurons, and inhibitory neurons with enhanced activity begin to inhibit other neurons and thus overall inhibitory neurons activity begins to decrease. This shows that rodent-specific biological information can help the network become more flexible and adaptable in the face of external input. From a time scale perspective, the reduction in synaptic density of excitatory neurons has no significant impact on network dynamics. This suggests that strong connections between excitatory neurons may compensate for the decrease in synaptic density in chronic stress models. But the results of the other three connection-changing experiments were not as significant as directly changing the connection strength between excitatory neurons, which had a more obvious positive impact on working memory performance. The specific results are shown in Supplementary Figures S3–S5. Compared with the homogeneous network, the overall working memory capacity of the heterogeneous network is higher, as shown in Figure 4B. We observed a strong correlation between increasing working memory capacity and increasing E-E connection strength in all three networks, as shown in Figure 4C. But we did not find this relationship in the I-E, II, and E-I connections. In the heterogeneous network, the impact of changes in these connections on working memory capacity was not significant. In Figures 4B, C, the working memory capacity of homogeneous networks was overall lower than that of heterogeneous networks in all cases, and the correlation between excitatory connections and working memory capacity was the lowest. This finding is consistent with Figures 3E, F. This fully demonstrates the importance of biological data for improving network working memory capacity.

Our network is a feedforward spike neuron network, in which neurons update spikes through postsynaptic currents. Specifically, we studied how the currents of other neurons affected by short-term synaptic plasticity are transmitted to the target neurons in the neuron model. This study aims to simulate different stress states by updating the network state through synaptic dynamics. The network consists of a single layer of spike neurons, and external stimuli are input to excitatory neurons through Poisson coding. Excitatory neurons stimulate the intra-group neurons and inter-group neurons connected to them. We use a short-term working memory task paradigm to record the spike discharge state of neurons updated at each time step. When the effect of changes in synaptic connections on neuronal discharge can be directly counted, the ability of neurons to maintain information changes dynamically, and the entire neuronal group as a whole shows a large amount of discharge or no discharge. The discharge of different neuronal groups constructs the concept of discrete values of working memory capacity, which is used to measure the task effect under current conditions.

The connection probability and number ratio of excitatory neurons and inhibitory neurons in the network were determined using published mouse in vivo anatomical data (Hass et al., 2016), and the electrophysiological parameters of pyramidal neurons and intermediate-type neurons in the second and third layers of the rodent prefrontal cortex were selected. The effects of acute stress on mouse pyramidal neurons were achieved by exposing animals to different stressors, such as: forced swimming, acute restraint stress, elevated platform stress, and corticosterone injection (Yuen et al., 2011). Changes in neuronal AMPA and NMDA currents were caused and incorporated into the modeling as a scaling factor for synaptic currents. The effects of chronic stress on neuronal structure in mice were investigated by daily restraint stress or exposure to repeated swimming stress (Strekalova and Steinbusch, 2010), which resulted in the reorganization of the apical dendritic structure of pyramidal neurons in the rat prefrontal area (Radley et al., 2006; Cook and Wellman, 2004), reducing the density and total length of apical dendritic spines (Radley et al., 2006). In order to simulate the changes in dendrite morphology of pyramidal neurons, the neuron model uses a multi-compartment model and introduces the newly published excitatory nerves in the mouse L2 layer ACC brain area (anterior cingulate cortex) measured using 3D electron microscopy (Loomba et al., 2022). The number of spine neurons and inhibitory neurons, among which inhibitory neurons are divided into MP (multipolar) and BP (bipolar) according to their morphology (Loomba et al., 2022). According to existing literature, synapses on mouse dendritic spines are mainly excitatory (Braitenberg and Schüz, 2013), while synapses on dendritic shafts are mainly inhibitory (Karimi et al., 2020), therefore it is necessary to distinguish the number of excitatory synapses and inhibitory synapses on different types of neurons according to the proportion of inhibitory inputs to the corresponding neurons (Loomba et al., 2022). The specific value table is shown in the Supplementary material.

Each population of identically preferred neurons consists of excitatory and inhibitory populations. Our model is based on Integrate-and-Fire spiking neuron (IF) (Abbott, 1999; Hertäg et al., 2012; Zeng et al., 2023):

Vdenotes the subthreshold membrane potential of neurons within the ith selective population, with membrane capacitance C_m, resting potential V_L, and synaptic current I_syn.

Synaptic currents include glutamatergic excitatory components mediated by AMPA and NMDA receptors, inhibitory components mediated by GABA receptors, dendritic currents, and background currents. External stimuli contribute currents through AMPA receptors. β_NMDA and β_AMPA are factors used to measure the multiples of current growth in biological experiments during acute stress conditions (actual acute stress current size/control group current size). The acute stress model in the biological model corresponds to the proportional factors of AMPA and NMDA current changes under three stress stimuli. Among them, external stimuli are input through I^{AMPA, ext} current. The total current is shown in the figure below:

The number of dendritic spines affects the reception of synaptic current (Froudist-Walsh et al., 2021). We established the relationship between the number of dendritic spines and synaptic current as follows:

Among them, EXN represents the average number of dendritic spines of excitatory neurons, BP represents the average number of dendritic spines of bipolar neurons, and MP represents the average number of dendritic spines of multipolar neurons. EXN_IN represents the estimated value of the inhibitory input fraction [i/(i+e)] on the dendrites of mouse excitatory neurons, BP_IN and MP_IN represent the inhibitory input fraction values of bipolar interneurons and multipolar interneurons, respectively. When simulating the scenario of chronic stress, the percentage of dendritic spines of excitatory neurons decreased, fac was 0.16 in chronic stress and 0 in other cases. The specific values above are from the experimental literature (Loomba et al., 2022). The subscript e represents excitatory neurons, and i represents the value of inhibitory neurons.

The NMDA current is calculated as shown below:

μ^{NMDA, rec} corresponds to the proportion of the number of NMDA receptor synapses of actual excitatory pyramidal neurons in the biological model, taking the average number of excitatory synapses. g^{NMDA, rec} is the synaptic conductance of the receptor NMDA in excitatory neurons, V_L is the reversal potential, s^{NMDA, rec}is the fraction of receptor NMDA open channels, and w_j is the synaptic weight. The factor u_j regulates excitatory synapses. γand β are constants.N is the number of neurons.

The AMPA current is calculated as shown below:

where gEAMPA,ext is the synaptic conductance of the receptor AMPA in excitatory neurons, and s^{AMPA, ext}is the fraction of receptor AMPA open channels that also receive external input. In most simulations, external EPSCs were mediated exclusively by AMPA receptors. In a few simulations, we introduced NMDA receptors in external inputs (Brunel and Wang, 2001).

μ^{AMPA, rec} corresponds to the proportion of the number of AMPA receptor synapses of actual excitatory pyramidal neurons in the biological model, taking the average number of excitatory synapses. g^{AMPA, rec} is the synaptic conductance of the receptor AMPA in excitatory neurons, and s^{AMPA, rec} is the fraction of receptor AMPA open channels. N is the number of neurons.

The GABA current is calculated as shown below:

where g^GABA is the synaptic conductance of the receptor GABA in inhibitory neurons, and s^GABA is the fraction of receptor GABA open channels. N is the number of neurons.

The dendritic currents current is calculated as shown below:

I^dend,exc is the excitatory current flowing into the dendrites, I^dend,inhis the inhibitory current flowing into the dendrites, and I^soma,dend is the current from the dendrites to the soma. c1 to c6 control the gain, translation, inversion point, and shape of the nonlinear function (Marlin and Carter, 2014).

After introducing a short-term synaptic plasticity mechanism for excitatory synapses, as follows (Mi et al., 2017; Froudist-Walsh et al., 2021; Mongillo et al., 2008):

Among them, U=0.15, τ^u = 1, 500ms, τ^x = 2ms, γ_AMPA, γ_NMDA are two factors used to regulate the magnitude of the two currents under acute stress. The coefficient factor γ_I is used to regulate the GABA current and is used to change the proportion of the GABA current under chronic stress.

To calculate the neuron time scale we calculated the decay time constant of the spike-count autocorrelation function for each unit during the fixation period (Murray et al., 2014). We calculated the neuron firing sequence for the first 5 seconds of the network, using the 50 ms Convert the neuronal firing sequence into a spike count matrix for a time window of unit length. Next, we calculated the Pearson correlation coefficient between the two time windows (corresponding columns of the spike count matrix) lagged Δ. Fit an exponential decay function (ρ̄) to the Pearson correlation coefficient at each lag value using Levenberg-Marquardt (Moré, 2006; Kim, 2021) nonlinear least squares method:

where a is the fitted amplitude and b is the offset. τ is the time scale of the neuron, measuring the decay rate of the autocorrelation coefficient.