# Quantifying uncertainty and sensitivity in an Alzheimer's disease model: A mathematical approach

Mitali Maji, Laurent Pujo-Menjouet, Subhas Khajanchi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Mitali Maji, Laurent Pujo-Menjouet, Subhas Khajanchi. Quantifying uncertainty and sensitivity in an Alzheimer's disease model: A mathematical approach. 2025. hal-05360190

HAL Id: hal-05360190

https://hal.science/hal-05360190v1

Preprint submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Quantifying uncertainty and sensitivity in an Alzheimer's disease model: A mathematical approach

Mitali Maji<fnref n="1" />*, Laurent Pujo-Menjouet<fnref n="2" />,†, Subhas Khajanchi³‡

## Abstract

To understand the dynamics of Alzheimer's disease, we formulate a generalized mathematical model based on three events: aggregation of disease-related proteins, activation of immune cells and initiation of inflammation. We incorporate functional forms in the model to represent the complex biological interactions between components related to Alzheimer's disease. We take explicit forms depending on the properties of functions in the model. We describe the system dynamics by locating biologically feasible steady states, determining stability properties and identifying the effective parameters. Parameters are estimated using two methods: biological literatures and data fitting. We perform sensitivity and uncertainty analyses to identify the most influential parameters. Partial Rank Correlation Coefficient and scatter plots are used to visualize global sensitivity. Our results reveal that lower activation rate and higher proliferation rate of microglia may contribute to a reduction in toxic protein aggregate levels, thus slowing the disease progression.

**Keywords:** Neuroinflammation, Amyloid beta oligomers, Tau protein, NFTs, Microglia

## 1 Introduction

One of the most prevalent forms of dementia is Alzheimer's disease (AD), which is a progressive neurological disorder that affects memory, thinking, personality and behavior. Approximately 55 million people worldwide are affected by the disease currently and the number of cases continues to rise due to a lack of effective cures and preventive measures. Mostly, AD ranks as the seventh leading causes of death globally [1].

AD is triggered by the abnormal accumulation of, on one hand, amyloid beta ($A\beta$) protein in the brain, forming extracellular $A\beta$ plaques, and on the other hand, tau ($\tau$) proteins aggregating intra-cellularly into neurofibrillary tangles (NFTs). This affects neuronal functioning and connectivity, leading to progressive loss of brain function and the affected brain areas include the hippocampus and cerebral cortex.

<footnote n="1">¹,³ Department of Mathematics, Presidency University, 86/1 College Street, Kolkata, India.</footnote>

<footnote n="2">² Universite Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, 69622 Villeurbanne, France</footnote>

<footnote>*Email: mitali.rs@presiuniv.ac.in</footnote>

<footnote>†Email: pujo@math.univ-lyon1.fr</footnote>

<footnote>‡Email: subhas.maths@presiuniv.ac.in</footnote>
AD is a multifactorial condition and the application of mathematical models represents just one facet of its comprehensive study. Developing mathematical models can greatly enhance our understanding of the fundamental causes of AD. Numerous mathematical models have already been investigated to understand disease progression, and many studies are currently ongoing. The work by Hao and Friedman [2] presents a mathematical model of AD involving the key elements, including neurons, astrocytes, microglia, macrophages, cytokines, Aβ aggregation and τ protein dynamics. The model simulates the effects of various drugs and suggests that combined therapy with tumor necrosis factor-α (TNF-α) inhibitors and anti-Aβ drugs could be more effective in slowing disease progression. In neurodegenerative diseases, protein-protein interactions play a critical role. Nonlocal continuous and network models have been suggested due to the importance of concentration variations in Aβ and τ proteins in the onset and development of AD [3, 4]. A temporal mathematical model involving the key components for AD-like microglia, astrocytes, neurons and Aβ was given by Puri and Li in their work [5]. In the study of Craft et al. [6], authors designed a mathematical compartmental model to explain the accumulation of the Aβ protein within the brain, cerebrospinal fluid (CSF) and plasma. To get a deeper understanding of the dynamics of oligomer formation and its connection to the onset of AD, Lindstrom et al. [7] created a mathematical model of the oligomerization process. Ciuperca et al. [8] developed a novel mathematical model for the development of AD, taking into account the interactions between Aβ monomers, oligomers, plaque, microglia and interleukins. A series of work has been done on AD with different types of mathematical approaches that are reflected in the review article by Maji and Khajanchi [9].

This work focuses on a mathematical model that involves the complex interactions between disease-related proteins and glial cell. Since the parameter values for the model were obtained from various sources, sensitivity analysis (SA) and uncertainty analysis (UA) have been conducted to find most sensitive parameters. SA and UA are two important methods to understand the correlation between the model's output and input. To enhance this understanding, we apply SA and UA to our complex model. UA and SA both are related to each other but they are not identical. UA generally examines the variability or variations of model output in the presence of a known distribution in input parameters. Hence, this uncertainty in the inputs limits in getting a specific conclusion to the model's output. These imperfections are measured with UA and treats the simulation result as probabilistic rather than deterministic [10]. One of the most well-known techniques to do this UA is Monte Carlo (MC) simulation which is a sample-based method. There are several sampling techniques to extract parameters from probabilistic distributions, including Simple MC sampling, Latin hypercube sampling (LHS) and Importance sampling. In this work, we use the LHS method in our simulation. The LHS method, a stratified sampling approach within MC framework, was first introduced by McKay in 1979 [11]. On the other hand, SA is a method that determines how changes in the input affect the output and identifies which input uncertainties have the most significant impact. SA involves performing multiple simulations by varying the input around a small value. Global SA assesses the impact of varying all input parameters simultaneously over their entire range of values [12, 13]. Global SA techniques are commonly implemented using MC simulations (sampling-based methods). A useful initial approach when applying sampling-based methods is to utilize scatter plots.

Previous models of AD focus on either Aβ or τ pathology separately, neglecting the complex interactions between these proteins and the immune response, such as microglial activation and cytokine dynamics. Few studies have been done on the role of neuroinflammation and the dual nature (protective and harmful) of microglial states, highlighting a notable gap in the existing research. The main objective of this study is to introduce and investigate a biological-based mathematical model for AD, which involves the most important disease-related proteins and immune responses.
Next, we improve the model’s predictive power by fitting the data to specific model components, which allows us to have more precise simulations, and to better understand the disease’s progression and potential treatment outcomes.

The paper is organized as follows. Section 2 presents the formulation of the mathematical model through a system of differential equations with the elements in functional forms. We provide some basic properties of the model. We determine the steady states and do the stability analysis in Section 3. Based on the properties of functions in the model, we put the explicit forms to the functions in Section 4. Section 5 covers the numerical simulation of the analytical results, along with parameter estimation, SA and UA. Finally, a brief discussion and conclusion are provided in Section 6.

# 2 The model

## 2.1 Biological Background

AD is a complex disease with multiple contributing factors. This disease can be characterized by the accumulation of protein, activation of glial cells, production of inflammation or oxidative stress, brain damage, neuronal loss, etc. [14, 15]. Here, we mainly focus on three important factors involved in AD: Aβ protein, τ protein and a glial cell. The following portion provide details descriptions of the contributions, background and interconnections of proteins, microglia and inflammation.

*   **Aβ protein:** Aβ protein is a small peptide derived from the larger amyloid precursor protein (APP). The Aβ protein is typically composed of 36 to 43 amino acids. Aβ42 is a specific form of the Aβ peptide, consisting of 42 amino acids. Aβ42 is considered particularly neurotoxic as it has high tendency to aggregate more readily than other forms like Aβ40 [16]. It can form soluble oligomers that contribute to the formation of Aβ plaques observed in AD [17]. Aβ oligomers are small, soluble aggregates of Aβ peptides, while Aβ plaques are large, insoluble aggregates of Aβ peptides that form in the extracellular spaces of the brain [18]. Aβ oligomers are considered the most toxic form of Aβ [19]. Numerous studies [6, 20] show that Aβ aggregation and its processes strongly impact AD. Aβ proteins and their polymer units, which include Aβ monomers, dimers, polymers, oligomers, fibrils and plaques, are considered contributing factors in the progression of AD. It is believed that Aβ plaques are toxic to the brain and are among the primary factors responsible for the advancement of the disease. But recent studies show that the Aβ oligomers are more toxic to the brain than plaques [21]. Hence, we give importance to the Aβ monomer and oligomers in this work.

*   **τ protein:** τ protein is also critical to the pathogenesis of AD. τ is a microtubule-associated protein that helps neurons to maintain their internal structure. When τ undergoes aberrant chemical changes, this causes τ proteins to separate from microtubules and form insoluble tangles inside neurons, known as neurofibrillary tangles (NFTs).

τ protein hyperphosphorylation is an important phase in the development of AD that is strongly related to the influence of Aβ oligomers [22]. Indeed, these oligomers interact with neuron surface receptors and activate intracellular signaling cascades and kinases such as glycogen synthase kinase-3 β (GSK-3β) and cyclin-dependent kinase-5 (CDK-5) [23, 24]. These kinases phosphorylate τ proteins. As a result, hyperphosphorylated τ clumps into insoluble fibrils that form NFTs.
* **Microglia and inflammation:** Rather than proteins, neuron-supporting cells, called glial cells, contribute significantly to disease progression through their roles in neuroinflammation and immune response. Several types of glial cells, including astrocytes, microglia and oligodendrocytes are implicated in the disease progression. Microglia are the most prevalent immune cells in the brain [25]. Cytokines are signaling proteins released by microglia.

