Complex non‐Markovian dynamics and the dual role of astrocytes in Alzheimer’s disease development and propagation

TL;DR

This study explores how astrocytes and memory effects influence Alzheimer's disease progression using a mathematical model.

Contribution

The paper introduces a fractional order model capturing astrocytes' dual role and memory effects in Alzheimer's.

Findings

Astrocytes can both clear and contribute to Aβ accumulation in Alzheimer's.

Memory effects slow the spread of toxic proteins in the brain.

Brain regions show varied Aβ and tau accumulation due to mixed tauopathy.

Abstract

Alzheimer’s disease (AD) is a common neurodegenerative disorder nowadays. Amyloid‐beta (Aβ) and tau proteins are among the main contributors to the AD progression. In AD, Aβ proteins clump together to form plaques and disrupt cell functions. On the other hand, the abnormal chemical change in the brain helps to build sticky tau tangles that block the neuron’s transport system. Astrocytes generally maintain a healthy balance in the brain by clearing the Aβ plaques (toxic Aβ). However, overactivated astrocytes release chemokines and cytokines in the presence of Aβ and react to pro‐inflammatory cytokines, further increasing the production of Aβ. In this study, we construct a mathematical model that can capture astrocytes’ dual behavior. Furthermore, we reveal that the disease progression depends on the current time instance and the disease’s earlier status, called the “memory effect,” making non‐Markovian processes an appropriate approach. We consider a fractional order network mathematical model to capture the influence of such memory effects on AD progression. We have integrated brain connectome data into the model and studied the memory effect, the dual role of astrocytes, and the brain’s neuronal damage. Based on the pathology, primary, secondary, and mixed tauopathies parameters are considered in the model. Due to the mixed tauopathy, different brain nodes or regions in the brain connectome accumulate different toxic concentrations of Aβ and tau proteins. Finally, we explain how the memory effect can slow down the propagation of such toxic proteins in the brain, decreasing the rate of neuronal damage.

Figures

The sensitivity analysis of different parameters on the nonfractional model of Equations (3) and (4).

Weighted adjacency matrix for the brain connectome data: limbic (L), frontal (F), basal ganglia (B), parietal (P), temporal (T) and occipital (O). In our integrated brain connectome data, each region contains one or more brain IDs, and they are listed in Table 2.

Solutions for the nonfractional model of Equation (6) in the brain connectome. (A) Solutions in each region and (B) the spatial average solution. The fixed parameter values are given in Table 1.

Spatial average solutions of toxic Aβ (dotted) and toxic τP (solid) for the nonfractional model of Equation (6) for different values of c 1 and μ: (red) μ = 0.1, (green) μ = 0.2 and (blue) μ = 0.3. (A) c 1 = 0.3 and (B) c 1 = 0.1.

Spatial average solutions of toxic amyloid‐beta (dotted), toxic tau protein (solid) and damage (dashed) for the fractional model of Equations (6) with (7) for different values of α with the other fixed parameter values of Table 1 over the brain connectome. In the plot, α = 1, α = 0.9, and α = 0.8 are represented by the red, green and blue curves, respectively.

Brain ID‐wise average toxic tau protein propagation (A) and the corresponding neural damage (B) for the nonfractional model. The y‐axis represents the brain‐ID‐wise average toxic tau protein concentrations.

Three‐dimensional views of the positions of the nodes for the integrated brain connectome data. (A) Axial view and (B) sagittal view. Different colors are used to indicate different brain regions.

Brain region‐wise average toxic tau protein (v˜) propagation for the nonfractional (A) and fractional with α = 0.8 (B) models. The y‐axis represents the average toxic tau protein concentrations.

Node‐wise damage propagation (q) for the traditional nonfractional (upper panel) and the fractional with α = 0.8 (lower panel) models in brain connectome. The dark red represents the high damage, and the light yellow represents the low damage.

Tables

TABLE 1. Synthetic parameter values [45].

Parameter	Value	Parameter	Value	Parameter	Value	Parameter	Value
a0	1.035	a1	1.38	a2	1.38	a˜1	0.828
b0	0.69	b1	1.38	b2	1.035	b˜1	0.552
b3	4.14	c0	1.0	c1	0.1	μ	0.1
ρ1	1.38	ρ2	0.138	ρ3	1.38	ρ4	0.014
k1	0.0001	k2	0.01	k3	0.1	k4	0.001

TABLE 2. Brain IDs associated with brain regions [71].

Brain region	Brain ID
Limbic	Rostralanteriorcingulate, posteriorcingulate, caudalanteriorcingulate, parahippocampal, isthmuscingulate, and entorhinal
Frontal	Frontalpole, lateralorbitofrontal, parsorbitalis, medialorbitofrontal, precentral, parstriangularis, parsopercularis, rostralmiddlefrontal, caudalmiddlefrontal, and superiorfrontal
Parietal	Postcentral, supramarginal, superiorparietal, inferiorparietal, precuneus, and paracentral
Basal ganglia	Left‐thalamus‐proper, left‐putamen, left‐caudate, left‐accumbens‐area, left‐pallidum, left‐amygdala, right‐thalamus‐proper, right‐putamen, RightCaudate, right‐accumbens‐area, right‐pallidum, and right‐amygdala
Brain stem	Brain‐stem
Occipital	Cuneus, pericalcarine, lateraloccipital, and lingual
Temporal	Middletemporal, left‐hippocampus, right‐hippocampus, temporalpole, inferiortemporal, bankssts, superiortemporal, transversetemporal, insula, and fusiform

