Coupled Neural–Glial Dynamics and the Role of Astrocytes in Alzheimer’s Disease

Abstract

Neurodegenerative diseases such as Alzheimer’s (AD) are associated with the propagation and aggregation of toxic proteins. In the case of AD, it was Alzheimer himself who showed the importance of both amyloid beta ($A\beta$) plaques and tau protein neurofibrillary tangles (NFTs) in what he called the “disease of forgetfulness”. The amyloid beta forms extracellular aggregates and plaques, whereas tau proteins are intracellular proteins that stabilize axons by cross-linking microtubules that can form largely messy tangles. On the other hand, astrocytes and microglial cells constantly clear these plaques and NFTs from the brain. Astrocytes transport nutrients from the blood to neurons. Activated astrocytes produce monocyte chemoattractant protein-1 (MCP-1), which attracts anti-inflammatory macrophages and clears $A\beta$. At the same time, the microglia cells are poorly phagocytic for $A\beta$ compared to proinflammatory and anti-inflammatory macrophages. In addition to such distinctive neuropathological features of AD as amyloid beta and tau proteins, neuroinflammation has to be brought into the picture as well. Taking advantage of a coupled mathematical modelling framework, we formulate a network model, accounting for the coupling between neurons and astroglia and integrating all three main neuropathological features with the brain connectome data. We provide details on the coupled dynamics involving cytokines, astrocytes, and microglia. Further, we apply the tumour necrosis factor alpha (TNF-$\alpha$) inhibitor and anti-$A\beta$ drug and analyze their influence on the brain cells, suggesting conditions under which the drug can prevent cell damage. The important role of astrocytes and TNF-$\alpha$ inhibitors in AD pathophysiology is emphasized, along with potentially promising pathways for developing new AD therapies.

1. Introduction

Alzheimer’s disease (AD) is one of the most common late-life dementias, with colossal social and economic impacts. The study by the Institute for Health Metrics and Evaluation published this year in [1] predicts 153 million people will be living with Alzheimer’s disease by 2050. While there are various medical products that help manage the symptoms of AD, as of today, there is only one drug officially approved by the FDA that was designed to treat a possible cause of this form of dementia, rather than the symptoms. Yet, the cost and controversy are limiting the use of this drug, known as aducanuman and marketed as aduhelm. Much of this controversy is related to whether or not the build-up of a protein called amyloid $\beta$ in the brain can be used as a biomarker. According to the AD “amyloid cascade hypothesis”, this build-up causes neurodegeneration, but the link between clearance of amyloid $\beta$ from the brain and deceleration of memory loss and cognitive decline requires further clarification. With over 99% of drugs developed for AD having failed in clinical trials, in addition to more traditional targets related to A$\beta$ and tau proteins, there is an increasing interest in the potential of TNF-$A\beta$ inhibition to prevent AD and improve cognitive function [2].

With underlying interconnections between the processes and factors mentioned above, computational experiments, based on mathematical modelling and computer simulations, can effectively supplement in vivo and in vitro research. In this paper, we present a multiscale model for the onset and evolution of AD that accounts for the diffusion and agglomeration of amyloid beta ($A\beta$) peptide (amyloid cascade hypothesis) and the spreading of the disease through neuron-to-neuron transmission (prionoid hypothesis). Indeed, to cover such diverse facets of AD in a single model, different spatial and temporal scales must be taken into account: microscopic spatial scales to describe the role of the neurons, macroscopic spatial and short temporal (minutes, hours) scales for the description of relevant diffusion processes in the brain, and large temporal scales (years, decades) for the description of the global development of AD. The way in which we combine distinct scales in a single model with brain connectome data assimilation forms the core and major novelty of the paper. Following closely the biomedical literature on AD, we briefly describe the processes that we shall include in our model. In the neurons and their interconnections, several microscopic phenomena take place. We know that $A\beta$ monomers are present in healthy individuals, and therefore, they are unlikely to be toxic. Furthermore, the $\tau$ monomer is non-toxic [3]. On the other hand, $A\beta$ oligomers are highly toxic, playing an important role in the process of cerebral damage, as postulated by the already mentioned amyloid cascade hypothesis.

In the analysis of neurodegenerative diseases, AD in particular, it is important to account for the coupling between neuronal and glial dynamics. Furthermore, given the importance of astrocytes (collectively known as astroglia) in amyloid production [4], several coupled models have been recently developed in this direction, describing Alzheimer’s $A\beta$ accumulation based on calcium-dependent exosome release from astrocytes [5]. It represents a shift from a more traditional view, considering astrocytes as non-excitable brain cells, to a deeper investigation of reactive functions of astrocytes (e.g., increasing the calcium concentration level in response to neurotransmitters and neuromodulators) and their synaptic communication with neurons and other brain cells via what is sometimes labelled as astrocytic networks. Hence, with the ready availability of the data from the brain connectome, derived from various AD mouse models and obtained with the help of transcriptomics and other technologies [6], it is enlightening to go beyond single astrocyte’s consideration and to develop network models allowing such data assimilation (see, e.g., [7,8,9,10] and references therein). This idea is pursued further in this paper.

