Spiking dynamics of single neuron is described using the Izhikevich model (Izhikevich, 2003). Due to its simplicity and computational efficient, this neuron model has been widely used to study population activities of neurons, e.g., synchronization and oscillation (Khoshkhou and Montakhab, 2018).

Mathematical expressions of the Izhikevich model are:

with the auxiliary after-spike resetting:

Here, V represents the transmembrane potential, while U denotes a membrane recovery variable that provides negative feedback to V. The indices (i, j) indicate the corresponding neuron. c is the resting potential, and a, b, d are dimensionless parameters. I_app refers to the applied currents to the respective neurons, which will be explained further below.

Synaptic currents that neurons receive is expressed as Gordleeva et al. (2021):

Here, N is the total number of presynaptic neurons, E_syn is the reversal potential for the synapse, V_pre denotes the membrane potential of the presynaptic neuron, and k_syn is slope of the synaptic activation function. The parameter g_syn describes the synaptic strength (Gordleeva et al., 2021):

Here (m, n) denotes the index of the corresponding astrocyte that modulates the synaptic currents of neuron (i, j). η represents the synaptic weight without astrocyte influence, while v_ca denotes the astrocyte-induced modulation of synaptic strength.

Specific values of these parameters are given in Table 1.

As mentioned above, astrocytes in our network primarily generate various calcium signals. Following the method used in the previous study (Gordleeva et al., 2021), a mean-field approach is employed to describe the emergence of Ca²⁺ signals, as shown in Equation 5 and Figure 1.

Here Ca²⁺ indicate the intracellular calcium concentration, h is the fraction of the activated IP₃ receptors on the ER surface. Detailed expressions of each flux are:

Here J_ER represents the Ca²⁺ flux from the ER to the cytoplasm through IP₃ receptors and ryanodine receptors (RyR). J_pump denotes the Ca²⁺ flux pumped back into the ER via the SERCA, while J_leak indicates the leakage flux from the ER to the cytoplasm. J_in and J_out describe the calcium exchange with extracellular space.

Activation dynamics of IP₃ is expressed as Gordleeva et al. (2021) and Wagner et al. (2004):

Here IP3* denotes the steady-state concentration of the IP₃ receptors, while τ_IP3 is the rate constant for IP₃ loss. λ is a control parameter used to regulate various elements in IP₃ dynamics. J_PLCδ describes the production of IP₃ by phospholipase Cδ (PLCδ), expressed as Gordleeva et al. (2021):

J_glu describes the production of IP₃ induced by glutamate in response to neuronal activities, which is modeled as Gordleeva et al. (2021):

Here A_glu represents the amplitude of glutamate contributing to the production of IP₃, while t_glu denotes the periods when the total level of glutamate from all synapses reaches a given threshold:

Here N_a represents the number of neurons connected to a single astrocyte, G_thr = 0.7 is the threshold, and [x] denotes the Iverson bracket. F_act = 0.5 indicates the fraction of synchronously spiking neurons within the total neuronal ensemble associated with the astrocyte. G refers to the amount of glutamate, which will be explained further below.

J_Kinase and J_Phosphatase are currents for IP₃ degradation, which are expressed as Wagner et al. (2004):

The currents df_Ca in Equation 5 and df_IP3 in Equation 7 represent the diffusion of Ca²⁺ ions and IP₃ molecules via gap junctions between astrocytes, expressed as follows (Gordleeva et al., 2021):

Here d_Ca and d_IP3 represent the diffusion rate of the Ca²⁺ and IP₃, respectively. Following a previous study (Gordleeva et al., 2021), we assume that each astrocyte is diffusively coupled only with its four nearest neighbors:

Specific values of these parameters are given in Table 2.

As described in Gordleeva et al. (2021), the network model for WM comprises three layers: the input layer, the neurons layer, and the astrocytes layer. The input layer consists of two types of stimuli: (1) image signals labeled {“0″, “1″, “2″, “3″, “4″, “5″, “6″, “7″, “8″, “9″} (Figure 2), which are responsible for memory training and testing, and (2) background noise signals, that generate low-rate spontaneous spikes.

The neurons layer contains 79 × 79 neurons, which receive stimuli from the input layer and project to the astrocytes layer. The third layer consists of 26 × 26 astrocytes, with each astrocyte projecting back to neurons in the second layer.