Full text

INTRODUCTION

1

Alzheimer’s disease (AD) is a neurological disorder that worsens with age and is incurable. It affects thinking, memory, and behavior. These cognitive declines may be so severe that they interfere with daily tasks. Dr. Alois Alzheimer first observed this disease in 1906 and described it as “a peculiar disease” [1]. He studied the brain of a lady who had passed away from an uncommon mental condition that included memory loss and linguistic issues. He found many abnormal clumps (amyloid‐beta plaques) and fiber bundle tangles (tau tangles), which are now considered as one of the main contributors to AD progression [1, 2, 3]. They block communication between nerve cells, disrupt many processes, and cause memory loss, difficulty in speaking and other cognitive declines.

AD is not a natural part of the aging process although the chance of developing it grows with age [4, 5]. The majority of AD’s patients are 65 years or older, although the disease can develop before that age, a condition known as early‐onset illnesses AD [6]. Early in the course of the disease, people with Alzheimer’s disease have little memory loss but as they become older, they progressively lose their ability to maintain focus on a conversation or recall familiar faces. AD patients can survive up to 20 years following diagnosis although the average is eight to ten years [7, 8]. AD does not yet have a treatment, except a few drugs, such as aducanumab, which can assist in slowing cognitive loss in the early stages [9, 10, 11]. Many researchers have been working worldwide to understand the disease in a better way and prevent it from spreading.

Toxic amyloid‐beta (Aβ) accumulation in the extracellular space is commonly regarded as one of the key initiators of the early start of AD [12, 13, 14, 15, 16]. This accumulation may happen due to its overproduction or the lower clearance rate [17]. Aβ comprises 39–43 amino acids with different biophysical states, and soluble Aβ_40_ and insoluble Aβ_42_ are the two major isoforms observed in the brain. In a healthy brain, over 90% concentration of Aβ is detected in the form of Aβ_40_, whereas less than 5% can be found in the form of Aβ_42_ [18, 19, 20]. Unfortunately, an AD‐affected brain has no ability to maintain this state of balance and higher levels of Aβ_42_ peptide lead plaques to develop, disrupting cell function. In addition, tau protein (τP) has a crucial function in AD [21, 22]. The normal τP creates a microtubule that helps transmit nutrients and other substances from one area of a nerve cell to another [23]. Abnormal chemical changes in the brain cause tau protein to separate from microtubules and attach to other τP [24, 25]. This causes the tau protein to form neurofibrillary tangles (misfolded and abnormally shaped) inside neurons and block the neuron’s transport system.

Researchers have focused on identifying toxic amyloid‐beta and tau protein concentrations at the early onset of AD. At present, the accumulation of these proteins has not been completely quantified using blood tests and cannot be observed on CT or MRI images. The FDA‐approved amyloid PET scan tracer can identify the existence of Alzheimer’s disease, but it cannot adequately monitor disease progression; therefore, it is only used in clinical studies [26, 27, 28]. In contrast, F‐18 flortaucipir is the first FDA‐approved tau PET scan tracer that aids with the progression of AD neurodegeneration [29, 30]. Along with these two proteins, many other factors influence AD progression. Substantial efforts have been made to identify the disease state based on different factors [31, 32].

A specific form of glial cell is found in the central nervous system (CNS), which serves as an immunodefense to the CNS. They control blood flow, transport mitochondria to neurons, and contribute to neuronal metabolism [33, 34, 35, 36]. Active astrocytes generally clean waste from the brain and protect neurons against illness [37, 38, 39, 40]. But, in the AD‐affected brain, they lose the ability to maintain a healthy balance and support AD progression [41]. In the early stages of AD, a sufficient amount of toxic amyloid‐beta mainly disrupts this healthy balance. In this circumstance, astrocytes cannot maintain the brain’s ionic equilibrium, particularly intracellular Ca^2+^ concentrations. As a result, NADPH oxidase (NOX) is activated, and neuronal death occurs due to oxidative stress [42, 43]. Many other detrimental repercussions of astrocyte overactivation occur, including apolipoprotein E (ApoE) and excessive glutamate production. ApoE4 is a neurotoxic isoform of ApoE that contributes to toxic Aβ deposition during the early stages of Alzheimer’s disease [6, 44].