Finally, microglia and neuroinflammation hold a significant place in the pathology of AD. There is still a question about the protectiveness and harmfulness of microglia. In AD, microglia have a dual role in disease progression, as they are crucial for clearing $A\beta$ deposits but also contribute by releasing cytotoxic mediators [26]. The presence of $A\beta$ accumulation and NFTs triggers the activation of microglia [27, 28]. In the early stage of AD, activated microglia behave like neuroprotective and remove harmful aggregates using a process known as phagocytosis [29]. This early reaction protects the brain by removing toxic $A\beta$. However, microglia do not efficiently clear $A\beta$, as a result, $A\beta$ protein again starts accumulating leading to prolonged activation of these immune cells. This long-term activation causes microglia to become pro-inflammatory, releasing cytokines such interleukin-1$\beta$ (IL-1$\beta$), TNF-$\alpha$ and interleukin-6 (IL-6) [29, 30]. These cytokines contribute to neuroinflammation, which may cause neuronal and synaptic damage, and also accelerate cognitive decline. On the other hand, these pro-inflammatory cytokines play a crucial role in AD by promoting the aggregation of $A\beta$ monomers. Neuroinflammatory cytokines can increase amyloid precursor protein (APP) levels and unregulate the activity of $\beta$-secretase and $\gamma$-secretase [31]. This leads to an increase of the production of $A\beta$ monomer [32].

## 2.2 Model development

The development of our model begins by identifying the key biological components related to AD, including $A\beta$, $\tau$ protein and immune responses such as microglial activation and cytokine production, degradation and interaction dynamics of these proteins as well as the role of neuroinflammation in disease progression. We include below the underlying mechanism into the mathematical model:

* $A\beta$ proteins simulate their own production, aggregation and clearance, eventually forming $A\beta$ oligomers,

* $\tau$ proteins misfold and potentially interact with $A\beta$ proteins,

* misfolding of $\tau$ proteins results in formation of NFTs,

* microglial activation is modeled as part of the immune response, where activated microglia influence the clearance of $A\beta$ and produce cytokines,

* microglia show both protective and neurotoxic role reflecting its complex behavior in AD progression, regardless the stages of AD progression,

* production and release of cytokines occur by activated microglia, which influence the inflammatory environment,

* cytokines help in the production of $A\beta$ monomer.

The above underlying mechanism in AD are schematically represented through a diagram given in Figure 1.
<figure><img src="image_3.png" /><figcaption>Figure 1: A schematic diagram for the progression of AD to reflect the model (1) that involves the following components: Aβ monomers, Aβ oligomers, τ protein, NFTs, microglia and cytokines. Explanation of the schematic diagram: 1. Aβ oligomers are aggregated from Aβ monomers, 2. τ proteins are phosphorylated by Aβ oligomers through activation of GSK-3, 3. The phosphorated τ proteins accumulate and form NFTs, 4. Microglia are activated by Aβ oligomers and NFTs, 5. Activated microglia produce cytokines, 6. Aβ monomers are produced by cytokines through increasing APP levels, neuronal stress due to Aβ oligomers is not explicitly considered, 7. Aβ oligomers are cleared from the brain by early-stage activated microglia.</figcaption></figure>

Based on these events, we formulate a system of differential equations to describe the interactions between these proteins and the immune system.

## 2.3 Model formulation

A system of six ordinary differential equations is employed to describe our proposed dynamics of AD. The considered model populations at time $t \ge 0$ are given by:

* $A_m(t)$ denotes the concentration of Aβ monomers;

* $A_β(t)$ represents the concentration of Aβ oligomers;

* $\tau_p(t)$ corresponds to the concentration of τ proteins;

* $N_F(t)$ signifies the concentration of NFTs;
* $M(t)$ indicates the density of microglia;

* $I(t)$ refers to the concentration of cytokines.

The mathematical model reads as