By now, it is well known that the neuronal and astroglial networks of the brain are innately interwoven, with astrocytes carrying out a multitude of functions in various brain processes, including homoeostasis and neurogenesis, with both positive and negative effects reported [11]. In particular, they are critical in defining the normal operation of the nervous system, but they could also actively contribute to the pathogenesis of AD and other neurodegenerative disorders [8]. As observed in experiments on astroglial atrophy at earlier stages of such neurodegenerative diseases such as Alzheimer’s, Parkinson’s, and various forms of dementia, they play a major role in them, leading to disruptions in synaptic connectivity, disbalance in neurotransmitter homoeostasis, and neuronal death through increased excitotoxicity [12,13]. They maintain their importance in the progression of these diseases at the later stages as well, in particular through their activation and contribution to the neuroinflammatory component of grey matter in pathological neurodegeneration. Given the significance of the contribution of neuroinflammation to Alzheimer’s disease (AD) progression [14], our better understanding and ultimately controlling of these coupled neuronal–astroglial networks become increasingly important, opening the door to developing future therapies. For this to happen, increasing attention is being paid to a relationship between the astrocytes’ effects in the brain and such fundamental processes as synaptic transmission, cognition, and myelination [11]. At the same time, conclusive experimental studies of the role of astrocytes remain extremely challenging, given that the multiple functionalities of these cells are dependent on numerous (and sometimes contradictory) factors during the disease progression. Researchers have shown that both microglia and astrocytes are very heterogeneous in their functions in the diseased brain [15,16,17]. Indeed, on the one hand, they can contribute to the clearance of $A\beta$ and limit the growing inflammation in the brain, while, on the other hand, they may neglect their metabolic role and release neurotoxins, contributing in this way to AD neurodegenerative processes. This leads to a situation where mathematical and computational models developed in a data-driven environment may very efficiently complement the progress made in the experimental domain.

While we briefly touch on other aspects, in this paper, our focus is mainly on the role of astrocytes in Alzheimer’s disease via their dynamic interactions with agglomerations of $A\beta$ peptides. Not only AD is typified by such agglomerations, along with activated glial cells, but also because $A\beta$ plaques trigger intracellular NFT formation, neuronal cell death, neuroinflammation, and gliosis, whereas reactive astrocytes in AD, surrounding these plaques, may additionally contribute to the overall amyloid burden in the brain by secreting $A\beta$ [4]. Indeed, today, we know that a reactive character of astrocytes in AD is usually expressed by intermediate filament proteins and cellular hypertrophy, as well as that these star-shaped glial cells can regulate synaptic communication and modulate brain network functions [7].

The rest of the paper is organised as follows. We develop a network mathematical model for brain connectome data assimilation in Section 2. With the help of brain connectome data, in Section 3, we provide details on two groups of computational experiments elucidating the role of cytokines and astrocytes in AD and giving further details on AD TNF-$\alpha$ inhibitor drugs, quantifying their influence on the reduction of neuronal damage. All numerical results, reported in this section, were obtained with our new network model. Several possible extensions of this work are discussed in Section 4, with concluding remarks given in Section 5.

2. AD Network Model for Brain Connectome Data Assimilation

In this section, we develop a network model based on the consideration originally presented in [18], where a PDE model on AD was discussed. Before going to the full network model in the brain connectome, we first define the diffusion and chemoattraction terms in a network [19,20]. Suppose the network graph $\mathcal{G}$ has V number of nodes and E number of edges. For $j,k=1,2,3,\dots ,V$, the elements of the adjacency matrix ${\mathbf{W}}^l$ corresponding to the graph $\mathcal{G}$ are

$$W_{jk}^l=\frac{n_{jk}}{l_{jk}^2},$$

(1)

where $n_{jk}$ is the mean fiber number and $l_{jk}^2$ is the mean length squared between the nodes j and k. We define a matrix ${\mathbf{L}}^l$ with entries

$$L_{jk}^l=(D_{jj}^l-W_{jk}^l),j,k=1,2,3,\dots ,V,$$

(2)

where $D_{jj}^l={\sum }_{k=1}^V{W}_{jk}^l$. Therefore, at each node j, we take the contribution of the diffusion term for a dummy variable denoted below as u in the following form

$${(\Delta u)}_j=-\sum _{k=1}^V{L}_{jk}^l{u}_k.$$

(3)

Similarly, at each node j, we find the chemoattraction term as

$${(\nabla ·(v\nabla u))}_j=\left(\sum _{k=1}^V{L}_{jk}^c{v}_k\right)\left(\sum _{k=1}^V{L}_{jk}^c{u}_k\right)-v_j\sum _{k=1}^V{L}_{jk}^l{u}_k,$$

(4)

where $L_{jk}^c=(D_{jj}^c-W_{jk}^c)$ with $W_{jk}^c=n_{jk}/l_{jk}$ and $D_{jj}^c={\sum }_{k=1}^V{W}_{jk}^c$.

Now, at the node j in the brain connectome, we are ready to define a network model for Alzheimer’s disease incorporating the astrocytes’ dynamics [18]. We use j as the node index in each of the upcoming equations. Suppose $N_j$ is the density of the living neurons and $N_0$ is the reference density of the neurons in brain cells. Inside the neurons, amyloid beta $A_{\beta }^i$ is constitutively produced from APP at a rate ${\lambda }_{\beta }^i$ and degraded at a rate $d_{A_{\beta }^i}$. In the early stage of disease progression, $A_{\beta }^i$ is overproduced by reactive oxygen species (ROS) factor R. Therefore, the dynamics of $A_{\beta }^i$ is given by