One of the main aspects of our current work is to analyze astrocytes’ dual role before and after AD. We developed a mathematical model that incorporates the involvement of astrocytes in AD along with the Aβ and τP interactions. In reaction kinetics, each of these proteins (Aβ and τP) follows a heterodimer model for interactions between proteins, with a coupling parameter between them [45]. We modify the exponential growth by logistic growth in the growth term for both the healthy proteins’ equation [46, 47]. We consider a logistic expansion in the astrocyte equation and presume they eliminate toxic amyloid‐beta [48, 49]. Furthermore, these toxic proteins damage the neurons in the brain connectome. The amount of neuronal damage is studied here by coupling toxic Aβ and τP, which also gives the disease status in the brain [45].

A key challenge in modeling complex biological systems is extracting meaningful insights from available data. To address this, researchers have explored various modeling approaches that enhance interpretability and predictive power. Fractional calculus offers significant advantages over traditional integer‐order models, as it naturally captures non‐Markovian dynamics [50, 51, 52]. Due to its global correlation, it can reflect the historical process of the systematic function and act as nonlocal interactions [53]. Regarding the data fitting, it has been observed that the fractional model has one more degree of freedom over the traditional integer‐order model [54]. Considering that a reaction–diffusion process could rely not only on the previous time instance’s concentrations but also on each of the past stages of concentrations with specified weights, which is further discussed in this study [52, 55, 56, 57, 58, 59].

Fractional‐order derivative models widen classical calculus by expanding differentiation to noninteger orders, including memory effects and long‐range interactions common in complex systems. They are specified using integral formulas, such as the Riemann–Liouville and Caputo derivatives, making them useful for simulating anomalous diffusion, viscoelastic materials, and biological processes [50, 58, 60]. Fractional operators, such as the fractional Laplacian, use integrals over entire domains to capture nonlocal behavior when extended to subsets of Euclidean space, often requiring specialized boundary conditions. This framework enables the formulation of fractional partial differential equations on restricted domains, which are commonly applied in fields, such as mathematical biology, image processing, and control theory, where classical integer‐order models may fall short in representing complex spatial‐temporal dynamics.

A time‐fractional reaction–diffusion equation is often used to explore the memory effects in AD processes, as it can capture the influence of past states on the present dynamics. Investigating such memory effects in AD patients is crucial, as the disease progression and recovery rates can vary significantly between individuals. For instance, the recovery rate of a patient who has been affected by the disease for 20 years may differ from that of a patient who has been diagnosed for only 10 years, highlighting the need for models that account for long‐term memory and individual variability in disease dynamics. There is a growing number of works on fractional models, which are applied to mathematical biology and other areas, such as fractional reaction–diffusion models in pattern formation and the dynamics of chemical kinetics in a heterogeneous setting [61], a typical measure of reaction time in thermally activated barrier‐crossing processes [62], and long‐range movement of certain organisms in the presence of a chemoattractant [63], gene expression [50], cell motions [51], long‐range memory [64], etc. With several additional practical uses, the Caputo derivative is among the best fractional operators for use in this kind of modeling [57, 65, 66]. In this work, we first construct a time‐fractional partial differential equation (PDE) model to describe the AD progression. Then, we develop the network model that aligns with the PDE model to integrate the data on brain connections and examine the damage dynamics associated with the influence on memory. Different tauopathies have been studied for the network model to incorporate different scenarios depending on the toxic Aβ and toxic τP. Furthermore, we have compared the disease progression in the absence and presence of memories on nodes and brain connectome regions.

The rest of this article is organized as follows. Section 2 presents comprehensive numerical simulations to investigate the dual role of astrocytes and the memory effect in AD progression. Section 3 discusses the key findings and explores future research directions. In Section 4, we develop temporal models for AD progression under two scenarios: with and without memory effects. Additionally, we analyze the equilibria and their stability properties. Furthermore, the temporal model is extended within a subset of Euclidean space and integrated into a network framework to incorporate brain connectome data.

RESULTS

2

This section presents the numerical results for the nonfractional and fractional models applied to the brain connectome. Before conducting numerical simulations, we introduce the synthetic parameter values used in the models, listed in Table 1, as estimated by Thompson et al. [45]. Additionally, we have performed a sensitivity analysis, which is essential for assessing how variations in parameter values affect the model’s behavior. This analysis helps identify key parameters that drive system dynamics, emphasizing those with a significant impact on the outputs, which may require precise estimation or experimental validation. Figure 1 illustrates the sensitivity analysis results for the nonfractional model of Equations (3) and (4) in MATERIALS AND METHODS. Pearson correlation coefficients were calculated for each parameter by generating one thousand uniform random samples with a 10% deviation from the values listed in Table 1 and evaluating the model’s solution at t = 200. In this case, the initial condition is chosen as u0=0.75,u˜0=0.0075,v0=0.5,v˜0=0.005, and q0=0. The bar plot reveals that parameter a 0 exhibits a strong positive correlation with u,v, and q, while it shows a strong negative correlation with v and a weaker correlation with both u and w variables. Sensitivity patterns for the other parameters are also illustrated in the figure. Overall, the parameters on the left side of the x‐axis label in Figure 1 exhibit stronger correlations with q. Therefore, they may serve as potential targets for AD‐modifying therapies.