$$ 
\left\{
\begin{aligned}
\frac{dA_m}{dt} &= A_m r(A_m) + A_m f(I) - h(A_m), \\
\frac{dA_\beta}{dt} &= h(A_m) - s(M)A_\beta - \delta_\beta A_\beta, \\
\frac{d\tau_p}{dt} &= \gamma_p A_\beta - \delta_p \tau_p - g(\tau_p), \\
\frac{dN_F}{dt} &= g(\tau_p) - \delta_N N_F, \\
\frac{dM}{dt} &= \gamma_M - \delta_M M + (f_M(A_\beta) + g_M(N_F))(\tilde{M} - M)M, \\
\frac{dI}{dt} &= h_I(A_\beta)M - \delta_I I,
\end{aligned}
\right.
\qquad (1) $$

with the initial conditions $A_m(0) \ge 0$, $A_\beta(0) \ge 0$, $\tau_p(0) \ge 0$, $N_F(0) \ge 0$, $M(0) \ge 0$ and $I(0) \ge 0$. The descriptions of parameters and functions in the above model (1) are given below:

* $r(A_m)$ is the growth rate of Aβ monomer;

* $f(I)$ describes the production activation of Aβ monomer by cytokines;

* $h(A_m)$ represents the formation of Aβ oligomers from Aβ monomers;

* $s(M)$ is the clearance rate of Aβ oligomers by activated microglia;

* $\gamma_p$ denotes the phosphorylation of $\tau$ protein by Aβ oligomers;

* $g(\tau_p)$ signifies the formation of NFTs from $\tau$ proteins;

* $\gamma_M$ is the production rate of microglia;

* $f_M(A_\beta)$ and $g_M(N_F)$ denote the activation of microglia by Aβ oligomers and NFTs, respectively. With this, it is considered that surrounding the neuron, the number of activated microglia cannot exceed a certain threshold $\tilde{M}$;

* $h_I(A_\beta)$ describes the production of cytokines that depends on the concentration $A_\beta$;

* $\delta_\beta$, $\delta_p$, $\delta_N$, $\delta_M$ and $\delta_I$ are the degeneration rates of Aβ oligomers, $\tau$ protein, NFTs, microglia and cytokines, respectively.

## 2.4 Model assumptions

To investigate the dynamics of system (1), we consider the subsequent biologically takeable assumptions delineating the characteristics of the functions mentioned in the system (1). All the functions $r$, $f$, $h$, $s$, $g$, $f_M$, $g_M$ and $h_I$ are supposed continuously differentiable functions.

* **Properties of the growth rate $r(A_m)$:** $r(A_m)$ needs to satisfy $r(0) = r_m > 0$ and $r(k_m) = 0$, where $k_m > 0$ is the carrying capacity of Aβ monomers in absence of any inflammation. Also, it is decreasing on $[0, k_m]$.
* **Properties of $f(I)$**: inflammation plays an important role in the activation of Aβ monomer production. We consider the function $f$ in such a way that if the population of cytokines becomes too large, the production of Aβ monomer cannot exceed infinity; $f$ is increasing, $f(0) = 0$ and $\lim_{I \to \infty} f(I) = \mu_I > 0$.

* **Properties of $h(A_m)$**: the expression $h(A_m)$ denotes the process of oligomerization, where monomers combine to form oligomers. In other words, the formation of oligomers needs at least two monomers to come together. Hence, we consider that $h$ is increasing and $h(0) = 0$. Moreover, $\lim_{A_m \to \infty} h(A_m) = +\infty$.

We can even determine the exact form of $h(A_m)$ by the monomer-oligomers reaction process. Let us assume, as in [33, 34, 35, 36], that oligomers are produced by the following $n-1$ reactions from monomers in such a way that

$$R_i + A_m \rightleftharpoons R_{i+1}, \quad (2)$$

for $i = 1, 2, ..., (n-1)$. $R_i$ represents the aggregates of i-monomers for $i = 2, ..., n$, where $n \in \mathbb{N}$, and the reaction rate of each equation is given by $r_i$. By taking all the equations at the equilibrium, we get $R_2 = r_2A_m^2$, $R_3 = r_3R_2A_m = r_3r_2A_m^3$, and so on with $R_n = \mu_m A_m^n$, where $\mu_m$ is the formation rate of oligomers and it is given by $\mu_m = r_1r_2...r_n > 0$. We set $h(A_m) = \mu_m A_m^n$ with $n > 1$.

* **Properties of $s(M)$**: microglia increases the clearance efficiency of Aβ oligomers up to some saturation limit, that is, $s$ is an increasing function, followed by $s(0) = 0$, $\lim_{M \to \infty} s(M) = \tilde{\mu}(> 0)$, a constant, and $s'$ is decreasing.

* **Properties of $g(\tau_p)$**: the phosphorylated $\tau$ proteins assemble into protofibrils, called paired helical filaments (PHFs), and then into NFTs. Using a similar concept [37, 38] to model the formation of protein monomers into oligomers, protofibrils, and fibrils (e.g., Aβ protein), we define the aggregation term as $g(\tau_p) = \mu_p \tau_p^m$, where $m > 1$ represents the number of $\tau_p$ monomers required to form fibrils and $\mu_p$ is the formation rate of fibrils.

* **Properties of $f_M(A_\beta)$**: as the concentration of Aβ oligomers increases, it implies more activation of microglia. We take $f_M$ in such a way that, it is increasing and $f_M(0) = 0$. We assume that the microglial response can accelerate or remain constant with increasing Aβ oligomers with no saturation limit. Hence, $f'_M$ is nondecreasing.

* **Properties of $g_M(N_F)$**: since NFTs are implicated in inflammation and immune response, more NFTs lead to greater microglial activation. $g_M$ is an increasing function with $g_M(0) = 0$. Here also we assume no biological limits to the activation of microglia is considered and $g'_M$ is nondecreasing.

* **Properties of $h_I(A_\beta)$**: the production of cytokines is influenced by the interaction between microglia and Aβ oligomers. The presence of Aβ oligomers stimulates microglia to produce cytokines up to a certain limit, that is, $h_I$ is increasing with $\lim_{A_\beta \to \infty} h_I(A_\beta) = \mu_\beta$. Moreover, $h_I(0) = 0$ and $h'_I$ is decreasing.

Let us give a positivity and boundedness results for the solutions, ensuring that cell populations remain positive and do not grow unboundedly.
**Definition [39, 40]:** Let us assume Y be a set in $\mathbb{R}^n$. A dynamical system is given by

$$ \dot{x} = f(t, x), \quad t \ge 0, $$

where $x(t) \in \mathbb{R}^n$ are the state variables, $t$ is the time variable, and $f(t, x)$ is a real-valued continuous function. The set Y is called a *positively invariant set* of the dynamical system if $x(0) \in Y$ implies $x(t) \in Y$ for all $t \ge 0$.

**Proposition [39, 40]:** Let the nonnegative quadrant $\mathbb{R}_+^n$ be represented as $\{x \in \mathbb{R}^n | x_i \ge 0, i = 1, 2, \dots, n\}$. Then $\mathbb{R}_+^n$ is an invariant set with respect to the above dynamical system if any point $x$ on the boundary of $\mathbb{R}_+^n$ holds the following condition

$$ f_{ij}(t_0, x) \ge 0, \quad j = 1, 2, \dots, k, $$

where $x_{ij} = 0$ for $j = 1, 2, \dots, k$.

**Statement 1:** A nonnegative set $\mathbb{R}_+^6 = \{(A_m, A_\beta, \tau_p, N_F, M, I) \in \mathbb{R}^6 | A_m \ge 0, A_\beta \ge 0, \tau_p \ge 0, N_F \ge 0, M \ge I \ge 0\}$ is an invariant set with respect to the dynamical system (1) if any point $(A_m, A_\beta, \tau_p, N_F, M, I)$ on the boundary of $\mathbb{R}_+^6$ holds the following condition

$$ f_i(0, A_m, A_\beta, \tau_p, N_F, M, I) \ge 0, $$

for $i = 1, 2, \dots, 6$, where $A_m = 0, A_\beta = 0, \tau_p = 0, N_F = 0, M = 0, I = 0$, respectively.

**Statement 2:** All the solution of the system (1) are bounded.

Proofs of positivity and boundedness for the solutions are given in Appendix A.

# 3 Steady states and their stability

In this section, we find biologically feasible steady states of the system (1) and discover the dynamical nature of the same. To obtain the equilibria, we need to solve the system of equations:

$$
\begin{align}
A_m r(A_m) + A_m f(I) - h(A_m) &= 0, \nonumber \\
h(A_m) - s(M)A_\beta - \delta_\beta A_\beta &= 0, \nonumber \\
\gamma_p A_\beta - \delta_p \tau_p - g(\tau_p) &= 0, \nonumber \\
g(\tau_p) - \delta_N N_F &= 0, \nonumber \\
\gamma_M - \delta_M M + (f_M(A_\beta) + g_M(N_F))(\tilde{M} - M)M &= 0, \nonumber \\
h_I(A_\beta)M - \delta_I I &= 0. \tag{3}
\end{align}
$$

## 3.1 Existence and stability of disease-free equilibrium

A straight forward calculation leads to the existence of a unique disease-free equilibrium $(0, 0, 0, 0, \frac{\gamma_M}{\delta_M}, 0)$ of the system (1).
For examining the stability nature of the system (1), we employ the Jacobian matrix

$$ 
\begin{pmatrix}
L_1 & 0 & 0 & 0 & 0 & A_m f'(I) \\
h'(A_m) & -s(M) - \delta_\beta & 0 & 0 & -s'(M)A_\beta & 0 \\
0 & \gamma_p & -\delta_p - g'(\tau_p) & 0 & 0 & 0 \\
0 & 0 & g'(\tau_p) & -\delta_N & 0 & 0 \\
0 & f'_M(A_\beta)(\tilde{M} - M)M & 0 & g'_M(N_F)(\tilde{M} - M)M & L_2 & 0 \\
0 & h'_I(A_\beta)M & 0 & 0 & h_I(A_\beta) & -\delta_I
\end{pmatrix},
$$

where $L_1 = A_m r'(A_m) + r(A_m) + f(I) - h'(A_m)$ and $L_2 = -\delta_M + (f_M(A_\beta) + g_M(N_F))(\tilde{M} - 2M)$.

The eigenvalues are $r_m, -s(\frac{\gamma_M}{\delta_M}) - \delta_\beta, -\delta_p, -\delta_N, -\delta_M$ and $-\delta_I$. Since, $r_m > 0$, this implies that the disease-free equilibrium is unstable.

## 3.2 Existence of disease equilibrium

**Theorem 1.** A disease equilibrium ($A_m^*, A_\beta^*, \tau_p^*, N_F^*, M^*, I^*$) exists in the system (1) if the following condition holds

$$ \tilde{M}\delta_M B(K) - 2\gamma_M B(K) + \delta_M \sqrt{(B(K)\tilde{M} - \delta_M)^2 + 4\gamma_M B(K)} > \delta_M^2, $$

for $\tau_p < K(> 0)$.

*Proof.* The proof is given in Appendix B. □

## 4 The mathematical model with its functional forms

In the earlier sections, we have examined the model in its general expression of nonlinearities. Further analysis of this general form is quite challenging. Hence, to investigate deeper the dynamical behavior of the model (1), we conduct direct calculations and explore numerical solutions by replacing the nonlinear functions in system (1) with explicit expressions. These specific choices are based on their biological relevance (Section 2) and are aligned with the properties outlined in Subsection 2.4. The forms of the functions in the mathematical model (1) can be taken as

$$ r(A_m) = \frac{r_m}{k_m} (k_m - A_m), \quad f(I) = \frac{\mu_I I}{\alpha_I + I}, \quad h(A_m) = \mu_m A_m^3, \quad s(M) = \frac{\mu_M M}{1 + \alpha_M M}, $$

$$ g(\tau_p) = \mu_p \tau_p^2, \quad f_M(A_\beta) = C_M A_\beta, \quad g_M(N_F) = C_N N_F, \quad h_I(A_\beta) = \frac{\mu_\beta A_\beta}{\alpha_\beta + A_\beta}. $$

Once we substitute the functional forms into the mathematical model (1), the resulting model
is as follows:

$$
\left\{
\begin{aligned}
\frac{dA_m}{dt} &= \frac{r_m A_m}{k_m} (k_m - A_m) + \frac{\mu_I A_m I}{\alpha_I + I} - \mu_m A_m^3, \\
\frac{dA_\beta}{dt} &= \mu_m A_m^3 - \frac{\mu_M A_\beta M}{\alpha_M M + 1} - \delta_\beta A_\beta, \\
\frac{d\tau_p}{dt} &= \gamma_p A_\beta - \delta_p \tau_p - \mu_p \tau_p^2, \\
\frac{dN_F}{dt} &= \mu_p \tau_p^2 - \delta_N N_F, \\
\frac{dM}{dt} &= \gamma_M - \delta_M M + (C_M A_\beta + C_N N_F)(\tilde{M} - M)M, \\
\frac{dI}{dt} &= \frac{\mu_\beta A_\beta M}{\alpha_\beta + A_\beta} - \delta_I I,
\end{aligned}
\right.
\qquad (4)
$$

with the initial condition $A_m(0) \ge 0, A_\beta(0) \ge 0, \tau_p(0) \ge 0, N_F(0) \ge 0, M(0) \ge 0$ and $I(0) \ge 0$.

All parameters of the model (4) are positive and their corresponding descriptions are given in Table 1.

<table><caption>Table 1: Descriptions of parameters of the model (4)</caption><thead><tr><th>Parameters</th><th>Descriptions</th><th>Baseline values</th><th>References</th><th>Ranges</th></tr></thead><tbody><tr><td>r<sub>m</sub></td><td>Growth rate of Aβ monomers</td><td>7.49326 h<sup>-1</sup></td><td>fitted</td><td>[5.620, 9.3665]</td></tr><tr><td>k<sub>m</sub></td><td>Carrying capacity of Aβ monomers</td><td>125.721 μM</td><td>fitted</td><td>[94.281, 157.161]</td></tr><tr><td>μ<sub>I</sub></td><td>Activation production rate of Aβ monomer by TNF-α</td><td>6.18925 h<sup>-1</sup></td><td>fitted</td><td>[4.642, 7.736]</td></tr><tr><td>α<sub>I</sub></td><td>Saturation coefficient of TNF-α</td><td>2.51104 μM</td><td>fitted</td><td>[1.883, 3.139]</td></tr><tr><td>μ<sub>m</sub></td><td>Conversion rate from Aβ monomers to Aβ oligomers</td><td>0.000618 h<sup>-1</sup>μM<sup>-2</sup></td><td>fitted</td><td>[0.000464, 0.000773]</td></tr><tr><td>μ<sub>M</sub></td><td>Clearance rate of Aβ oligomers due to microglia</td><td>5.01579 cells.mm<sup>-3</sup>h<sup>-1</sup></td><td>fitted</td><td>[3.762, 6.269]</td></tr><tr><td>α<sub>M</sub></td><td>Saturation coefficient of microglia</td><td>8.50706 × 10<sup>-6</sup> cells.mm<sup>-3</sup></td><td>fitted</td><td>[6.380 × 10<sup>-6</sup>, 1.0635 × 10<sup>-6</sup>]</td></tr><tr><td>δ<sub>β</sub></td><td>Degeneration rate of Aβ oligomers</td><td>0.3465 h<sup>-1</sup></td><td>estimated</td><td>[0.2599, 0.4330]</td></tr><tr><td>γ<sub>p</sub></td><td>Phosphorylation rate of τ by Aβ oligomers</td><td>0.00020373 μMh<sup>-1</sup>nM<sup>-1</sup></td><td>fitted</td><td>[0.0001528, 0.0002546]</td></tr><tr><td>δ<sub>p</sub></td><td>Degeneration rate of τ protein</td><td>0.0012 h<sup>-1</sup></td><td>estimated</td><td>[0.0009, 0.0015]</td></tr><tr><td>μ<sub>p</sub></td><td>Conversion rate from τ monomers to NFTs</td><td>0.170351 h<sup>-1</sup>μM<sup>-1</sup></td><td>fitted</td><td>[0.1278, 0.2129]</td></tr><tr><td>δ<sub>N</sub></td><td>Degeneration rate of NFTs</td><td>6.49667 h<sup>-1</sup></td><td>fitted</td><td>[4.8725, 8.1208]</td></tr><tr><td>γ<sub>M</sub></td><td>Proliferation rate of microglia</td><td>0.2532 cells.mm<sup>-3</sup>h<sup>-1</sup></td><td>estimated</td><td>[0.1899, 0.3165]</td></tr><tr><td>δ<sub>M</sub></td><td>Degeneration rate of microglia</td><td>0.0000634 h<sup>-1</sup></td><td>estimated</td><td>[0.0000476, 0.0000792]</td></tr><tr><td>C<sub>M</sub></td><td>Activation rate of microglia by Aβ oligomers</td><td>0.000095 h<sup>-1</sup></td><td>[2]</td><td>[0.0000713, 0.0001188]</td></tr><tr><td>C<sub>N</sub></td><td>Activation rate of microglia by NFTs</td><td>0.00833 h<sup>-1</sup></td><td>[2]</td><td>[0.00625, 0.01042]</td></tr><tr><td>M̃</td><td>Maximum capacity of microglia</td><td>1 cells.mm<sup>-3</sup></td><td>[8]</td><td>[0.75, 1.25]</td></tr><tr><td>μ<sub>β</sub></td><td>Production rate of TNF-α by Aβ oligomers</td><td>0.0012 μM mm<sup>3</sup>h<sup>-1</sup> cells<sup>-1</sup></td><td>[2]</td><td>[0.0009, 0.0015]</td></tr><tr><td>α<sub>β</sub></td><td>Saturation coefficient of Aβ oligomers</td><td>2.68854 nM</td><td>fitted</td><td>[2.016, 3.361]</td></tr><tr><td>δ<sub>I</sub></td><td>Degeneration rate of TNF-α</td><td>9.041 h<sup>-1</sup></td><td>estimated</td><td>[6.781, 11.302]</td></tr></tbody></table>

We have already determined the steady states for the generalized mathematical model (4) which is given in Section 3. We get two equilibria: disease-free equilibrium and disease equilibrium. As mentioned previously, the disease-free equilibrium $(0, 0, 0, 0, \frac{\gamma_M}{\delta_M}, 0)$ is unstable. The disease equilibrium can be represented as

$$A_{\beta} = \frac{\delta_{p}\tau_{p} + \mu_{p}\tau_{p}^{2}}{\gamma_{p}}, \quad N_{F} = \frac{\mu_{p}\tau_{p}^{2}}{\delta_{N}},$$

$$M = \frac{\tilde{M}\alpha - \delta_M + \sqrt{(\tilde{M}\alpha - \delta_M)^2 + 4\alpha\gamma_M}}{2\alpha}, \quad \text{where } \alpha = C_M A_\beta + C_N N_F,$$

and

$$I = \frac{\mu_{\beta} A_{\beta} M}{(\alpha_{\beta} + A_{\beta}) \delta_{I}}.$$
We clearly observe that $A_β$, $N_F$, $M$ and $I$ can be defined in terms of $A_m$ and $τ_p$, where $A_m$ and $τ_p$ are the solution of the below system:

$$ \tilde{f}(A_m, \tau_p) = \frac{r_m}{k_m} (k_m - A_m) + \frac{\mu_I I}{\alpha_I + I} - \mu_m A_m^2 = 0, \qquad (5) $$

$$ \tilde{g}(A_m, \tau_p) = \mu_m A_m^3 - \frac{\mu_M A_β M}{\alpha_M M + 1} - \delta_β A_β = 0. \qquad (6) $$

Finding an explicit analytical solution for this system is somewhat challenging due to its complexity. To ensure the existence of equilibria, we have given a graphical representation in numerical section (Figure 3).

# 5 Numerical simulation

Parameters and their corresponding values play a significant role in the field of mathematical biology. By having the intention to get biological equivalent parameter values, in this section, we estimate them using two methods: (i) directly from the literature and (ii) through data fitting based on published data. Using these parameter values, we validate analytical results further. Combining the UA with Partial Rank Correlation Coefficient (PRCC), we check the sensitivity of our model output with respect to the model parameters. We apply the standard LHS-PRCC scheme and draw the scatter plots with correlation index, PRCC, and also provide the respected p-values.

## 5.1 Literature-derived parameter estimation

In this subsection, we gather data from various scientific or clinical studies and papers to determine the values of parameters involved in the mathematical model. The parameters of the model are estimated as follows:

**Estimation of $\delta_M$**: The half life of microglia is approximately 15 months [41, 42], which is equal to $1.0957 \times 10^4$ hours ($h$). Hence the degeneration rate of microglia can be calculated as

$$ \delta_M = \frac{\ln 2}{1.0957 \times 10^4} h^{-1} = 0.633 \times 10^{-4} h^{-1}. $$

**Estimation of $\gamma_M$**: The density of microglia is 3000-9000 cells.mm$^{-3}$ [43]. For healthy brain, the steady state of microglia from the fifth equation of the model (4) is $\gamma_M - \delta_M M = 0$. Let us assume the density as 4000 cells.mm$^{-3}$. Hence, the growth rate of microglia can be calculated as $\gamma_M = \delta_M M = 4000 \times 0.633 \times 10^{-4} cells.mm^{-3}.h^{-1} = 0.2532 cells.mm^{-3}.h^{-1}$.

**Estimation of $\delta_I$**: The half life of TNF-$\alpha$ is 4.6 min [44, 45]. The degeneration rate of TNF-$\alpha$ is

$$ \frac{\ln 2}{0.076667} = 9.04104 h^{-1}. $$

**Estimation of $\delta_β$**: The majority of $A_β$ oligomers has half life 1-5 $h$ [36]. Let us take, the half life is 2 $h$, hence the degeneration rate is calculated as

$$ \frac{\ln 2}{2} = 0.3465 h^{-1}. $$
**Estimation of $\delta_p$**: The half life of $\tau$ protein in the brain is 552 h [46]. The degeneration rate can be determined by

$$ \frac{\ln 2}{552} = 0.0012 h^{-1}. $$

## 5.2 Data-driven parameter estimation

This subsection focuses on estimating the rest of the parameters through data fitting. The key factors in AD are the misfolding and the aggregation of proteins into oligomers and fibrillar structures. Here, we concentrate on fitting the data for Aβ oligomers concentration and $\tau$ protein concentration across time. We try to find good estimate of the parameters $k_m, r_m, \mu_I, \alpha_I, \mu_m, \mu_M, \alpha_M, \gamma_p, \mu_p, \delta_N$ and $\alpha_\beta$ of the mathematical model (4).

We fit the variables to the experimental data using a nonlinear least-squares optimization approach. This approach was selected due to the inherent nonlinearity in our model. The fitting process was conducted using NonlinearModelFit algorithm in Mathematica software. This function optimizes the model parameters by minimizing the sum of squared residuals between the observed data and the model's predictions. Although the number of observed data points was fewer compared to the complexity of the model and parameters, we tried to achieve the best possible fit. In such a scenario, the fitting process becomes a little challenging.

We have adapted the Aβ oligomers concentration data and $\tau$ protein concentration data from the papers by Dear et al. [36, 47]. The data points displayed in Figure 2(a) for Aβ oligomer concentration are sourced from Fig. 6(a) (page 6242) in [47], while the data points shown in Figure 2(b) for $\tau$ protein concentration are obtained from Fig. 8(a) (page 6244) in [47]. The data used in the study was obtained from *in vitro* experiments, conducted in controlled laboratory conditions using experimental techniques such as single-molecule fluorescence resonance energy transfer (FRET) and induced aggregation reactions. In Figure 2, the red data points indicate the observed concentrations of Aβ oligomers and $\tau$ protein over time, and the curves represent the model's fit to these experimental results. The estimated parameter values are given in Table 1 columnized as baseline values.
<figure><img src="image_4.png" /><figcaption>Figure 2: Data points are used to estimate the parameters $k_m, r_m, \mu_I, \alpha_I, \mu_m, \mu_M, \alpha_M, \gamma_p, \mu_p, \delta_N$ and $\alpha_\beta$ in the mathematical model (4). The figure (a) describes the data and corresponding best fit curve for the population of Aβ oligomers and (b) describes the data and corresponding best fit curve for the population of $\tau$ protein. Source of data: [36, 47]</figcaption></figure>

## 5.3 Model validation and primary analysis

Figure (3) shows the existence of disease-free equilibrium and the disease equilibrium for the parameter values given in Table 1. In Figure 3, two interaction points are shown in the $(\tau_p, A_m)$ plane, where $\tau_p$ and $A_m$ denote the $\tau$ protein and Aβ monomer population, respectively. This confirms the existence of two equilibria (0, 0, 0, 0, 3993.69, 0) and (72.0064, 0.793138, 0.0274771, 0.0000198, 57.986, 0.00175) of the system (4). We have taken a set of initial points from where all the trajectories are starting as $A_m(0) = 70, A_\beta(0) = 0.79, \tau_p(0) = 0.02, N_F(0) = 0.00001, M(0) = 57, I(0) = 0.001$.

<figure><img src="image_5.png" /><figcaption>Figure 3: Graphical representation of the existence of equilibria of the system (4).</figcaption></figure>

In Figure 4, we have plotted the temporal dynamics of the population of Aβ protein oligomers
and NFTs, taking into account the influence of the immune response (microglia population) on them. For plotting the figure, we have taken two cases: higher activation rates with lower proliferation rate of microglia and lower activation rates with higher proliferation rate of microglia. In the plots, the red curves represent the first case ($C_M = 0.0001188$, $C_F = 0.01042$, $\gamma_M = 0.1899$), whereas the black curves define the second case ($C_M = 0.0000713$, $C_N = 0.00625$, $\gamma_M = 0.3165$). The other parameters values are taken from Table 1. When the proliferation rate of microglia is high and the activation rates of microglia are low, the both of Aβ oligomers and NFTs population decreases.

<figure><img src="image_6.png" /><figcaption>Figure 4: Temporal dynamics of Aβ oligomers and NFTs at different activation and proliferation rate of microglia.</figcaption></figure>

## 5.4 Uncertainty and sensitivity analysis

To capture the most sensitive parameters of this complex mathematical model (4), we find the Partial Rank Correlation Coefficient (PRCC). To get a clearer understanding of the parameters sensitivity, we examine scatter plots depicting the relationships between model variables and individual parameters. Remind that PRCC is a statistical technique used to measure the sensitivity of model outputs to changes in input [48, 49]. Generally, PRCC values range from -1 to 1, with values closer to -1 or 1 indicating stronger negative or positive correlations, respectively. A scatter plot is a dot graph used to display the relationship between two variables. Patterns in the distribution of points are useful for identifying correlations between the variables [50, 51].

We have taken the baseline parameters values from Table 1 and the range of the parameters are considered in 25% intervals of the baseline values. For generating PRCC and scatter plots, we used an inspired MATLAB code provided by Marino et al. [52].

We perform PRCC analysis at four different time points to investigate which input variables (model parameters) have a great influence on the output variables (model populations). The four different time points that have been taken, are $t = 100 h$, $200 h$, $300 h$ and $400 h$, and for doing this sensitivity analysis the sample size is set to 2000. In our model, we are mainly interested in investigating the parameter influence on two important populations Aβ oligomers ($A_β$) and NFTs ($N_F$).
<figure><img src="image_7.png" /><figcaption>Figure 5: PRCC plot for the Aβ oligomers at four different times $t= 100, 200, 300, 400$ h.</figcaption></figure>

Figure 5 shows the PRCC plot for the protein population of Aβ oligomers at different time points $t= 100$ h, $200$ h, $300$ h and $400$ h. These four time points have been picturised in four slots 1, 2, 3 and 4, respectively. A color chart is attached which shows the color indicter and corresponding parameters. In Figure 5, PRCC results show how the parameters $r_m, k_m, \mu_m, \mu_M, \gamma_M, C_M$ and $\tilde{M}$ are consistently significant and most important parameters for the population of Aβ oligomers. It suggests that the dynamics of the population of Aβ oligomers is strongly affected by the variation of parameters $r_m, k_m, \mu_m, \mu_M, \gamma_M, C_M$ and $\tilde{M}$. PRCC scatter plots are also produced for the Aβ oligomers with respect to these parameters.

<figure><img src="image_8.png" /><figcaption>Figure 6: Scatter plot showing the relationship between Aβ oligomer population and the parameters: PRCC and the corresponding p-values of the parameters (a) $r_m$, (b) $k_m$, (c) $\mu_m$, (d) $\mu_M$, (e) $\gamma_M$, (f) $C_M$, (g) $\tilde{M}$ for the population Aβ oligomers at time $t = 400$ h.</figcaption></figure>

In Figure 6, the axes indicate the residuals of the linear regression between the rank-transformed
values of the Aβ population against the rank-transformed values of each of the parameters. Here, we consider that the parameters are statistically significant with the significance value 0.05. The growth rate $r_m$ and carrying capacity $k_m$ of Aβ monomer have a strong positive correlation (PRCC 0.89934 and 0.94453) with the Aβ population (Figure 6(a) and (b)). This suggests that as the rates at which Aβ monomers produce increases, the number of Aβ oligomers also increases significantly. We can see in Figure 6(c), the parameter $\mu_m$, conversion rate from Aβ monomers to oligomers, has a moderate positive correlation with Aβ population as the PRCC value is 0.24667 and p-value 7.6169e-29 < 0.05. The parameters with large absolute PRCC values, $PRCC > |0.5|$, with corresponding p-values $\le 0.05$, are considered to be the most important parameters [54, 55]. The Aβ population has the strong negative correlation with the parameters $\mu_M$ (the clearance rate of Aβ oligomers) and $\gamma_M$ (proliferation rate of microglia), with PRCC values -0.91705 and -0.81704, respectively (Figure 6(d) and (e)). This indicates that a decrease in the parameters $\mu_M$ and $\gamma_M$, governed by the microglia population, is associated with an increase in the Aβ oligomer population. The parameter $C_M$, activation rate of microglia by Aβ oligomers, has a positive correlation with the Aβ population with the PRCC value 0.64061 (Figure 6(f)). While the correlation is positive, it is not as strong as in some other cases, indicating that $C_M$ plays a significant but not dominant role. The parameter $\tilde{M}$ has a weak negative correlation with the Aβ oligomers population as the PRCC value is -0.06385 and the p-value is 0.0044694 < 0.05 (Figure 6(g)). That suggests statistical significance but the strength of the correlation is relatively low. This weak correlation suggests that variations in the maximum capacity of microglia, $\tilde{M}$ have a minimal effect on the population of Aβ oligomers. Other parameters $\mu_I$, $\alpha_I$, $\alpha_M$, $\delta_\beta$, $\gamma_p$, $\delta_p$, $\mu_p$, $\delta_N$, $\delta_M$, $C_N$, $\mu_\beta$, $\alpha_\beta$ and $\delta_I$ have very minimal PRCC values and p-values are greater than 0.05. This suggests that changes in these parameters have an almost negligible influence on the population of Aβ oligomers.

<figure><img src="image_9.png" /><figcaption>Figure 7: PRCC plot for the neurofibrillary tangles at four different times $t= 100, 200, 300, 400$ h.</figcaption></figure>

Figure 7 shows the PRCC plot for the protein population of neurofibrillary tangles ($N_F$) at different time points $t= 100 h, 200 h, 300 h$ and $400 h$. The consistently significant parameters for $N_F$ population are $r_m, k_m, \mu_m, \mu_M, \gamma_p, \delta_p, \mu_p, \delta_N, \gamma_M, C_M$ and $\tilde{M}$. Hence, the variation of these parameters mainly impacts the protein population $N_F$. For a graphical visualization, PRCC scatter plot are also given for $N_F$ with relative to the significant parameters.
<figure><img src="image_10.png" /><figcaption>Figure 8: Scatter plot showing the relationship between NFTs population and the parameters: PRCC and the corresponding p-values of the parameters (a) $r_m$, (b) $k_m$, (c) $\mu_m$, (d) $\mu_M$, (e) $\gamma_p$, (f) $\delta_p$, (g) $\mu_p$, (h) $\delta_N$, (i) $\gamma_M$, (j) $C_M$, (k) $\tilde{M}$ for the population NFTs population at time $t = 400h$.</figcaption></figure>

Here, we have taken the statistically significance value at 0.05. The growth rate $r_m$ and carrying capacity $k_m$ of A$\beta$ monomer have a strong positive correlation (PRCC 0.87318 and 0.9263, respectively) with the $N_F$ population (Figure 8(a) and 8(b)). This strong correlation indicates that variations in two parameters $r_m$ and $k_m$, are crucial in affecting the NFTs population. Remind that parameters with large absolute PRCC values, PRCC > |0.5|, with corresponding p-values < 0.05, are considered to be the most important parameters [54, 55]. Besides $r_m$ and $k_m$, the parameters $\mu_m$, $\mu_M$, $\gamma_p$, $\delta_N$, $\gamma_M$, $\delta_p$, $\mu_p$, $\tilde{M}$ and $C_M$ are also significant and have impact on the NFTs population (see Figure 8). Rather than the significant parameters, there are many parameters that have very less impact on the NFTs population, which is negligible. The parameters $\mu_I$, $\alpha_I$, $\alpha_M$, $\delta_\beta$, $\delta_M$, $C_N$, $\mu_\beta$, $\alpha_\beta$ and $\delta_I$ has less PRCC value and p-value > 0.05.

# 6 Discussion and conclusion

In this paper, we have formulated a mathematical model based on the events: aggregation of proteins, activation of immune cells and initiation of inflammation. Aggregation of proteins (A$\beta$ plaque and NFTs) has been recognized as a hallmark in the progression of AD since its discovery. Having a critical relevance to AD, our model includes two key disease-related proteins along with their aggregated forms. Microglia, which have an intricate involvement in AD, are known to regulate both phagocytosis and inflammatory responses [53]. These immune cells contribute to neuroinflammation
and can either worsen or lessen disease progression depending on their activation state. Recogniz-
ing the significance of microglia and inflammation in the progression of AD, we have included both
microglia and cytokines as components of our model. We believe that our model generalizing the
existing ones, provides flexibility and is able to capture realistic biological dynamics.

After formulating the mathematical model, we have conducted an extensive review of both
biological and mathematical literature to identify the relevant properties of the functional forms.
This exploration allowed us to choose functional forms that not only represent the underlying
biological processes but also follow mathematical principles.

We have investigated the fundamental properties of our mathematical model and analyzed its
dynamical behavior by identifying the equilibrium points. Once the equilibria are determined,
stability analysis have been performed to understand the system's behavior near these points. We
get that the disease-free equilibrium is always unstable, implying that once pathological processes
begin, the system cannot recover to a completely normal state. This is consistent with the largely
irreversible character of AD, while early-stage protective mechanisms may reverse the progression.

Due to the complexity of the model, determining the exact form of the disease equilibrium
was challenging. However, through mathematical calculations and numerical simulation, we have
demonstrated the existence of a disease equilibrium for the system.

When we formulated the generalized mathematical model, some basic properties of the functions
were assumed. Based on the biological significant assumptions, we have considered the explicit forms
of the model.

In numerical solutions, we have done the parameter estimation through two methods: from
biological literature and data fitting. To improve the model's predictive power, we have performed
the data fitting on Aβ oligomer and τ protein datasets to calibrate some of the model parameters.
Considering the estimated values as baseline values of the parameters, we have identified key pa-
rameters that significantly influence the model's output through the SA and UA. This allows us to
choose which factors we should focus on for more accurate predictions and a deeper understanding
of the system's dynamics. To understand the scenario globally, we have taken the help of PRCC
determination and scatter plots. In PRCC plot, we compared the situations at four different time
points 100 h, 200 h, 300 h and 400 h.

The production or growth of Aβ monomer, that is, the growth rate r_m and the carrying capacity k_m has a positive impact on both the Aβ oligomer and NFTs population. Since in both the cases, we have seen that the PRCC values are > 0.8 implies a strong positive correlation between the parameters and model populations. Also observing the scatter plots, we conclude that increase in growth factors of Aβ monomer, leads to increase in both Aβ oligomer and NFTs population.

In the Figure 8(e), we note that the phosphorylation rate $\gamma_p$ has a strong positive correlation
with NFTs population. As the rate of phosphorylation ($\gamma_p$) increases, there is a corresponding rise in
the NFTs population. This fact aligns with the well-established biological fact that NFTs formation
is heavily influenced by the phosphorylation of $\tau$ protein. The increased phosphorylation leads to
more $\tau$ proteins becoming dysfunctional and forming NFTs, which helps in AD progression. Hence,
the figure highlights the critical role of phosphorylation in NFT formation and AD advancement.

The two parameters $\mu_M$ and $\gamma_M$, related to microglia, have a strong negative correlation with both the Aβ oligomers and NFTs. From this, we can conclude that as the saturation and proliferation of microglia increase, the levels of Aβ oligomers and NFTs tend to decrease. This suggests that microglial activity may play a significant role in reducing the accumulation of these pathological biomarkers.

In this work, the application of SA related to UA is employed to identify the parameters with
the highest sensitivity regarding the progression of AD. PRCC helps researchers understand which
parameters are most critical in their models, guiding further research, data collection efforts and model refinement.

**Limitations of this work:** the model faces limitations, such as the exclusion of nucleation, elongation and fragmentation processes, which may limit its representation of Aβ aggregation. Another limitation is that our model does not explicitly include dissociation of monomers from oligomers. Additionally, only microglia has been considered as an immune system in our model, might not fully capture the complexity of neuroinflammation in AD. The absence of spatial dynamics in the model may limit its relevance to the progression of diseases in specific regions. We have not gone through the complete τ propagation mechanisms.

**Future directions:** for future directions, several extensions can be explored to enhance the current model of AD. Spatial dynamics can be incorporated into the model. Additionally, including nucleation, elongation and fragmentation processes could provide a more comprehensive view of Aβ aggregation and its impact on disease outcomes. To finding treatment strategies, optimal control theory can be used.

# Appendix A: Positivity and boundedness of the model

## Proof of statement 1

*Proof.* Let us consider

$$ 
\begin{aligned}
f_1(0, A_m, A_\beta, \tau_p, N_F, M, I) &= A_m r(A_m) + A_m f(I) - h(A_m), \\
f_2(0, A_m, A_\beta, \tau_p, N_F, M, I) &= h(A_m) - s(M)A_\beta - \delta_\beta A_\beta, \\
f_3(0, A_m, A_\beta, \tau_p, N_F, M, I) &= \gamma_p A_\beta - \delta_p \tau_p - g(\tau_p), \\
f_4(0, A_m, A_\beta, \tau_p, N_F, M, I) &= g(\tau_p) - \delta_N N_F, \\
f_5(0, A_m, A_\beta, \tau_p, N_F, M, I) &= \gamma_M - \delta_M M + (f_M(A_\beta) + g_M(N_F))(\tilde{M} - M)M, \\
f_6(0, A_m, A_\beta, \tau_p, N_F, M, I) &= h_I(A_\beta)M - \delta_I I.
\end{aligned} 
\quad (7) $$

Now, for $A_m = 0$, $f_1 = 0$, for $A_\beta = 0$, $f_2 = h(A_m) \ge 0$, for $\tau_p = 0$, $f_3 = \gamma_p A_\beta \ge 0$, for $N_F = 0$, $f_4 = g(\tau_p) \ge 0$, for $M = 0$, $f_5 = \gamma_M \ge 0$, for $I = 0$, $f_6 = h_I(A_\beta)M \ge 0$.

This proves that $\mathbb{R}_+^6$ is an invariant set with respect to the dynamical system (1). Hence for, $A_m(0) \ge 0$, $A_\beta(0) \ge 0$, $\tau_p(0) \ge 0$, $N_F(0) \ge 0$, $M(0) \ge 0$, $I(0) \ge 0$, then $A_m(t) \ge 0$, $A_\beta(t) \ge 0$, $\tau_p(t) \ge 0$, $N_F(t) \ge 0$, $M(t) \ge 0$, $I(t) \ge 0$.

□

## Proof of statement 2

*Proof.* The solutions of the system have the following properties:

* The maximum value of the term $f(I)$ is $\mu_I$. The first equation of model (1) can be written as $\frac{dA_m}{dt} \le A_m r(A_m) + A_m f(I)$. In absence of inflammation, the maximum value of $A_m$ will be $k_m$, the carrying capacity for Aβ monomer population. But, here $f(I)$ positively influences the production of Aβ monomer, it could allow the population to reach a higher maximum size than the carrying capacity. Moreover, $\frac{dA_m}{dt} \le A_m(r(A_m) + \mu_I)$. Hence the maximum value of
$A_m$ could be $(A_m)_{\max} = b k_m$, where $b > 0$ is a constant. We can determine the constant $b$ depending on the specific form of the expressions.

* For the variable $A_\beta$, $\frac{dA_\beta}{dt} \le h(bk_m) - \delta_\beta A_\beta$ which implies the maximum value of $A_\beta$ is $(A_\beta)_{\max} = \max\left\{\frac{h(bk_m)}{\delta_\beta}, A_\beta(0)\right\}$.

* In the third equation, $\frac{d\tau_p}{dt} \le \gamma_p(A_\beta)_{\max} - \delta_p \tau_p$ implies $\tau_p \le (\tau_p)_{\max} = \max\left\{\frac{\gamma_p(A_\beta)_{\max}}{\delta_p}, \tau_p(0)\right\} = K(say)$.

* For the variable $N_F$, $\frac{dN_F}{dt} \le g((\tau_p)_{\max}) - \delta_N N_F$ implies $N_F \le (N_F)_{\max} = \max\left\{\frac{g((\tau_p)_{\max})}{\delta_N}, N_F(0)\right\}$.

* In the fifth equation, $\frac{dM}{dt} \le \gamma_M - \delta_M M + ((A_\beta)_{\max} + (N_F)_{\max})\frac{\tilde{M}^2}{4}$ implies $M \le M_{\max} = \max\left\{\frac{((A_\beta)_{\max} + (N_F)_{\max})\tilde{M}^2 + 4\gamma_M}{4\delta_M}, M(0)\right\}$.

* The maximum value of the term $h_I(A_\beta)$ is $\mu_I$. Hence, $\frac{dI}{dt} \le \mu_I M_{\max} - \delta_I I$ implies $I \le I_{\max} = \max\left\{\frac{\mu_I M_{\max}}{\delta_I}, I(0)\right\}$.

□

# 7 Appendix B: Existence of disease equilibrium

## Proof of Theorem 1

The disease equilibrium is given by $(A_m^*, A_\beta^*, \tau_p^*, N_F^*, M^*, I^*)$ where

$$A_\beta^* = \frac{\delta_p \tau_p^* + g(\tau_p^*)}{\gamma_p}, \quad N_F^* = \frac{g(\tau_p^*)}{\delta_N}, \quad I^* = \frac{h_I(A_\beta^*) M^*}{\delta_I}.$$

To find the value of $M^*$, we solve the fifth equation of the system (3)

$$B(\tau_p^*) M^{*2} - (B(\tau_p^*) \tilde{M} - \delta_M) M^* - \gamma_M = 0, \quad (8)$$

where

$$B(\tau_p^*) = f_M(A_\beta^*) + g_M(N_F^*) = f_M\left(\frac{\delta_p \tau_p^* + g(\tau_p^*)}{\gamma_p}\right) + g_M\left(\frac{g(\tau_p^*)}{\delta_N}\right).$$

It is clear that equation (8) has always one positive root regardless any conditions. The solution is

$$M^* = \frac{(B(\tau_p^*) \tilde{M} - \delta_M) + \sqrt{(B(\tau_p^*) \tilde{M} - \delta_M)^2 + 4B(\tau_p^*) \gamma_M}}{2B(\tau_p^*)}.$$

Now to find the disease steady states of the variables $A_m$ and $\tau_p$, we need to solve the first two equations

$$r(A_m^*) + f(I^*) - h(A_m^{*\frac{n-1}{n}}) = 0, \quad (9)$$

$$h(A_m^*) - s(M^*) A_\beta^* - \delta_\beta A_\beta^* = 0. \quad (10)$$
Putting the values of $A_{\beta}^{*}$, $M^{*}$ and $I^{*}$, the equations (9) and (10) read as

$$r(A_{m}^{*}) + f\left(\frac{h_{I}\left(\frac{\delta_{p}\tau_{p}^{*}+g(\tau_{p}^{*})}{\gamma_{p}}\right)\left((B(\tau_{p}^{*})\tilde{M}-\delta_{M})+\sqrt{(B(\tau_{p}^{*})\tilde{M}-\delta_{M})^{2}+4B(\tau_{p}^{*})\gamma_{M}}\right)}{2B(\tau_{p}^{*})\delta_{I}}\right) - h(A_{m}^{*\frac{n-1}{n}}) = 0, (11)$$

$$h(A_m^*) - \left( s \left( \frac{(B(\tau_p^*)\tilde{M} - \delta_M) + \sqrt{(B(\tau_p^*)\tilde{M} - \delta_M)^2 + 4B(\tau_p^*)\gamma_M}}{2B(\tau_p^*)} \right) + \delta_\beta \right) \frac{\delta_p \tau_p^* + g(\tau_p^*)}{\gamma_p} = 0, (12)$$

Let us assume an auxiliary function $F(\tau_p)$ of the equation (12)

$$F(\tau_p) = h(A_m) - \left( s \left( \frac{(B(\tau_p)\tilde{M} - \delta_M) + \sqrt{(B(\tau_p)\tilde{M} - \delta_M)^2 + 4B(\tau_p)\gamma_M}}{2B(\tau_p)} \right) + \delta_\beta \right) \frac{\delta_p \tau_p + g(\tau_p)}{\gamma_p} \quad (13)$$

for $A_m > 0$. For every $A_m > 0$, $F(0) = h(A_m) > 0$ and $\lim_{\tau_p \to \infty} F(\tau_p) = -\infty$. Hence, if $F(\tau_p)$ is a decreasing function, that is, $F'(\tau_p) < 0$ then for each $A_m$, there is exactly one positive steady state for $\tau_p$. Now differentiating (13) with respect to $\tau_p$, we have

$$F'(\tau_p) = -\frac{1}{\gamma_p} \left( (g'(\tau_p) + \delta_p) \left( s \left( \frac{\sqrt{(B(\tau_p)\tilde{M} - \delta_M)^2 + 4\gamma_M B(\tau_p) + B(\tau_p)\tilde{M} - \delta_M}}{2B(\tau_p)} \right) + \delta_\beta \right) \right) \\ -\frac{1}{\gamma_p} \frac{\left( \delta_M \sqrt{(B(\tau_p)\tilde{M} - \delta_M)^2 + 4\gamma_M B(\tau_p) + \tilde{M}\delta_M B(\tau_p) - 2\gamma_M B(\tau_p) - \delta_M^2} \right)}{2B(\tau_p)^2 \sqrt{(B(\tau_p)\tilde{M} - \delta_M)^2 + 4\gamma_M B(\tau_p)}} \\ \times B'(\tau_p) (g(\tau_p) + \delta_p \tau_p) s' \left( \frac{\sqrt{(B(\tau_p)\tilde{M} - \delta_M)^2 + 4\gamma_M B(\tau_p) + B(\tau_p)\tilde{M} - \delta_M}}{2B(\tau_p)} \right).$$