The architecture of synaptic connections among neurons is random; specifically, each neuron connects to N = 40 local postsynaptic target neurons via excitatory chemical synapses, with these neurons chosen randomly from the neurons layer. In contrast, the synaptic connections among astrocytes are deterministic, with each astrocyte connecting to its four nearest neighbors via gap junctions (Gordleeva et al., 2021).

Additionally, each astrocyte connects with N_a = 16 (i.e., 4 × 4) neurons through reciprocal excitatory chemical synapses. Thus, the synaptic strength in Equation 4 has a Ca²⁺ dependent component v_Ca, which can be expressed as:

Here vCa* denotes the strength of astrocyte-induced modulation of synaptic weight, Θ(x) is the Heaviside step-function, and [Ca2+]thr is the threshold.

In Equation 10, the amount of glutamate G is characterized as Gordleeva et al. (2021):

Here α_glu represents the glutamate clearance constant, while k_glu indicates the efficacy of the release. Similar to Equation 17, Θ(x) is the Heaviside step-function.

Specific values of relevant parameters are given in Table 3.

Our stimulation protocol follows the delayed match-to-sample task (DMS), which is commonly used in experimental studies of memory formation and recall. During the training phase, an image labeled with a specific digit is presented for 200 ms (500 ms-700 ms), followed by a 700 ms break (700 ms-1,400 ms). After this, two non-match images are displayed, each for 150 ms (1,400 ms-1,550 ms and 1,800 ms-1,950 ms), with a 250 ms break in between (1,550 ms-1,800 ms and 1,950 ms-2,200 ms). Finally, the matching stimulus appears for 150 ms (2,200 ms-2,350 ms) to test the network's memory recall capability. The total duration of the simulation experiment is 3,000 ms. Detailed stimulus information during training and testing phases can refer to Gordleeva et al. (2021). A demonstration of a sample stimulus current and the corresponding neuronal activities during the training and test phases is presented in Supplementary Figure 17.

To quantify the amount of memory that can be stored, we use a measure calculated based on the similarity between a recalled image and the sample image. The similarity ranges from 0 to 1, with larger values indicating better WM performance.

Here w = 1 ms represents the time step, thr denotes the spiking threshold of the neuron, P is the set of pixels belonging to the sample image, and W and H are the network dimensions. I is the indicator function, and T_P is the set of frames within the tracking range of pattern P. For more details (see Gordleeva et al., 2021).

The effect of Aβ on calcium dynamics is reflected in its influence on the currents J_in and J_ER. Following a previous study (Latulippe et al., 2018), J_in and J_ER with the addition of Aβ are expressed as:

Here A_β represents a fixed level of Aβ concentration.

Previous experimental results indicate that high concentrations of Ca²⁺ can promote the release probability of neurotransmitters, thereby strengthening synaptic connections between neurons (Neher and Sakaba, 2008). However, excessive levels of Ca²⁺ can disrupt synaptic connections by inducing irreversible excitotoxic injury (Vermma et al., 2022). In our model, Equation 21 is employed to prevent exaggerated Ca²⁺ concentrations.

Another experiment demonstrated that the accumulation of Aβ negatively affects Ca²⁺ levels (Toglia et al., 2018). Therefore, Equation 22 is introduced to simulate this variation.

Here A_β represents a fixed level of Aβ concentration. Based on the value of w_b, we have vCa*=wb * vCa*.

All simulations were performed using the Matlab software (Matlab R2021a), and the fourth-order Runge-Kutta algorithm was employed to calculate the values of different variables with a time integration step of 0.1 ms.

In the absence of Aβ (indicating no AD symptoms), the WM network performs effectively. Results shown in Figure 3 illustrate that the training stimulus “1” can be successfully recalled during the testing period after two non-match stimuli (“0” and “7”). Quantitative results in terms of similarity and peak frequency indicate that the similarity of neuronal responses to the sample stimulus and the match stimulus reaches ~0.9195. Additionally, the peak frequency of neurons responding to the match stimulus is comparable to that of neurons responding to the sample stimulus.

The dotted pink line in Figure 3B is determined by evaluating each stimulus from the set {“0″, “2″, “3″, “4″, “5″, “6″, “7″, “8″, “9″} as the match stimulus and calculating the maximum peak frequency among these stimuli (see Supplementary Figure 18). If the peak frequency of neurons under the match stimulus “1” exceeds this line, we consider neuronal activity under “1” to be distinguishable from the activities elicited by the other stimuli.