The sensitivity analysis of different parameters on the nonfractional model of Equations (3) and (4).

We use the brain connectome data, which is available freely at BrainGraph.org—the network of the brain. These data give information on a network with nodes and edges in various brain areas, allowing us to investigate the brain’s spatiotemporal behavior. In this brain graph data, each node corresponds to a tiny area (1−1.5 cm^2^) of the gray matter, called the region of interest (ROI). An edge may be connected to two nodes if a diffusion‐MRI‐based procedure discovers fibers of axons going between those two nodes in the brain’s white matter [67, 68]. We have integrated the brain connectome data into our computational environment (MATLAB) and extracted the corresponding Laplacian for the real data. The network data consists of N = 1015 nodes and E = 16,280 edges. The number of fibers in the integrated data varies between 1 and 4966.5 with an average of 39.33 fibers per edge. Most of these fibers are located between the superior parietal and precuneus regions [69]. Additionally, the average fiber length (l _ ij _) ranges from 10.270 mm to 83.003 mm, with an average of 30.089 mm. We use the Laplacian to derive the numerical solution for the network model to handle various scenarios. The fourth‐order Runge–Kutta method is applied to integrate the resulting system with a time step of dt = 0.01, and the results remain consistent for smaller time steps. Furthermore, the predictor corrector method is employed to solve the fractional model [70]. We have computed the numerical results based on our code implemented in C‐language. We have used SHARCNET supercomputers to run multiple jobs simultaneously, which helped us efficiently analyze different aspects of the model.

Figure 2 depicts the weighted adjacency matrix for the considered network data. In the plot, blue to red colors represent the minimum to maximum strength of the connection between the nodes. The integrated brain connectome data contains one or more brain IDs in each region, listed in Table 2. The brain stem region is also in the diagram between the frontal and basal ganglia, but it is not visible because it contains only one node. This figure depicts the relationship between the brain connectome’s nodes and regions. This matrix determines the spatiotemporal behavior of the brain connectome network model. According to the parameter values, both models (nonfractional and fractional models) share the feasible homogeneous steady states. In the homogeneous steady‐state, the concentration of toxic amyloid‐beta could be governed by the concentration of toxic tau protein. In this case, it is called secondary tauopathy; otherwise, it is a primary tauopathy. We will discuss both cases in the coming subsection.

Weighted adjacency matrix for the brain connectome data: limbic (L), frontal (F), basal ganglia (B), parietal (P), temporal (T) and occipital (O). In our integrated brain connectome data, each region contains one or more brain IDs, and they are listed in Table 2.

Before moving to the numerical simulations, we first mention the initial conditions for each variable in the network model. In the brain connectome, the initial seeding sites for the toxic amyloid‐beta are the temporobasal and frontomedial regions, and the toxic tau proteins are the transentorhinal and locus coeruleus associated regions [45, 47]. For the seeding locations, we add tiny toxic concentrations of 0.25% and 0.38% in toxic tau protein (v˜) and toxic amyloid‐beta (u˜), respectively. Due to these small perturbations, the toxic concentrations propagate all over the brain connectome and spread AD. On the other hand, we consider healthy concentrations for both amyloid‐beta (u) and tau proteins (v) and a small concentration for the astrocytes (w). Some other perturbations of these initial concentrations can change the initial propagation profiles of the concentrations, but the final results (long‐term behaviors) are the same. These concentrations are uniform on each node in the brain network. We set the initial condition q=0 for the damage equation to each node. As toxic loads propagate over the brain connectome, they damage the neurons in the brain.

Primary and secondary tauopathies

2.1

As in our previous study, we have shown that the evolution profiles of both toxic loads remain the same for primary and secondary tauopathies in the absence of memory (α=1) and astrocytes [47]. We also observe that the profiles remain consistent for both tauopathies in the presence of astrocytes. Therefore, without any loss of generality, we present the results using the parameter values associated with the secondary tauopathy. Table 1 provides a synthetic parameter set specific to the secondary tauopathy. For this parameter set, the nontrivial equilibrium point E ∗ = (0.596, 0.154, 0.33, 0.14, 0.1) is locally asymptotically stable, and a numerical solution of the system is shown in Figure 3. Additionally, we have explored a more general scenario (mixed tauopathy), where nonuniform parameter values are assigned to different nodes in the brain connectome.

Solutions for the nonfractional model of Equation (6) in the brain connectome. (A) Solutions in each region and (B) the spatial average solution. The fixed parameter values are given in Table 1.

Dual role of astrocytes

2.1.1