From our assumptions, we have the properties that $g$, $s$, $f_M$ and $g_M$ are positive and $g'$, $s'$, $f'_M$ and $g'_M$ are increasing, decreasing, nondecreasing and nondecreasing, respectively. It follows that $F'(\tau_p)$ is decreasing when the condition holds:

$$\tilde{M}\delta_{M}B(\tau_{p}) - 2\gamma_{M}B(\tau_{p}) + \delta_{M}\sqrt{(B(\tau_{p})\tilde{M} - \delta_{M})^{2} + 4\gamma_{M}B(\tau_{p})} > \delta_{M}^{2}. \quad (14)$$

Since $F'(\tau_p)$ is decreasing and negative, then $F(\tau_p)$ is decreasing for any $\tau_p > 0$.

We have $\tau_p < K(> 0)$, hence, the condition (14) can be written as

$$\tilde{M}\delta_{M}B(K) - 2\gamma_{M}B(K) + \delta_{M}\sqrt{(B(K)\tilde{M} - \delta_{M})^{2} + 4\gamma_{M}B(K)} > \delta_{M}^{2}. \quad (15)$$
The condition 15 ensures the existence of disease equilibrium within the system. This means it determines whether a stable state can be achieved where all the population persists.

**ETHICAL STATEMENT:**