Results presented in Figures 4A, B illustrate the spiking sequence of a stimulus-specific neuron along with the amount of glutamate (G) it releases, demonstrating comparable spiking activities for both the sample and match stimuli. Figures 4C, D display the Ca²⁺ signal and IP₃ levels in an astrocyte that is that is synaptically connected to the stimulus-specific neuron. It is evident that once the glutamate release from the neuron reaches a certain threshold (“0.7”, as indicated in Equation 10), the IP₃ level rises rapidly. This continuous increase in IP₃ concentration gradually elevates the Ca²⁺ levels, ultimately enhancing neuronal responses. These findings highlight the critical role of Ca²⁺-mediated modulation in sustaining working memory, aligning with previous research (Mongillo et al., 2008).

Accumulation of Aβ has been linked to memory loss in AD in numerous studies (Hampel et al., 2021; Ma and Klann, 2012; Poling et al., 2008). In our research, we simulate Aβ accumulation by varying the A_β parameter and observe changes in memory performance. Figure 5 illustrates how similarity and peak frequency vary with respect to Aβ levels, showing that increased Aβ significantly impairs WM performance. Specifically, the peak frequency decreases from over 200 Hz (A_β = 0) to 93.87 Hz (A_β = 1.6), while similarity drops from 0.9195 (A_β = 0) to 0.541 (A_β = 1.6), indicating a substantial loss of memory recall fidelity.

In the following sections, we will separately analyze how memory performance is disrupted and the potential for restoration at A_β = 1.2 (mild AD) and A_β = 1.6 (severe AD).

The results in Section 3.2 indicate that excessive levels of Ca²⁺ induced by Aβ play a crucial role in disrupting the memory performance of WM network. In this section, we investigate whether downregulation of IP₃ can partially restore WM under pathological conditions by reducing intracellular calcium levels, as IP₃ positively influences intracellular calcium concentration (see Equations 5, 6).

In Figure 10, we analyze two cases under the condition of A_β = 1.2 (mild AD): (1) When λ = 0.1 (see Equation 7), similarity increases from 0.5832 to 0.7752, a rise of ~32.9%, and peak frequency increases from 98.98 Hz to 115.36 Hz, reflecting a 16.5% increase; (2) When λ = 0.05, similarity rises from 0.5832 to 0.7808, representing a 33.9% increase, while peak frequency increases from 98.98 Hz to 117.74 Hz, showing a 19.0% increase. These findings suggest that downregulation of IP₃ activation can significantly enhance WM performance under mild AD conditions.

However, in Figure 11, under the condition of A_β = 1.6 (severe AD), both λ = 0.1 and λ = 0.05 show minimal effects on similarity and peak frequency. This indicates that downregulating IP₃ activation does not improve WM performance in severe AD, underscoring the necessity for a dual-target modulation strategy to achieve full recovery in severe AD patients.

The effects of upregulating SERCA activation on WM performance are illustrated in Figure 12. In mild AD conditions, when γ = 1.25, similarity increases from 0.5832 to 0.6132, reflecting a ~5.1% increase, and peak frequency rises from 98.98 Hz to 101.58 Hz, a ~2.6% increase. However, when γ = 2.0, similarity increases from 0.5832 to 0.7710, representing a ~32.2% increase, while peak frequency rises from 98.98 Hz to 114.06 Hz, a ~ 15.2% increase (see Supplementary Figure 19). The corresponding variations in the Ca²⁺ signal for γ = 1.25 is shown in Figure 12B.

Moreover, simultaneously altering the activations of IP₃ and SERCA significantly enhances WM performance. Specifically, similarity rises from 0.5832 to 0.8105, a ~39.0% increase, and peak frequency increases from 98.98 Hz to 124.06 Hz, reflecting a ~25.3% increase. These findings suggest that concurrently downregulating IP₃ activation while upregulating SERCA activation plays a crucial role in restoring WM under mild AD conditions.

In Figures 12E–H, the model results for severe AD conditions indicate that upregulating SERCA activation has minimal impact on WM performance. Additionally, simultaneous changes to the activations of IP₃ and SERCA show only slight improvements in WM performance. These results imply that simply regulating intracellular calcium levels is insufficient to effectively restore WM performance in severe AD conditions.

To evaluate the general applicability of our network model, we conducted two additional experiments by changing the sample stimulus to two other digit stimuli: “2” and “7”.