$$\frac{d{A}_{\beta j}^i}{dt}=({\lambda }_{\beta }^i(1+R)-d_{A_{\beta }^i}{A}_{\beta j}^i)\frac{N_j}{N_0}.$$

(5)

The density of extracellular amyloid beta peptides ($A_{\beta }^o$), depends on different factors, such as neuronal death, microglias, astrocytes, etc. The equation for $A_{\beta }^o$ is given by

$$\begin{array}{cc}\hfill \frac{d{A}_{\beta j}^o}{dt}=& A_{\beta j}^i\left|\frac{d{N}_j}{dt}\right|+{\lambda }_N\frac{N_j}{N_0}+{\lambda }_A\frac{A_j}{A_0}\hfill \\ \hfill & -(d_{A_{\beta }^o\widehat{M}}({\widehat{M}}_{1j}+\theta {\widehat{M}}_{2j})+d_{A_{\beta }^{o}M}(M_{1j}+\theta M_{2j}))\frac{A_{\beta j}^o}{A_{\beta j}^o+{\overline{K}}_{A_{\beta }^o}},\hfill \end{array}$$

(6)

where $A_0$ is the reference astrocyte cell density and ${\overline{K}}_{A_{\beta }^o}$ is a Michaelis–Menten coefficient [21]. The first term on the right-hand side of (6) is the contribution due to neuronal death. The second and third terms of (6) are the growths released from amyloid precursor protein (APP) [22] and astrocytes [23], respectively. The last multiplying factor is the clearance of $A_{\beta }^o$ by peripheral macrophages ${\widehat{M}}_1$ and ${\widehat{M}}_2$ and the activated microglias $M_1$ and $M_2$. Here, $0\le \theta <1$ as ${\widehat{M}}_1$ and $M_1$ are more effective in clearing the extracellular $A\beta$ compared to ${\widehat{M}}_2$ and $M_2$. APP on live neurons shed $A\beta$ peptides in both the intracellular and extracellular space [18]. We assumed that most of the $A_{\beta }^o$ is produced from dead neurons, so the production from the live neurons is neglected.

The second most critical factor in AD is the tau protein. Suppose that the tau protein is constitutively produced, and the degradation rates are ${\lambda }_{\tau 0}$ and $d_{\tau }$, respectively. Due to the abnormal concentrations of $A\beta$, i.e., when the production of $A_{\beta }^i$ exceeds a threshold, say $A_{\beta c}^i$, glycogen synthase kinase-type 3 (GSK-3) becomes activated, and it mediates the hyperphosphorylation of tau proteins. Suppose $d_{\tau }$ is the degradation rate of tau proteins due to ROSs. Then, the rate of change of tau protein is given by

$$\frac{d{\tau }_j}{dt}=({\lambda }_{\tau 0}+{\lambda }_{\tau }R-d_{\tau }{\tau }_j)\frac{N_j}{N_0}.$$

(7)

Inside the neurons, NFTs form from the hyperphosphorylation of tau proteins and are released into the extracellular space after the death of the neurons [24,25,26,27]. The equations for NFTs inside the neurons and the extracellular space are given by

$$\frac{d{F}_{ij}}{dt}=({\lambda }_F{\tau }_j-d_{F_i}{F}_{ij})\frac{N_j}{N_0},$$

(8)

$$\frac{d{F}_{oj}}{dt}=F_{ij}\left|\frac{d{N}_j}{dt}\right|-d_{F_o}{F}_{oj},$$

(9)

respectively.

Due to NFTs’ formation in the brain cell, microtubules are depolymerised and destructed, leading to neuron death [24,25,26,27]. Not only NFTs, proinflammatory and anti-inflammatory cytokines are also responsible for neuronal death in the brain. Including these factors in the dynamics for N, we obtain

$$\frac{d{N}_j}{dt}=-d_{NF}\frac{F_{ij}}{F_{ij}+K_{F_i}}{N}_j-d_{NT}\frac{T_{\alpha j}}{T_{\alpha j}+K_{T_{\alpha }}}\frac{1}{1+\gamma I_{10j}{K}_{I_{10}}}{N}_j,$$

(10)

where $T_{\alpha }$ and $I_{10}$ denote the proinflammatory and anti-inflammatory cytokines, respectively.

Astrocytes are primarily activated by the proinflammatory cytokines $T_{\alpha }$ [28], but they are also activated by the extracellular amyloid beta $A_{\beta }^o$ [23]. Therefore, the equation for astrocytes is given by

$$\frac{d{A}_j}{dt}={\lambda }_{A{A}_{\beta }^o}{A}_{\beta j}^o+{\lambda }_{A{T}_{\alpha }}{T}_{\alpha j}-d_A{A}_j.$$

(11)

Microglias and peripheral macrophages clear the NFTs in the extracellular space and keep neurons healthy. Therefore, the dynamics of neuronal death is given by