**Funding:** This study of Mitali Maji is supported by University Grants Commission (UGC) [NTA Ref.: 211610045129], Govt. of India. Subhas Khajanchi acknowledges the financial support from the Department of Science and Technology (DST), Govt. of India, under the Scheme "Fund for Improvement of S&T Infrastructure (FIST)" [File No. SR/FST/MS-I/2019/41].

**Conflict of Interest:** The authors declares that they have no conflict of interest.

**Informed consent:** Not Applicable.

**Author contribution:** All the authors equally contributed to this work.

**Data Availability Statement:** All the used data are included in the manuscript.

**References**

[1] The top 10 causes of death. World Health Organization. https://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death.

[2] Hao W, Friedman A. Mathematical model on Alzheimer's disease. BMC System Biology. 10(1) (2016) 1-18.

[3] Pal S, Melnik R. Nonlocal models in the analysis of brain neurodegenerative protein dynamics with application to Alzheimer's disease. Scientific Report. 12(1) (2022) 7328.

[4] Thompson TB, Chaggar P, Kuhl E, Goriely A, Alzheimer's disease neuroimaging initiative. protein-protein interactions in neurodegenerative diseases: A conspiracy theory. PLoS Computational Biology. 16(10) (2020) e1008267.

[5] Puri IK, Li L. Mathematical modeling for the pathogenesis of Alzheimer's disease. PloS One. 5(12) (2010) e15176.