As previously stated, astrocytes have a dual function in AD transmission. They work to remove harmful amyloid‐beta and maintain a healthy equilibrium in the brain. However, due to the accumulation of toxic amyloid‐beta, astrocytes become overactive and contribute to disease progression rather than brain protection. Therefore, two scenarios can occur depending on the concentrations of astrocytes present in the brain cells: they can manage the toxic amyloid‐beta, but sometimes they cannot. Here, we capture both cases through our considered network mathematical model in the absence of memory. In our model, the parameter c 1 represents the brain cells’ maximum concentration (carrying capacity) of astrocytes. We consider two different carrying capacities for the astrocytes in the network model, and the average toxic density propagations over time are shown in Figure 4. For c 1 = 0.3, increasing the clearance rate μ reduces the toxic burden on the brain connectome (Figure 4A). This shows that astrocytes can control the brain’s toxic loads. On the other hand, c 1 = 0.1, they fail to manage the proper equilibrium in the brain connectome and encourage a rise in toxic loads (Figure 4B). The nontrivial equilibrium point E ∗ is locally asymptotically stable in both cases.

Spatial average solutions of toxic Aβ (dotted) and toxic τP (solid) for the nonfractional model of Equation (6) for different values of c 1 and μ: (red) μ = 0.1, (green) μ = 0.2 and (blue) μ = 0.3. (A) c 1 = 0.3 and (B) c 1 = 0.1.

Memory effect

2.1.2

Once memory effects become significant, the Markovian framework does not adequately describe the underlying complex dynamic processes behind the progression of neurodegenerative diseases. Here, we analyze the memory effect of AD progression in the brain connectome. The model of Equation (6) has a memory for 0 < α < 1 and memoryless for α → 1, and hence the underlying processes are non‐Markovian for 0 < α < 1. Figure 5 depicts both the toxic propagation over the brain connectome for no‐memory and with memory. In the figure, we plot the spatial average of toxic amyloid‐beta and toxic tau protein. In both cases (α = 0.9 and α = 0.8), the nontrivial equilibrium point E ∗ satisfies the conditions for locally asymptotically stable, as mentioned in Section 4.1. Here, the evolution time of the toxic loads for the fractional model is higher compared to the nonfractional model. Furthermore, with an increase in the memory effect (decreasing the value of α), the evolution time of toxic loads also increases (Figure 5).

Spatial average solutions of toxic amyloid‐beta (dotted), toxic tau protein (solid) and damage (dashed) for the fractional model of Equations (6) with (7) for different values of α with the other fixed parameter values of Table 1 over the brain connectome. In the plot, α = 1, α = 0.9, and α = 0.8 are represented by the red, green and blue curves, respectively.

Neuronal damage

2.1.3

Following the model of Equation (7), the neuronal damage depends on the toxic concentrations present in the brain connectome; hence, the total brain damage depends on the evolutional time (the time required to converge to the stable steady state) of toxic loads. We have mentioned the parameter values directly associated with the neuronal damage in Table 1. These parameter values give us the influence of toxic tau proteins on neural damage and the presence of toxic amyloid‐beta [72, 73, 74, 75]. We plot the spatial average of the damage in Figure 5, and it validates the dependency. It has been observed that the damage converges to its equilibrium point q ∗ = 1 for both fractional and nonfractional models, but the case of the fractional model takes a longer time than the nonfractional model. Overall, memory has a pronounced effect on AD progression.

Mixed tauopathy

2.2

Here, we focus on disease progression for nonuniform parameters over the brain connectome. This is more realistic than the uniform parameters as the presence of heterogeneous density of the ingredients in the brain (e.g., proteins, chemical ions, etc.). We consider the parameter values of b 2 and b 3 from Table 1 in all the brain identities (IDs) except some regions mentioned in references [45, 47]. The methodology of getting these values is mentioned in reference [45]. In the network model, a combination of primary and secondary tauopathies, known as mixed tauopathy, arises because of the nonuniform parameter values in the brain connectome. This causes different stable coexisting steady‐states in the network model, and we divide these into two parts: region ID and region‐wise disease progression.

Region ID‐wise AD progression

2.2.1

The integrated brain connectome data contains 49 brain IDs, each with one or more nodes. We calculate the average concentration of the toxic amyloid‐beta for each brain ID by the formula [71]:

$$[eqn]$$

where Rd is the set of all nodes in that brain ID, and nd denotes the total number of elements in Rd. We use the same formula for the toxic tau proteins and damage dynamics. For the nonfractional model, we observe uniform average concentrations of toxic amyloid‐beta (not shown here) and nonuniform average concentrations of toxic tau proteins along the brain IDs (Figure 6A). This happens due to the direct involvement of the nonuniform parameters b 2 and b 3 in the healthy and toxic tau protein equations. Furthermore, the damage propagation profiles for each brain ID are different (Figure 6B). According to the integrated brain connectome data, the maximum concentrations of toxic amyloid‐beta accumulate in the region ID precuneus, followed by the region IDs left‐putamen, right‐putamen, entorhinal, etc. The damage dynamics show that these region IDs are affected the most at the initial stage of AD progression.

Brain ID‐wise average toxic tau protein propagation (A) and the corresponding neural damage (B) for the nonfractional model. The y‐axis represents the brain‐ID‐wise average toxic tau protein concentrations.

Region‐wise AD progression

2.2.2