First, our network model includes only excitatory neurons, overlooking the crucial role of inhibitory interneurons, which are essential for synchronizing network rhythms and preventing overexcitability. Future studies should incorporate inhibitory interneurons to enhance the realism of the memory model. Second, the astrocyte model is based on mean-field approximations, averaging calcium responses across compartments (Bazargani and Attwell, 2016). This approach simplifies the intricate spatial signaling seen in actual astrocytes, particularly in presynaptic domains. Integrating spatially distinct, cell-type-specific mechanisms would enhance the model's fidelity and biological relevance for future laboratory research.

In recent decades, the pathogenesis of AD has received considerable attention, leading to several prominent hypotheses: (1) Amyloid Hypothesis (Aβ): This hypothesis suggests that neuronal damage in AD patients arises from abnormalities in the processing of amyloid precursor protein (APP) and the subsequent accumulation of Aβ (Hampel et al., 2021). (2) Cholinergic Hypothesis: Supported by experimental observations, this hypothesis highlights a significant reduction in choline acetyltransferase, the enzyme responsible for synthesizing acetylcholine (ACh), in the amygdala, cortex, and hippocampus of postmortem brains from AD patients (Chen et al., 2022). (3) Glutamate Toxicity Hypothesis: This hypothesis is based on evidence of a marked reduction in the binding of 1-[3H] glutamate in the brains of individuals with AD (Maragos et al., 1987). (4) Tau Hypothesis: This hypothesis is grounded in the observation of aggregates of misfolded tau proteins present in the brains of AD patients (Frost et al., 2009).

While these four hypotheses aim to elucidate the pathophysiological basis of AD, it is important to recognize that the triggers associated with these hypotheses likely do not operate in isolation and may interact synergistically (Patow et al., 2023). In this study, we focused solely on AD conditions induced by Aβ accumulation. Future research should also examine AD conditions triggered by other significant factors.

Plasticity is a well-documented phenomenon at nearly all synapses in the brain, encompassing various types such as short-term plasticity, long-term plasticity (including long-term potentiation and long-term depression), Hebbian learning, and spike-timing-dependent plasticity (STDP). Moreover, synaptic plasticity is believed to play a critical role in the formation and storage of WM (Froudist-Walsh et al., 2018; Huang and Wei, 2021). In this study, however, the performance of the WM network is driven solely by intracellular calcium dynamics, without incorporating mechanisms of synaptic plasticity. Therefore, a potential extension of our network model could involve integrating plasticity mechanisms to explore changes in synaptic weights among cells.

Our network model operates at the cellular level, while the specific mechanisms of A_β function predominantly occur at the genetic level. For example, Hao and Friedman developed a mathematical model that outlines the Aβ aggregation process in detail, including production, clearance, and degradation (Hao and Friedman, 2016). Mustafa et al. created a two-compartment diffusion model to examine the effects of Aβ on the acetylcholine neurocycle, incorporating a differential equation to characterize changes in Aβ levels (Mustafa et al., 2021). Similarly, Chamberland et al. proposed a multiscale model consisting of 19 ordinary differential equations to analyze the progression of AD, introducing equations for both intracellular and extracellular Aβ (Chamberland et al., 2024). Models that focus on Aβ at the genetic level may offer a more accurate representation of its dynamic variations. Therefore, future modeling studies on Aβ-related AD generation and progression should incorporate more genetic-level components.

Memory impairments are common symptoms observed in the brains of AD patients. Alongside deficits in WM, other types of memory, such as short-term, episodic, and semantic memories (Jahn, 2013), are also affected during the progression of AD. Furthermore, various memory-related symptoms, including sleep disturbances (Pathmanathan et al., 2025) and accelerated long-term forgetting (Stamate et al., 2020), have been observed. Recent years have seen the development of computational models addressing these symptoms. For example, Horn et al. constructed a biologically motivated Hopfield model to investigate how the interplay between synaptic deletion and compensation influences memory deterioration patterns (Horn et al., 1993). Razi et al. built a network model of coupled cortical columns to explore the propagation modes between sleep and wakefulness (Razi et al., 2021). Additionally, Murre et al. introduced a mathematical model to analyze the dynamics of forgetting and amnesia, applying their findings to experimental observations from mice, rats and monkeys (Murre et al., 2013). Given this context, future computational studies could integrate these model explorations with our spiking network model to investigate potential therapeutic targets for these symptoms under AD conditions.