$$\begin{array}{cc}\hfill \frac{d{N}_{dj}}{dt}=& d_{NF}\frac{F_{ij}}{F_{ij}+K_{F_i}}{N}_j+d_{NT}\frac{T_{\alpha j}}{T_{\alpha j}+K_{T_{\alpha }}}\frac{1}{1+\gamma I_{10j}{K}_{I_{10}}}{N}_j\hfill \\ \hfill & -d_{N_{d}M}(M_{1j}+M_{2j})\frac{N_{dj}}{N_{dj}+{\overline{K}}_{N_d}}-d_{N_d\widehat{M}}({\widehat{M}}_{1j}+{\widehat{M}}_{2j})\frac{N_{dj}}{N_{dj}+{\overline{K}}_{N_d}}.\hfill \end{array}$$

(12)

Amyloid beta oligomers are soluble, and they diffuse in the brain tissue [29,30,31]. Incorporating the diffusion of the oligomers in the network model along with its production (from $A_{\beta }^o$) and degradation, we obtain

$$\frac{d{A}_{oj}}{dt}=-D_{A_o}\sum _{k=1}^V{L}_{jk}^l{A}_{ok}+{\lambda }_{A_o}{A}_{\beta j}^o-d_{A_o}{A}_{oj},$$

(13)

where $D_{A_o}$ is the diffusion coefficient.

In the AD-affected brain, dying neurons produce nonhistone chromatin-associated protein (HMGB-1), and it diffuses in the brain cells [32,33]. The reaction–diffusion equation for the PDE-based model of [18] is simplified into the ODE in the network as follows:

$$\frac{d{H}_j}{dt}=-D_H\sum _{k=1}^V{L}_{jk}^l{H}_k+{\lambda }_H{N}_{dj}-d_H{H}_j.$$

(14)

Microglias travel in the brain cell [34]. Activated microglias are chemoattracted to the cytokines’ high mobility group box 1 (HMGB-1). Furthermore, microglias are activated by the extracellular NFTs and soluble oligomers. The $M_1$ and $M_2$ phenotypes are characterised by the proinflammatory and anti-inflammatory signals from $T_{\alpha }$ and $I_{10}$, respectively. These two types of microglias satisfy the following equations

$$\begin{array}{cc}\hfill \frac{d{M}_{1j}}{dt}=& M_{1j}\sum _{k=1}^V{L}_{jk}^l{H}_k-\left(\sum _{k=1}^V{L}_{jk}^c{M}_{1k}\right)\left(\sum _{k=1}^V{L}_{jk}^c{H}_k\right)-{\lambda }_{M_1{T}_{\beta }}\frac{T_{\beta j}}{T_{\beta j}+K_{T_{\beta }}}{M}_{1j}\hfill \\ \hfill & -d_{M_1}{M}_{1j}+M_G^0({\lambda }_{MF}\frac{F_{oj}}{F_{oj}+K_{F_o}}+{\lambda }_{MA}\frac{A_{oj}}{A_{oj}+K_{A_o}})\frac{\beta {ϵ}_1}{\beta {ϵ}_1+{ϵ}_2},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \frac{d{M}_{2j}}{dt}=& M_{2j}\sum _{k=1}^V{L}_{jk}^l{H}_k-\left(\sum _{k=1}^V{L}_{jk}^c{M}_{2k}\right)\left(\sum _{k=1}^V{L}_{jk}^c{H}_k\right)+{\lambda }_{M_1{T}_{\beta }}\frac{T_{\beta j}}{T_{\beta j}+K_{T_{\beta }}}{M}_{1j}\hfill \\ \hfill & -d_{M_2}{M}_{2j}+M_G^0({\lambda }_{MF}\frac{F_{oj}}{F_{oj}+K_{F_o}}+{\lambda }_{MA}\frac{A_{oj}}{A_{oj}+K_{A_o}})\frac{{ϵ}_2}{\beta {ϵ}_1+{ϵ}_2},\hfill \end{array}$$

(16)

where ${ϵ}_1=T_{\alpha j}/(T_{\alpha j}+K_{T_{\alpha }})$, ${ϵ}_2=I_{10j}/(I_{10j}+K_{I_{10}})$. The parameter $\beta$ is the ratio of the proinflammatory and anti-inflammatory environment, and it determines the relative strengths of $T_{\alpha }$ and $I_{10}$. Here, the ratios $\beta {ϵ}_1/(\beta {ϵ}_1+{ϵ}_2)$ and ${ϵ}_2/(\beta {ϵ}_1+{ϵ}_2)$ in the right-hand sides measure the activated microglias becoming $M_1$ and $M_2$ macrophages, respectively.

Depending on the relative concentrations of $T_{\alpha }$ and $I_{10}$, the incoming macrophages are divided into two phenotypes ${\widehat{M}}_1$ and ${\widehat{M}}_2$ [35]. Furthermore, the phenotype of macrophages ${\widehat{M}}_1$ change to the macrophages ${\widehat{M}}_2$ under the signal $T_{\beta }$. Therefore, the peripheral macrophages satisfy the following equations:

$$\begin{array}{cc}\hfill \frac{d{\widehat{M}}_{1j}}{dt}=& {\widehat{M}}_{1j}\sum _{k=1}^V{L}_{jk}^l{A}_{ok}-\left(\sum _{k=1}^V{L}_{jk}^c{\widehat{M}}_{1k}\right)\left(\sum _{k=1}^V{L}_{jk}^c{A}_{ok}\right)-{\lambda }_{{\widehat{M}}_1{T}_{\beta }}\frac{T_{\beta j}}{T_{\beta j}+K_{T_{\beta }}}{\widehat{M}}_{1j}\hfill \\ \hfill & -d_{{\widehat{M}}_1}{\widehat{M}}_{1j}+\alpha \left(P_j\right)(M_0-{\widehat{M}}_j)\frac{\beta {ϵ}_1}{\beta {ϵ}_1+{ϵ}_2},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \frac{d{\widehat{M}}_{2j}}{dt}=& {\widehat{M}}_{2j}\sum _{k=1}^V{L}_{jk}^l{A}_{ok}-\left(\sum _{k=1}^V{L}_{jk}^c{\widehat{M}}_{2k}\right)\left(\sum _{k=1}^V{L}_{jk}^c{A}_{ok}\right)+{\lambda }_{{\widehat{M}}_1{T}_{\beta }}\frac{T_{\beta j}}{T_{\beta j}+K_{T_{\beta }}}{\widehat{M}}_{1j}\hfill \\ \hfill & -d_{{\widehat{M}}_2}{\widehat{M}}_{2j}+\alpha \left(P_j\right)(M_0-{\widehat{M}}_j)\frac{{ϵ}_2}{\beta {ϵ}_1+{ϵ}_2},\hfill \end{array}$$

(18)

where ${\widehat{M}}_j={\widehat{M}}_{1j}+{\widehat{M}}_{2j}$ and $\alpha \left(P_j\right)=\alpha P_j/(P_j+K_P)$.

$T_{\alpha }$ is produced by proinflammatory macrophages $M_1$ and ${\widehat{M}}_1$. $T_{\beta }$ and $I_{10}$ are produced by $M_2$ and ${\widehat{M}}_2$. Therefore, the equations for $T_{\alpha }$, $T_{\beta }$ and $I_{10}$ are in the form

$$\begin{array}{cc}\hfill \frac{d{T}_{\alpha j}}{dt}=& -D_{T_{\alpha }}\sum _{k=1}^V{L}_{jk}^l{T}_{\alpha k}+{\lambda }_{T_{\alpha }{M}_1}{M}_{1j}+{\lambda }_{T_{\alpha }{\widehat{M}}_1}{\widehat{M}}_{1j}-d_{T_{\alpha }}{T}_{\alpha j},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \frac{d{T}_{\beta j}}{dt}=& -D_{T_{\beta }}\sum _{k=1}^V{L}_{jk}^l{T}_{\beta k}+{\lambda }_{T_{\beta }M}{M}_{2j}+{\lambda }_{T_{\beta }\widehat{M}}{\widehat{M}}_{2j}-d_{T_{\beta }}{T}_{\beta j},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{d{I}_{10j}}{dt}=& -D_{I_{10}}\sum _{k=1}^V{L}_{jk}^l{I}_{10k}+{\lambda }_{I_{10}M}{M}_{2j}+{\lambda }_{I_{10}\widehat{M}}{\widehat{M}}_{2j}-d_{I_{10}}{I}_{10j}.\hfill \end{array}$$

(21)

Activated astrocytes and microglias produce monocyte chemoattractant protein-1 (MCP-1) [36,37], and it is assumed to be of the $M_2$ phenotype. Hence,

$$\frac{d{P}_j}{dt}=-D_P\sum _{k=1}^V{L}_{jk}^l{P}_k+{\lambda }_{PA}{A}_j+{\lambda }_{P{M}_2}{M}_{2j}-d_P{P}_j.$$

(22)

We have used the same estimated parameter values developed by Hao and Friedman [18] for the network model (5)–(22). Further details are provided in

Table 1. Parameter values.

Parameter	Value	Parameter	Value	Parameter	Value
$d_{A_{\beta }^i}$	9.51	$d_{A_{\beta }^i}$	9.51	$K_{M_1}$	0.03
$d_{F_i}$	$2.77\times {10}^{-3}$	$d_{F_o}$	$2.77\times {10}^{-4}$	$K_{M_2}$	0.017
$d_{\tau }$	$0.277$	$d_N$	$1.9\times {10}^{-4}$	$K_{{\widehat{M}}_1}$	0.04
$d_{NF}$	$3.4\times {10}^{-4}$	$d_{NT}$	$1.7\times {10}^{-4}$	$K_{{\widehat{M}}_2}$	0.007
$d_{N_{d}M}$	$0.06$	$d_{N_d\widehat{M}}$	$0.02$	$K_{F_i}$	$3.36\times {10}^{-10}$
$d_A$	$1.2\times {10}^{-3}$	$d_{A_o}$	$0.951$	$K_{F_o}$	$2.58\times {10}^{-11}$
$d_{M_1}$	$0.015$	$d_{M_2}$	$0.015$	$K_M$	0.47
$d_{{\widehat{M}}_1}$	$0.015$	$d_{{\widehat{M}}_2}$	$0.015$	$K_{I_{10}}$	$2.5\times {10}^{-6}$
$d_H$	$58.71$	$d_{I_{10}}$	$16.64$	$K_{T_{\beta }}$	$2.5\times {10}^{-7}$
$d_{T_{\alpha }}$	$55.45$	$d_{T_{\beta }}$	333	$K_{T_{\alpha }}$	$4\times {10}^{-5}$
$d_P$	$1.73$	$A_0$	0.14	$K_{\widehat{M}}$	0.47
$R_0$	6	$M_0$	0.05	$K_{A_o}$	${10}^{-7}$
$M_G^0$	0.47	$N_0$	0.14	$K_P$	$6\times {10}^{-9}$
${\overline{K}}_{A_{\beta }^o}$	$7\times {10}^{-3}$	${\overline{K}}_{N_d}$	${10}^{-3}$