We focus on the evolution of the toxic load distributions and their damage profile in seven brain regions (brain stem, frontal, temporal, limbic, basal ganglia, parietal, and occipital), and each region containing one or more brain IDs. The integrated brain connectome data contains Cartesian coordinates for all the nodes in three‐dimensional space and their brain IDs. We plot them according to their regions (mentioned in Table 2) in Figure 7, and in the plot, different colors of the nodes belong to different regions. We have used the same color codes in Figures 7 and 8 for the brain regions.

Three‐dimensional views of the positions of the nodes for the integrated brain connectome data. (A) Axial view and (B) sagittal view. Different colors are used to indicate different brain regions.

Brain region‐wise average toxic tau protein (v˜) propagation for the nonfractional (A) and fractional with α = 0.8 (B) models. The y‐axis represents the average toxic tau protein concentrations.

We apply the Equation (1) to find the average toxic loads for all regions. In this case, the summation is taken over the nodes belonging to the respective regions. We plot the toxic load corresponding to the tau protein for each region in Figure 8. The toxic load for the tau protein converges to different levels for different nodes due to the heterogeneous parameter values in the tau protein equation. According to the integrated data, the occipital region accumulates the most toxic concentration, followed by the parietal, basal ganglia, limbic, temporal, frontal, and brain stem. Moreover, the toxic propagation profile for each region is different. Some regions accumulate more toxic tau protein concentration after the initial progression of the disease but settle down to a comparatively lower concentration for a longer time, for example, basal ganglia, parietal, and limbic. For the other regions, there is not much accumulation in the concentration after the start of the disease; rather, they slowly accumulate the toxic loads and help in disease progression. Figure 8A and B show the toxic tau protein propagation in regions for the traditional nonfractional model and fractional model with α = 0.8, respectively (other parameters are mentioned in the caption). This comparison demonstrates that the memory effect reduces the propagation speed in brain regions. We have observed the fractional model for other values of α (< 1), and the propagation speed decreases with decreasing values of α.

We apply the same Equation (1) to find the regions’ average damage profile. The damage profile for each region is different because of the nonuniform distributions of the toxic tau proteins. The region corresponding to the maximum toxic concentration is damaged first, then the region with the second‐highest concentration, and so on. Figure 9 shows the node‐wise neuronal damage propagation for nonfractional and fractional models. In the results, the damage dynamics are shown till t = 400 (nondimensional time), but we have observed that the required time to damage each region in the brain for the fractional model is higher than the nonfractional model. Hence, the memory effect takes longer to damage the brain cells.

Node‐wise damage propagation (q) for the traditional nonfractional (upper panel) and the fractional with α = 0.8 (lower panel) models in brain connectome. The dark red represents the high damage, and the light yellow represents the low damage.

CONCLUSIONS

3

This work uses a modified heterodimer model to explain interactions between two proteins Aβ and τP. We have incorporated the coupled dynamics of astrocytes dynamics into the modified model and studied the dual role of astroglia before and after AD. Furthermore, we have studied the memory effect in AD progression, which is highly relevant to disease progression. Most of these investigations have been carried out by considering the heterogeneous parameter values, and it is a more realistic synthetic parameter set‐up.

Depending on the activated or deactivated astrocytes, the considered network model shows a dual behavior in disease progression. The density of toxic amyloid‐beta increases as the astrocyte clearance rate increases. For lesser astrocyte densities, however, the reverse scenario occurs. As a result, if enough astrocytes are present in the brain, they can reduce or postpone AD progression; otherwise, they help in AD progression. Furthermore, the fractional differential derivative framework presented here helps to model the memory effect on AD progression. We have shown that an increase in memory (by decreasing the parameter values α) causes a delay in the toxic density propagations in the brain. As a result, it slows down AD progression in the brain, giving experimentalists more freedom in terms of parameters to fit their data appropriately.

We have studied the network model for the parameter values where primary and secondary tauopathy conditions are satisfied in distinct brain areas. This causes a non‐homogeneous distribution of toxic tau proteins in the brain. In addition, the network model demonstrates that nodes with high connectivity have a higher chance of getting the infection and evolving into hubs for spreading disease. Furthermore, different neuronal damage profiles are shown on different brain IDs and in different brain regions. Hence, heterogeneous parameter values in the network model capture a realistic scenario of AD progression in the brain [76]. These nonuniform parameter values in the parameters involved in the amyloid‐beta could be a good extension of this work. The coupling of astrocytes to Aβ and τP represents an advancement in this direction, and one could use this model in an experimental configuration to improve data fitting. Along with the memory, considering heterogeneous parameter values corresponding to amyloid‐beta or astrocytes in different brain IDs or regions is an important avenue for future research on this model. To a greater extent, neurodegenerative diseases involve complex and multiscale processes with multiple levels of biological framework, ranging from molecular and cellular to systemic and even societal. The presented work and the developed methodology allow us to reveal new trends and additional features of the underlying processes. More refined views on the complex dynamics of neurodegenerative diseases are expected with the subsequent incorporation of other scales (e.g., molecular‐to‐cellular and macro‐to‐micro) into the coupled biological framework.