[6] Craft DL, Wein LM, Selkoe DJ. A mathematical model of the impact of novel treatments on the Aβ burden in the Alzheimer's brain, CSF and plasma. Bulletin of Mathematical Biology 64(5) (2002) 1011-1031.

[7] Lindstrom MR, Chavez MB, Gross-Sable EA, Hayden EY, Teplow DB. From reaction kinetics to dementia: A simple dimer model of Alzheimer's disease etiology. PLoS Computational Biology 17(7) (2021) e1009114.

[8] Ciuperca I, Pujo-Menjouet L, Matar-Tine L, Torres N, Volpert V. A qualitative analysis of a Aβ-monomer model with inflammation processes for Alzheimer's disease. Royal Society Open Science 11 (2024) 231536.

[9] Maji M, Khajanchi S. Mathematical models on Alzheimer's Disease and its treatment: A review. Physics of Life Reviews. 52 (2025) 207-244.

[10] Uncertainty and Sensitivity Analysis. https://designbuilder.co.uk/helpv7.0/Content/UASA.htm
[11] McKay MD, Beckman RJ, Conover WJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 42(1) (2000) 55-61.

[12] Tian W. A review of sensitivity analysis methods in building energy analysis. Renewable and Sustainable Energy Reviews. 20 (2013) 411-6.