and

Table 2. Parameter values.

Parameter	Value	Parameter	Value	Parameter	Value
$D_{A_o}$	$4.32\times {10}^{-2}$	$D_H$	$8.11\times {10}^{-2}$	$D_P$	$2\times {10}^{-1}$
$D_{T_{\alpha }}$	$6.55\times {10}^{-2}$	$D_{T_{\beta }}$	$6.55\times {10}^{-2}$	$D_{I_{10}}$	$6.04\times {10}^{-2}$
${\lambda }_{\beta }^i$	$9.51\times {10}^{-6}$	${\lambda }_N$	$8\times {10}^{-9}$	${\lambda }_A$	$8\times {10}^{-10}$
${\lambda }_{{\tau }_0}$	$8.1\times {10}^{-11}$	${\lambda }_{\tau }$	$1.35\times {10}^{-11}$	${\lambda }_F$	$1.662\times {10}^{-3}$
${\lambda }_{PA}$	$6.6\times {10}^{-8}$	${\lambda }_{P{M}_2}$	$1.32\times {10}^{-7}$	${\lambda }_{A_o}$	$5\times {10}^{-2}$
${\lambda }_H$	$3\times {10}^{-5}$	${\lambda }_{MF}$	$2\times {10}^{-2}$	${\lambda }_{MA}$	$2.3\times {10}^{-3}$
${\lambda }_{M_1{T}_{\beta }}$	$6\times {10}^{-3}$	${\lambda }_{{\widehat{M}}_1{t}_{\beta }}$	$6\times {10}^{-4}$	${\lambda }_{T_{\beta }M}$	$1.5\times {10}^{-2}$
${\lambda }_{T_{\beta }\widehat{M}}$	$1.5\times {10}^{-2}$	${\lambda }_{T_{\alpha }{M}_1}$	$3\times {10}^{-2}$	${\lambda }_{T_{\alpha }{\widehat{M}}_1}$	$3\times {10}^{-2}$
${\lambda }_{I_{10}{M}_2}$	$6.67\times {10}^{-3}$	${\lambda }_{I_{10}{\widehat{M}}_2}$	$6.67\times {10}^{-3}$	$\theta$	$0.9$
$\alpha$	5	$\beta$	10	$\gamma$	1

with ${\lambda }_{A{A}_{\beta }^o}=1.793$ and ${\lambda }_{A{T}_{\alpha }}=1.54$. We have assumed that these parameter values are uniform for all the nodes in the brain connectome.

3. Numerical Results Based on the Network Model

In this section, we report two groups of computational experiments, focusing on (a) the role of cytokines and astrocytes and (b) the importance of TNF-$\alpha$ inhibitors in reducing neuronal damage in AD. We considered a high-resolution brain connectome structure consisting of $V=1015$ vertices and $E=$ 70,892 edges; the data source is available for the patients’ connectome data at https://braingraph.org (accessed on 20 April 2022) [38]. The network model developed in Section 2 was implemented by using the C programming language and Matlab. We simulated the network model (5)–(22) for each node $j=1,2,\dots ,V$ with uniform initial conditions for all the nodes [18,39]: $A_{\beta }^i={10}^{-6},A_{\beta }^o={10}^{-8},\tau =1.37\times {10}^{-10},F_i=3.36\times {10}^{-10},F_o=3.36\times {10}^{-11},N=0.14,A=0.14,M_1=M_2=0.02,{\widehat{M}}_1={\widehat{M}}_2=N_d=0,H=1.3\times {10}^{-11},T_{\beta }={10}^{-6},T_{\alpha }=2\times {10}^{-5},I_{10}={10}^{-5},P=5\times {10}^{-9}$.

Table 3. The variables of the model and their functions.

Variable	Function	Variable	Function
$A_{\beta }^i$	Amyloid beta inside neurons	$A_{\beta }^o$	Amyloid beta outside neurons
$\tau$	hyperphosphorylated tau protein	$F_i$	Neurofibrillary tangle inside neurons
$F_o$	Neurofibrillary tangle outside neurons	N	Live neurons
A	Astrocytes	$N_d$	Dead neurons
$A_o$	Amyloid beta oligomer	H	High mobility group box 1
$M_1$	Proinflammatory microglias	$M_2$	Anti-inflammatory microglias
${\widehat{M}}_1$	Peripheral proinflammatory macrophages	${\widehat{M}}_2$	Peripheral anti-inflammatory macrophages
$T_{\alpha }$	Tumour necrosis factor alpha	$T_{\beta }$	Transforming growth factor beta
$I_{10}$	Interleukin 10	P	Monocyte chemoattractant protein-1