Validating the considered fractional‐order derivative model requires robust datasets and experimental frameworks focusing on Aβ and tau proteins. Some of the potential datasets and experimental frameworks can be used for future validation, and most of them have controlled access. NIAGADS is a collaboration between the National Institute on Aging and the University of Pennsylvania that saves and distributes genetics and genomics data from AD, associated dementias, and aging research to qualified researchers worldwide. The Alzheimer’s Disease Neuroimaging Initiative (ADNI) provides complete data to research the course of Aβ and tau pathology, including imaging, biomarker, and clinical evaluations. The AMP‐AD Knowledge Portal provides free use of multi‐omic data from AD’s patients, including gene expression profiles for Aβ and tau. Analyzing postmortem brain tissue from AD’s patients can also give the model real‐world applicability.

A promising future direction involves the relationship between the memory effect at the molecular level and its phenotypic manifestations at the cellular and tissue levels [77, 78]. For example, biological mechanisms such as synaptic plasticity, protein misfolding, or feedback loops in tau and Aβ propagation may explain the memory‐like features in disease development. Data‐driven simulation of network‐based modeling approaches can offer deeper insights into the spatial and temporal patterns of tau deposition [79, 80, 81]. In addition, data from neurobiological processes such as cellular signaling pathways, neuroinflammation, and neural connections, which are implicated in the spread of various disorders, may represent the model’s history‐dependent progression. The issue of higher‐order interactions in complex networks such as human brain connectome is a very nontrivial task [82, 83]. Nevertheless, focusing on identifying critical fractional orders that best fit experimental data and performing sensitivity analyses to highlight key regulatory factors can give more accurate models for predicting disease progression and developing targeted therapeutic strategies [84].

MATERIALS AND METHODS

4

Alzheimer’s disease strongly correlates with glial cells called astrocytes found in the central nervous system [85]. These astrocytes play a dual role in healthy and AD‐affected brains, and in capturing such dynamics, we introduce an equation corresponding to astrocytes and modify the temporal model defined in Refs. [45, 47] as

$$[eqn]$$

where the subscript t denotes the first‐order ordinary derivative with respect to t, with the initial conditions given by u(0)=u0,u˜(0)=u˜0,v(0)=v0,v˜(0)=v˜0, and w(0)=w0. Here u and v are the healthy densities of Aβ and τP, respectively, and u˜ and v˜ are the toxic densities Aβ and τP, respectively. The parameters a 0 and a 1 denote the mean production and clearance rates of healthy Aβ, respectively, while b 0 and b 1 represent the mean production and clearance rates of healthy τP. The terms a˜1 and b˜1 describe the mean clearance rates of toxic forms of these proteins. The parameters a 2 and b 2 correspond to the mean conversion rates from healthy to toxic proteins. The coupling between the two proteins Aβ and τP is captured by the parameter b 3. Finally, the variable w represents the concentration of activated astrocytes, with c 0 as the production rate and c 0 c 1 indicating the saturation point. The parameter μ is responsible for the dual role of astrocytes. If w > u˜, then astrocytes clear e the concentrations of the toxic amyloid‐beta; otherwise, it helps to increase the toxic concentrations.

The reaction terms on the right side of the Equation (2) determine the substance concentration for any time t > 0. Practically, it means that an individual who has had dementia for 20 years has the same chance of clearing Alzheimer’s as someone who had dementia less than 10 years ago. It is an assumption based on the Markovian process, which is not generally valid. Non‐Markovian processes have been playing an increasingly important role in studying living systems [50, 86, 87], and neuroscience research is no exception where such processes have to be incorporated in state‐of‐the‐art models of neurodegenerative diseases. Clearly, the concentrations of the substances indicated above rely not only on the current time incident t but also on the weighted average concentrations of the pastime range, say [t _ p _, t] for t _ p _ < t. This is commonly referred to as the memory effect [52, 53, 56, 57, 59, 64, 66, 88]. The weight distribution relies on the power of the elapsed time, that is, (t‐t _ p _), and follows the power law correlation function [55, 57]. We can select t _ p _ = 0 without loss of generality. Now, incorporating these into the mathematical Equation (2), we obtain the fractional order differential equations as

$$[eqn]$$

where Dtαz(t) stands for the Caputo fractional derivative, defined as

$$[eqn]$$

and z′ denotes the first‐order ordinary derivative of z. Here in the modified fractional differential Equation (3), the influence of memory decreases when α → 1, and the system tends toward a memoryless system [57, 60]. In addition, the accumulations of toxic amyloid‐beta and tau proteins cause neuronal damage. We consider the memory effect in such neuronal damage equation by modeling it by the following equation:

$$[eqn]$$