[13] Hamby DM. A comparison of sensitivity analysis techniques. Health Physics. 68(2) (1995) 195-204.

[14] Al-Ghraiybah NF, Wang J, Alkhalifa AE, Roberts AB, Raj R, Yang E, Kaddoumi A. Glial cell-mediated neuroinflammation in Alzheimer's disease. International Journal of Molecular Sciences. 23(18) (2022) 10572.

[15] Tonnies E., Trushina E. Oxidative stress, synaptic dysfunction, and Alzheimer's disease. Journal of Alzheimer's Disease. 57(4) (2017) 1105-21.

[16] Sandberg A, Berenjeno-Correa E, Rodriguez RC, Axenhus M, Weiss SS, Batenburg K, Hoozemans JJ, Tjernberg LO, Scheper W. Aβ42 oligomer-specific antibody ALZ-201 reduces the neurotoxicity of Alzheimer's disease brain extracts. Alzheimer's Research & Therapy, 14(1) (2022) 196.

[17] Song C, Zhang T, Zhang Y. Conformational essentials responsible for neurotoxicity of Aβ42 aggregates revealed by antibodies against oligomeric Aβ42. Molecules, 27(19) (2022) 6751.

[18] Chen GF, Xu TH, Yan Y, Zhou YR, Jiang Y, Melcher K, Xu HE. Amyloid beta: structure, biology and structure-based therapeutic development. Acta Pharmacologica Sinica, 38(9) (2017) 1205-1235.