lists all the variables used in the model, and the units of all these variables are given in g/mL. We used the value of R as

$$R=\left\{\begin{array}{cc}{R}_{0}t/100\hfill & 0\le t\le 100,\hfill \\ R_0\hfill & t>100.\hfill \end{array}\right.$$

3.1. Computational Experiments on the Role of Cytokines and Astrocytes in AD

After integrating the brain connectome data, we computed all the components involved in Equations (1)–(4). We plot the average densities of twelve variables for the network model in Figure 1. These variables have an influence on the model, but we mainly focused on the astrocyte and microglia variables and the other variables directly associated with these two. Extracellular $A\beta$ and the pro-inflammatory cytokines $T_{\alpha }$ play a crucial role in the growth of astrocytes in the brain (see Equation (11)).

The density change in the astrocytes or the microglia causes a change in the MCP-1’s density, following from Equation (22). Figure 2 shows the density change in astrocytes and MCP-1 by changing the growth parameters in astrocytes, while the other concentrations do not change as such. As additional factors have been taken into account, here we obtained a different result as compared to [18]. As time progresses, the densities of extracellular $A\beta$, $A\beta$-oligomers, astrocytes, and MCP-1 decrease.

3.2. Computational Experiments on AD Drugs

Next, our attention is drawn to soluble inflammatory cytokines. It is well known that soluble cytokine receptors regulate inflammatory and immune events by functioning as antagonists of cytokine signalling. Among various biologic medical products that, by interfering with tumour necrosis factor (TNF), are used to treat autoimmune diseases, caused by an overactive immune response, etanercept (marketed as enbrel) is quite popular. In its essence, it is an amalgamation protein, produced by recombinant DNA, that fuses the TNF receptor to the constant end of the immunoglobulin G1 antibody. Etanercept acts as a TNF-$\alpha$ inhibitor, where TNF-$\alpha$ is considered to be the key regulator of the inflammatory/immune response in many organ systems. We used the amount of etanercept as a parameter in a series of computational experiments that will be discussed here.

Before going to the numerical results involving etanercept, we provide our motivation for focusing on TNF-$\alpha$ and further insight into its leading role in AD pathophysiology. First, we recall that three distinctive neuropathological features of AD are: (a) extracellular deposits of A$\beta$ peptides assembled in plaques, (b) intraneuronal accumulation of hyperphosphorylated tau proteins forming tangles, and (c) chronic inflammation [40]. While (a) and (b) we addressed earlier, here, we note that as far as (c) is concerned, the pro-inflammatory cytokine TNF-$\alpha$ plays a critical role. Moreover, with existing evidence, indicating that TNF-$\alpha$ signalling frequently makes pathologies related to (a) and (b) worse, a growing interest is seen in modulating this signalling and developing anti-TNF-$\alpha$ AD therapies, allowing improved cognitive performance. Compared to other approaches in developing AD treatments [41], TNF-$\alpha$ inhibitors have been consistently rated favourably. While the full range of pathogenetic mechanisms underlying neuronal death and dysfunction in AD remain unclear, most recent analyses convincingly imply that TNF-mediated neuroinflammation is linked to AD neuronal necroptosis [42]. Furthermore, given the pathophysiological importance of the entire TNF-TNFR1/2 system, more and more attention is currently paid to its other components as well, in particular to tumour necrosis factor receptor 2 (TNFR2) of the cytokines, which promotes neuronal survival downstream. While both TNFR1 and TNFR2 can induce pro-inflammatory activities, it is TNFR2 that can also elicit strong anti-inflammatory activities and has protective effects. Recent studies (e.g., [43]) indicate that the TNF pathway can contribute to resilience in AD. The latter concept is important in understanding heterogeneity in cognitive and behavioural phenotypes of AD, which requires involvement not only of A$\beta$ and tau proteins, but other molecular factors as well. This leads, among other things, to the investigation of genetic variants of the TNFR2 pathway as a marker of resilience and the TNFR2 pathway itself as a target for developing new AD therapies [43].

As we already mentioned earlier, debates over AD drugs continue today, with the first drug able to remove amyloid approved only in 2021 (it is also the first new AD drug approved since 2003). The controversy around this new drug, aducanumab, is effectively centred on AD biomarkers and whether the extent of amyloid plaques can be considered as one of them because some scientists believe that they are more like a side-effect of the disease process. In the meantime, this controversy has generated a burst of new research activities and the development of another drug, known as donanemab, which is currently in late-stage clinical trials. Considered to be an important advance in amyloid pathology, it is expected to be able to treat early symptoms of AD. In the meantime, scientists are in agreement that new treatments, drugs, and therapies are urgently needed [44], and mathematical modelling and computational experiments will be playing an increasingly important role in these new developments.

In what follows, we account for the fact that many clinical trials of drugs aimed at preventing or clearing the $A\beta$ and tau protein pathology have failed to demonstrate efficacy and that one of the possible treatments could be based on TNF-$\alpha$ inhibitors (suggested also in [18]). For the treatment, we first ran the model for what corresponds to 300 days in order to ensure that AD has been diagnosed, and then, we applied continuous treatment by the drug from Day 300 until the end of 10 years. In this case, we can replace Equation (19) with

$$\frac{d{T}_{\alpha j}}{dt}=-D_{T_{\alpha }}\sum _{k=1}^V{L}_{jk}^l{T}_{\alpha k}+{\lambda }_{T_{\alpha }{M}_1}{M}_{1j}+{\lambda }_{T_{\alpha }{\widehat{M}}_1}{\widehat{M}}_{1j}-d_{T_{\alpha }}{T}_{\alpha j}-f{T}_{\alpha j},$$

(23)

where f is proportional to the amount of etanercept. We simulated this equation with the full network model for the brain connectome and plotted the result in Figure 3. For this set of computations, we took $f=10d_{T_{\alpha }}$ along with ${\lambda }_{A{A}_{\beta }^o}=1.793$ and ${\lambda }_{A{T}_{\alpha }}=1.4$. After applying the drug, we observed reduced neuronal damage in the brain (see the middle-top sub-figure in Figure 3).

The drug aducanumab is considered to be one of the most effective in clearing $A\beta$. In this case, we replaced Equation (6) with

$$\begin{array}{cc}\hfill \frac{d{A}_{\beta j}^o}{dt}=& A_{\beta j}^i\left|\frac{d{N}_j}{dt}\right|+{\lambda }_N\frac{N_j}{N_0}+{\lambda }_A\frac{A_j}{A_0}\hfill \\ \hfill & -(d_{A_{\beta }^o\widehat{M}}({\widehat{M}}_{1j}+\theta {\widehat{M}}_{2j})+d_{A_{\beta }^{o}M}(M_{1j}+\theta M_{2j})(1+g))\frac{A_{\beta j}^o}{A_{\beta j}^o+{\overline{K}}_{A_{\beta }^o}},\hfill \end{array}$$

(24)

where g is proportional to the amount of the dosing level of the drug aducanumab. We simulated the entire network model (5)–(22) with (24) with $g=10$, ${\lambda }_{A{A}_{\beta }^o}=1.793$, and ${\lambda }_{A{T}_{\alpha }}=1.54$. Figure 4 depicts the effect of the drug in $A\beta$ aggregation in the advancement of time. It has a pronounced effect only on the extracellular $A\beta$ concentrations, not on neuronal death. However, the TNF-$\alpha$ inhibitors have a strong effect of reducing neuronal death. These results agree with those obtained in [18].

4. Discussion and Future Directions

There are several recently developed models, dealing with astrocytes, that consider averaged characteristics such as the collective exosomal release rate in astrocytes. Such models have an advantage of their generalisations to account for temperature effects, which is an important consideration, given the recent discovery of an intrinsic connection between the temperature dependence of exosome release and $A\beta$ neurotoxicity [5]. Among other things, such models can describe the synapse and astrocyte couplings and allow replicating typical calcium oscillations in astrocytes under the influence of $A\beta$. Therefore, it would be instructive to extend the network model proposed here to account for thermal effects.

Tau proteins play a more prominent role than the amyloid hypothesis suggests. The $\tau$Ps are usually considered as secondary agents in the disease even though: (i) other $\tau P$-related diseases (tauopathies), such as frontotemporal lobar degeneration, are mostly dominated by tau spreading; (ii) brain atrophy in AD is directly correlated with large concentrations of NFT; (iii) the $\tau P$ distribution determines disease staging, and lowering $\tau P$ levels prevent neuronal loss; (iv) $\tau P$ reduces neural activity and is the main factor associated with cognitive decline. This motivates another possible extension of the developed model by combining the tau protein dynamics more precisely [19,45], rather than considering only a linear contribution of the tau protein. In addition, the dynamics of the variables in different regions in the brain connectome would give a better understanding of the disease progression [45]. Following this recent study, we note that although we chose here uniform parameter values all over the brain connectome, this can be extended to different parameter values in respective regions according to the clinical data.

Regarding the route connected with TNF inhibitors, in addition to the possible extensions already mentioned in the previous section, we will also mention that one of the existing difficulties lies with the fact that classical biologic TNF-$\alpha$ inhibitor macromolecules cannot cross the blood–brain barrier [46,47]. This requires the development of blood–brain-barrier-penetrating TNF-$\alpha$ inhibitors, and from a modelling point of view, further extensions of the models developed here may be needed to qualitatively estimate this factor.

5. Conclusions

We constructed a network model to study neurodegenerative disorders in the brain connectome, focusing on Alzheimer’s disease. The developed model can capture the concentrations of the variables in different regions of the brain connectome, which could not be identified by earlier-developed simple PDE-based models. All three distinctive neuropathological features of AD, including amyloid beta and tau proteins, as well as neuroinflammation were considered in the network model for brain connectome data assimilation. Special attention was given to the role of cytokines and astrocytes, as well as to the influence of anti-$A\beta$ and TNF-$\alpha$ inhibitor drugs in AD pathophysiology. We showed that etanercept has good efficacy in most of the aspects, including neuronal death, while aducanumab has a good efficacy only in reducing the aggregation of extracellular amyloid beta. Among other applications, one may choose the developed methodology to address the diffusion and chemoattraction challenges by evaluating the corresponding term’s contributions in the network model. Finally, potentially promising pathways for developing new AD therapies were also discussed.