with a nonnegative initial condition q (0) = q 0. The case q = 0 signifies a healthy state, that is, neurons are properly functioning, and q = 1 implies an unhealthy or no longer functioning state [45]. In studying neurodegenerative diseases, the development of coupled dynamic models plays a critical role. Different aspects of such coupled models, including those at the neuron‐glial level and the toxic amyloid‐beta dynamics accounting for astrocytes, have been studied in recent papers [39, 53, 69, 89, 90, 91, 92]. The present work is a new step in further refining such coupled models where the Markovian assumption, which cannot be justified in the general dynamic studies of neurodegenerative diseases, is removed. Before proceeding to the analysis of such refined models, we note that such models degenerate into the Markovian case once, in the neural damage equation presented above, the fractional derivative Dtαq is replaced by the ordinary derivative dq/dt [45, 47].

Analysis for the homogeneous system

4.1

Here, we analyze the time‐varying behavior of both the fractional and nonfractional models. First, we describe the equilibria of the nonfractional model of Equation (2) and their stability behaviors. These equilibria for Equation (2) correspond to its time‐independent solutions and can be determined by solving it while setting the time derivatives to zero. In addition, they depend on the parameter values, and we have calculated them numerically. Moreover, each equilibrium point’s stability is determined by the nature of all the eigenvalues of the Jacobian matrix calculated at that point. For any equilibrium point E∗=(u∗,u˜∗,v∗,v˜∗,w∗), the Jacobian matrix of the Equation (2) is given by

$$[eqn]$$

where a11=a0−2a1u∗−a2u˜∗, a12=−a2u∗, a13=a14=a15=0, a21=a2u˜∗, a22=−a˜1+a2u∗−μ(w∗−2u˜∗), a23=a24=0, a25=−μu˜∗, a31=0, a32=−b3v∗v˜∗, a33=b0−2b1v∗−b2v˜∗−b3u˜∗v˜∗, a34=−b2v∗−b3u˜∗v∗, a35=a41=0, a42=b3v∗v˜∗, a43=b2v˜∗+b3u˜∗v˜∗, a44=−b˜1+b2v∗+b3u˜∗v∗, a45=a51=a52=a53=a54=0, and a55=c0−2w∗/c1.

If the real components of all the eigenvalues of J∗ are negative, then the equilibrium point E ∗ is stable; otherwise, it is unstable. In addition, the nonfractional damage equation has only one equilibrium point q ∗ = 1, which is stable. Furthermore, all the equilibrium points for the traditional nonfractional model are also the equilibrium points for the fractional model, but their stability behaviors are not the same for both models. For the fractional model with fixed α, an equilibrium point E ∗ is stable if all the eigenvalues λ _ i _ (i = 1, … ,5) of J∗ satisfy |arg (λ _ i _)| > απ/2; otherwise, it is unstable [66].

Network model in the brain network

4.2

Before going to the brain connectome network model, we extend the temporal Equation (3) into the reaction–diffusion model in a subset of the Euclidean space. This spatial extension is crucial in understanding the spatiotemporal evolution of Aβ and τP in the brain connectome. Indeed, several vivo and vitro studies indicated that the tau protein aggregates and can propagate along synapsis [93]. A spatiotemporal extension of the fractional Equation (3) in a general continuous spatial domain Ω⊂R3 is given by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The first term on the right‐hand side of the first four equations accounts for the random movement of concentrations in the domain Ω. It is assumed that the density of astrocytes is homogeneous in the domain Ω. Here, D1, D˜1, D2, and D˜2 are the diffusion tensors which describe each protein’s spreading. We consider the same damage Equation (4) in this spatial extension, and hereafter, the damage q also depends on the spatial location, that is, q (x, t), x ∈ Ω. The astrocytes also affect the dynamics of the neurons, as they are implicitly involved through toxic amyloid beta.

The main goal is to study the disease progression within the brain connectome. The modified Equation (5) is defined in a continuous domain Ω. Now, we develop a network mathematical model that correlates with the Equation (5) so that we can integrate the brain connectome data [45, 47]. Suppose G is the network brain data with N nodes and E edges. We generate the adjacency matrix A for the graph G, which enables us to build the graph’s Laplacian. The (i, j) (i, j = 1, 2, 3, …, N) element of the matrix A is defined as follows:

$$[eqn]$$

where lij2 represents the mean length squared between the nodes i and j and n _ ij _ is the mean fiber number. Let us now define the elements of the Laplacian matrix L as

$$[eqn]$$

where Dii=∑j=1NAij are the elements of the diagonal weighted‐degree matrix and ρ is the diffusion coefficient. This Laplacian matrix is used to construct a network model for graph G. We employ superscript notations for their respective Laplacian matrices to differentiate the diffusibility of each protein in the brain. In this case, the dynamics of each node j (j = 1, 2, 3, …, N) can be obtained by:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

the corresponding damage equation can be obtained by the fractional differential equation:

$$[eqn]$$

with nonnegative initial conditions. The equilibria of the homogeneous system correspond to the network model’s homogeneous stationary steady‐states of Equation (6).

AUTHOR CONTRIBUTIONS

Swadesh Pal: Conceptualization; investigation; formal analysis; resources; visualization; writing—original draft. Roderick Melnik: Conceptualization; project administration; resources; supervision; writing—review and editing.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

ETHICS STATEMENT

This article does not contain any studies with human or animal materials performed by any of the authors.