[19] Goure WF, Krafft GA, Jerecic J, Hefti F. Targeting the proper amyloid-beta neuronal toxins: a path forward for Alzheimer's disease immunotherapeutics. Alzheimer's Research & Therapy, 6 (2014) 1-15.

[20] Chen CY, Tseng YH, Ward JP. A mathematical model demonstrating the role of interstitial fluid flow on the clearance and accumulation of amyloid β in the brain. Mathematical Bioscience 317 (2019) 108258.

[21] Huang YR, Liu RT. The toxicity and polymorphism of β-amyloid oligomers. International Journal of Molecular Sciences. 21(12) (2020) 4477.

[22] De Felice FG, Wu D, Lambert MP, Fernandez SJ, Velasco PT, Lacor PN, Bigio EH, Jerecic J, Acton PJ, Shughrue PJ, Chen-Dodson E. Alzheimer's disease-type neuronal tau hyperphosphorylation induced by Aβ oligomers. Neurobiology of Aging. 29(9) (2008) 1334-47.

[23] Zhang H, Wei W, Zhao M, Ma L, Jiang X, Pei H, Cao Y, Li H. Interaction between Aβ and tau in the pathogenesis of Alzheimer's disease. International Journal of Biological Sciences. 17(9) (2021) 2181.

[24] Rankin CA, Sun Q, Gamblin TC. Tau phosphorylation by GSK-3β promotes tangle-like filament morphology. Molecular Neurodegeneration. 2 (2007) 1-4.
[25] Morimoto K, Nakajima K. Role of the immune system in the development of the central nervous system. Frontiers in Neuroscience, 13 (2019) 916.

[26] Doens D, Fernandez PL. Microglia receptors and their implications in the response to amyloid β for Alzheimer's disease pathogenesis. Journal of Neuroinflammation, 11 (2014) 1-14.

[27] Salvadores N, Moreno-Gonzalez I, Gamez N, Quiroz G, Vegas-Gomez L, Escandon M, Jimenez S, Vitorica J, Gutierrez A, Soto C, Court FA. Aβ oligomers trigger necroptosis-mediated neurodegeneration via microglia activation in Alzheimer's disease. Acta Neuropathologica Communications. 10(1) (2022) 31.

[28] Ohm DT, Fought AJ, Martersteck A, Coventry C, Sridhar J, Gefen T, Weintraub S, Bigio E, Mesulam MM, Rogalski E, Geula C. Accumulation of neurofibrillary tangles and activated microglia is associated with lower neuron densities in the aphasic variant of Alzheimer's disease. Brain Pathology. 31(1) (2021) 189-204.

[29] Zhang G, Wang Z, Hu H, Zhao M, Sun L. Microglia in Alzheimer's disease: a target for therapeutic intervention. Frontiers in Cellular Neuroscience. 15 (2021) 749587.

[30] Magni P, Ruscica M, Dozio E, Rizzi E, Beretta G, Facino RM. Parthenolide inhibits the LPS-induced secretion of IL-6 and TNF-α and NF-κB nuclear translocation in BV-2 microglia. Phytotherapy Research. 26(9) (2012) 1405-9.

[31] Meraz-Rios MA, Toral-Rios D, Franco-Bocanegra D, Villeda-Hernandez J, Campos-Pena V. Inflammatory process in Alzheimer's Disease. Frontiers in Integrative Neuroscience. 7 (2013) 59.

[32] Alasmari F, Alshammari MA, Alasmari AF, Alanazi WA, Alhazzani K. Neuroinflammatory cytokines induce amyloid beta neurotoxicity through modulating amyloid precursor protein levels/metabolism. BioMed Research International. 2018(1) (2018) 3087475.

[33] Dear AJ, Thacker D, Wennmalm S, Ortigosa-Pascual L, Andrzejewska EA, Meisl G, Linse S, Knowles TP. Aβ Oligomer Dissociation Is Catalyzed by Fibril Surfaces. ACS Chemical Neuroscience. (2024).

[34] Maji M Khajanchi S. Roles of astrocytes and prions in Alzheimer's disease: insights from mathematical modeling. Journal of Biological Physics, 50(2) (2024) 149-179.

[35] Portet S, Arino J. An in vivo intermediate filament assembly model. Mathematical Biosciences & Engineering. 6(1) (2008) 117-34.

[36] Dear AJ, Michaels TC, Meisl G, Klenerman D, Wu S, Perrett S, Linse S, Dobson CM, Knowles TP. Kinetic diversity of amyloid oligomers. Proceedings of the National Academy of Sciences. 117(22) (2020) 12087-94.

[37] Townsend D, Fullwood NJ, Yates EA, Middleton DA. Aggregation kinetics and filament structure of a tau fragment are influenced by the sulfation pattern of the cofactor heparin. Biochemistry. 59(41) (2020) 4003-14.

[38] Chen CY, Tseng YH, Ward JP. A mathematical model on the propagation of tau pathology in neurodegenerative diseases. Journal of Mathematical Biology. 89(1) (2024) 4.
[39] Song Y. Positive invariance condition for continuous dynamical systems based on Nagumo Theorem. arXiv preprint arXiv:2207.05429. (2022).

[40] Nagumo M. Über die lage der integralkurven gewöhnlicher differential gleichungen. Proceedings of the Physico-Mathematical Society of Japan. 3rd Series. 124 (1942) 551–9.

[41] Fuger P, Hefendehl JK, Veeraraghavalu K, Wendeln AC, Schlosser C, Obermuller U, Wegenast-Braun BM, Neher JJ, Martus P, Kohsaka S, Thunemann M. Microglia turnover with aging and in an Alzheimer’s model via long-term in vivo single-cell imaging. Nature Neuroscience. 20(10) (2017) 1371–6.

[42] Vidal-Itriago A, Radford RA, Aramideh JA, Maurel C, Scherer NM, Don EK, Lee A, Chung RS, Graeber MB, Morsch M. Microglia morphophysiological diversity and its implications for the CNS. Frontiers in Immunology. 13 (2022) 997786.

[43] Keller D, Ero C, Markram H. Cell densities in the mouse brain: a systematic review. Frontiers in neuroanatomy. 12 (2018) 83.

[44] Simo R, Barbosa-Desongles A, Lecube A, Hernandez C, Selva DM. Potential role of tumor necrosis factor-α in down regulating sex hormone-binding globulin. Diabetes. 61(2) (2012) 372–82.

[45] Lai X, Hao W, Friedman A. TNF-α inhibitor reduces drug-resistance to anti-PD-1: A mathematical model. Plos One. 15(4) (2020) e0231499.

[46] Hier DB, Azizi S, Thimgan MS, Wunsch DC. Tau kinetics in Alzheimer’s disease. Frontiers in Aging Neuroscience. 14 (2022) 1055170.

[47] Dear AJ, Meisl, G, aric A, Michaels TC, Kjaergaard M, Linse S., Knowles TP. Identification of on-and off-pathway oligomers in amyloid fibril formation. Chemical Science, 11(24) (2020) 6236–6247

[48] Zi Z. Sensitivity analysis approaches applied to systems biology models. IET systems biology. 5(6) (2011) 336–46.

[49] Savatorova V. Exploring parameter sensitivity analysis in mathematical modeling with ordinary differential equations. CODEE Journal. 16(1) (2023) 4.

[50] Hastie T, Tibshirani R, Friedman J. The elements of statistical learning: data mining, inference, and prediction. (2017). Springer.

[51] Iooss B, Lemaitre P. A review on global sensitivity analysis methods. Uncertainty management in simulation-optimization of complex systems: algorithms and applications. (2015) 101–22.

[52] Marino S, Hogue IB, Ray CJ, Kirschner DE. A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of Theoretical Biology. 254(1) (2008) 178–96.

[53] Miao J, Ma H, Yang Y, Liao Y, Lin C, Zheng J, Yu M, Lan J. Microglia in Alzheimer’s disease: Pathogenesis, mechanisms, and therapeutic potentials. Frontiers in Aging Neuroscience. 15 (2023) 1201982.
[54] Nagaraja S, Chen L, DiPietro LA, Reifman J, Mitrophanov AY. Predictive approach identifies molecular targets and interventions to restore angiogenesis in wounds with delayed healing. Frontiers in Physiology. 10 (2019) 636.

[55] Taylor R. Interpretation of the correlation coefficient: a basic review. Journal of Diagnostic Medical Sonography, 6(1) (2019) 35-39.