# Spatial pattern analysis of an Aβ-monomer model with inflammation processes for Alzheimer's disease

Maxime Estavoyer, Nicolas Torres, Julien Blohm, Malay Banerjee, Laurent

Pujo-Menjouet

► To cite this version:

Maxime Estavoyer, Nicolas Torres, Julien Blohm, Malay Banerjee, Laurent Pujo-Menjouet. Spatial pattern analysis of an Aβ-monomer model with inflammation processes for Alzheimer's disease. Mathematical Modelling of Natural Phenomena, 2025, 10.1051/mmnp/2025022. hal-04703797v2

HAL Id: hal-04703797

https://inria.hal.science/hal-04703797v2

Submitted on 10 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Spatial pattern analysis of a Aβ-monomer model with
inflammation processes for Alzheimer’s disease

Maxime Estavoyer*

Nicolas Torres†

Julien Blohm‡

Malay Banerjee§

Laurent Pujo-Menjouet‡

Abstract

We study the emergence of spatial patterns for a system of reaction-diffusion equations, modeling the progression of Alzheimer's disease through the interaction of Aβ-monomers, oligomers, microglial cells, and interleukins with neurons. In our work, these spatial patterns stand for inert amyloid plaques, which are extracellular deposits of Aβ-proteins and a characteristic feature of this neurodegenerative disease. Using linear analysis and numerical simulations, we show the existence of spatially heterogeneous solutions and exhibit a wide variety of possible spatially-dependent solutions: time-oscillating, low-amplitude, and high-amplitude patterns. Moreover, we carry out an extensive analysis of high-amplitude patterns in the one- and two-dimensional domains. In particular, we study the stability of branches of heterogeneous steady states through bifurcation diagrams and their selection. From this numerical bifurcation analysis, we develop some conjectures concerning the influence of inflammation and microglial cells in the formation of amyloid plaques. These findings offer insights into potential anti-inflammatory treatments that might be used to mitigate the progression of Alzheimer's disease and the emergence of inert amyloid plaques.

**Keywords:** Bifurcation analysis, Mathematical modeling, Turing patterns.

# 1 Introduction

## Biological background and state of the art

Alzheimer's disease (AD) is a progressive neurodegenerative disorder that affects memory, thinking, and behavior. The pathogenesis of AD is multifactorial and involves several risk factors such as age, genetics, head injuries, vascular diseases, infections, and environment. One of the main hypotheses proposed to explain the disease's pathogenesis is the misfolding and aggregation of amyloid β-proteins in the nervous system.

Amyloid-beta (Aβ) is a peptide that forms the primary component of the amyloid plaques found in the brains of Alzheimer's disease patients (see figure 1). Aβ peptides are intrinsically disordered in their monomeric form and assemble into stable structures such as oligomers and fibrils via a nucleation-dependent pathway. Many studies point to a key role for these Aβ-aggregates in neurotoxicity, leading to progressive, irreversible neuronal damage (see, for

<footnote>*Inria, and Université de Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5208, Institut Camille Jordan, F-69603 Villeurbanne, France.</footnote>

<footnote>†Universidad de Granada, Departamento de Matemática Aplicada, Granada, Andalusia, Spain.</footnote>

<footnote>‡Université Claude Bernard Lyon 1, CNRS, Centrale Lyon, INSA Lyon, Université Jean Monnet, ICJ UMR5208, Inria, 69622 Villeurbanne, France.</footnote>

<footnote>§Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur - 208016, India.</footnote>
example, [1, 2, 3, 4, 5]). Understanding the role of Aβ aggregates and the mechanisms of their
formation and degradation could provide valuable insights into the pathogenesis of Alzheimer’s
disease and potential therapeutic strategies.

<figure><img src="image_2.png" /><figcaption>Figure 1: **A.** A picture of Aβ amyloid plaques in the brain.<fnref n="1" /> **B.** Stationary patterns of our mathematical model.</figcaption></figure>

Neurons produce Aβ-monomers that almost instantaneously start to polymerize into proto-
oligomers. In this aggregation process, proto-oligomers can polymerize or depolymerize, and
once they reach a critical size, they become stable in the form of Aβ-oligomers. These latter are
assumed to be completely stable in the sense that neither polymerization nor depolymerization
is possible for Aβ-oligomers in equilibrium [6, 7]. This mechanism involving Aβ-oligomers is
known as the amyloid cascade hypothesis, and there is a consensus that it is a key factor in AD
progression.

The role of microglia, the primary immune cells in the brain, in the progression of Alzheimer’s
disease is also significant. They induce an inflammatory reaction through a chemical cascade in
microglial cells, releasing interleukins [8, 9, 10, 11]. These interleukins then activate an increase
in Aβ-monomer production from neurons. However, if the concentration of Aβ-oligomers is
high enough, a stress reaction called the unfolded protein response (UPR) [4] is triggered, which
leads to a decrease in Aβ-monomer production. At the same time, the remaining oligomers
diffuse into the neuronal environment. In this context, two opposing mechanisms of stimulation
and inhibition determine whether AD persists. Moreover, microglial cells also play a role in
the transfer of oligomers to amyloid plaques, where these oligomers become inert elements (no
diffusion, polymerization, or de-polymerization) [12, 13, 14].

From a mathematical point of view, previous works on partial differential equations on AD
include the works of Andrade *et al.* [15] on a space-dependent model for the polymerization of
Aβ-proteins, Ciuperca *et al.* [16] which takes into account the formation of Aβ-oligomers and
fibrils through a continuous size model based on Lifshitz-Slyozov equations and Hao *et al.* [17]
on a multi-component model taking into account astrocytes, microglia cells and peripheral cell
population, besides the Aβ-oligomers. Other mathematical models related to amyloid plaque
formation, or more generally to Alzheimer's disease, can be found in [18, 19, 20, 21, 22, 23, 24,

<footnote n="1">¹Source: Copyright © 2011 Michael Bonert, MD, FRCPC (https://commons.wikimedia.org/wiki>User:Nephron, https://experts.mcmaster.ca/display/bonertm). You are free to share and adapt this image as per the CC BY-SA 3.0 (https://creativecommons.org/licenses/by-sa/3.0/legalcode).</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the bi-monomeric model (1) introduced in the article [27]</figcaption></figure>

25, 26].

## Model formulation

In this article, we study the bi-monomeric model introduced in the work of [27], which is given by the following reaction-diffusion system in a bounded domain $\Omega$,

$$ egin{cases} 
\partial_t u = \nu_2 \Delta u + r_1 m^2 - \gamma(M)u - \tau_0 u, \\
\partial_t u_p = \gamma(M)u - \tau_p u_p, \\
\partial_t m = \nu_1 \Delta m + S(I, u) - dm - r_2 u m - r_1 m^2, \\
\partial_t M = D_1 \Delta M - \alpha \nabla \cdot (M \nabla u) + G(u)(\hat{M} - M)M - \sigma M + \lambda_M, \\
\partial_t I = D_I \Delta I + F(u)M - \tau_3 I, 
\end{cases} 
\quad (1) $$

The term $u$ corresponds to the concentration of the oligomers, while the term $m$ represents the concentration of the monomers. In this model, they simplify the assumption that oligomers correspond to the aggregate of two monomers, as in [27]. The term $u_p$ corresponds to the concentration of oligomers in amyloid plaques and is the main focus of our study. Finally, the terms $M$ and $I$ correspond to the concentration of microglial cells and interleukins respectively. The model is illustrated in figure 2.

To each component of the system, the authors associate a natural degradation term with a specific rate. Monomers, free oligomers, and interleukins are also assumed to diffuse with their own diffusion coefficients. In contrast, oligomers contained in amyloid plaques are considered inert and thus immobile. Microglial cells, however, are assumed to exhibit biased motion towards regions of high oligomer concentration, modeled through a combination of diffusion and chemotaxis. The other reaction terms are detailed below, equation by equation:

* In system (1), the equation for monomers is given by

$$ \partial_t m = \underbrace{\nu_1 \Delta m}_{\text{Diffusion}} + \underbrace{S(I, u)}_{\text{Production}} \underbrace{-r_2 u m - r_1 m^2}_{\text{Polymerization}} \underbrace{-dm}_{\text{Degradation}} $$
with the function $S(I, u)$ representing the production of monomers by neurons, defined by $S(I, u) := I\tau_S/(1 + Cu^n)$. The authors assume that monomer production is proportional to the local concentration of interleukins $I$, but that the actual amount produced also depends on the local concentration of oligomers. When the oligomer concentration is low ($u \approx 0$), production is essentially proportional to $I$; however, as oligomer concentration increases, the source term $S(I, u)$ tends to zero, modeling the *Unfolded Protein Response* (UPR). The two polymerization terms respectively correspond to the cases $u+m \to u$ and $m+m \to u$, based on the simplifying assumption previously mentioned.

* The equation for free oligomers is written as:

$$ \partial_t u = \underbrace{\nu_2 \Delta u}_{\text{Diffusion}} - \underbrace{\gamma(M) u}_{u \to u_p} + \underbrace{r_1 m^2}_{\text{Polymerization}} - \underbrace{\tau_0 u}_{\text{Degradation}} $$

where $\gamma(M)u$ denotes the recruitment of free oligomers into amyloid plaques, with the recruitment rate given by a function $\gamma$ that depends on the local concentration of microglial cells. The higher the microglial concentration, the greater the recruitment. The authors use a Michaelis-Menten type function: $\gamma(M) := \gamma_0 + \gamma_1 M/(1 + \gamma_2 M)$.

* The equation describing the concentration of plaque oligomers is:

$$ \partial_t u_p = \underbrace{\gamma(M)u}_{u \to u_p} - \underbrace{\tau_p u_p}_{\text{Degradation}} $$

* The equation governing the concentration of microglial cells is:

$$ \partial_t M = \underbrace{D_1 \Delta M}_{\text{Diffusion}} - \underbrace{\alpha \nabla \cdot (M \nabla u)}_{\text{Chemotaxis}} + \underbrace{G(u)(\hat{M} - M)M}_{\text{Proliferation}} + \underbrace{\lambda_M}_{\text{Production}} - \underbrace{\sigma M}_{\text{Degradation}} $$

The generation of new microglial cells arises either from proliferation or from a constant source term $\lambda_M$. Proliferation is assumed to follow logistic growth, with a proliferation rate given by a Michaelis-Menten type function $G(u)$. The motion of microglial cells is not purely Brownian: the authors assume a bias toward regions with high free oligomer concentration, modeled through the chemotactic term with coefficient $\alpha$, which will be particularly studied in this article.

* The equation for the concentration of interleukins is:

$$ \partial_t I = \underbrace{D_I \Delta I}_{\text{Diffusion}} + \underbrace{F(u)M}_{\text{Production}} - \underbrace{\tau_3 I}_{\text{Degradation}} $$

where the term $F(u)M$ represents the creation of new interleukins by microglial cells. This source term also depends on the local concentration of free oligomers through the function $F$, defined by the Michaelis-Menten expression: $F(u) := \tau_1 u/(1 + \tau_2 u)$. This dependence models the activation of microglial cells by free oligomers, which in turn leads to interleukin production. The parameters used in this model are listed in Table 1.

To simplify notation, we define the vector $V(x, t)$ as follows

$$ V(t, x) = (u(t, x), u_p(t, x), m(t, x), M(t, x), I(t, x)). $$

We study this problem in the one-dimensional case, i.e. $\Omega = [0, L]$, and also in the two-dimensional case, i.e. $\Omega = [0, L] \times [0, L]$, where $L > 0$ represents the size of the domain. Concerning the boundary conditions, we assume Neumann boundary conditions, i.e.

$$ \partial_{\mathbf{n}} V(y) = 0, \quad \text{for all } y \in \partial\Omega, \text{ where } \mathbf{n} \text{ denotes the unit normal to the boundary } \partial\Omega. $$
<table><caption>Table 1: Parameter values for the numerical simulations of the system (1).</caption><thead><tr><th>Parameter</th><th>Units</th><th>Description</th></tr></thead><tbody><tr><td>r<sub>1</sub></td><td>L (mol)<sup>-1</sup>(months)<sup>-1</sup></td><td>Bi-monomeric polymerization rate</td></tr><tr><td>r<sub>2</sub></td><td>L (mol)<sup>-1</sup>(months)<sup>-1</sup></td><td>Polymerization rate of monomers attaching to oligomers</td></tr><tr><td>d</td><td>(month)<sup>-1</sup></td><td>Degradation rate of monomers</td></tr><tr><td>γ<sub>0</sub></td><td>(month)<sup>-1</sup></td><td>Recruitment rate of oligomers to the amyloid plaques</td></tr><tr><td>τ<sub>1</sub></td><td>L (mol)<sup>-1</sup>(months)<sup>-1</sup></td><td>Production rate coefficient of interleukins induced by oligomers</td></tr><tr><td>τ<sub>2</sub></td><td>L (mol)<sup>-1</sup></td><td>Saturation constant for oligomer-induced interleukin production</td></tr><tr><td>τ<sub>3</sub></td><td>(months)<sup>-1</sup></td><td>Degradation rate of interleukins</td></tr><tr><td>τ<sub>p</sub></td><td>(months)<sup>-1</sup></td><td>Degradation rate of oligomers in the amyloid plaques</td></tr><tr><td>τ<sub>S</sub></td><td>(months)<sup>-1</sup></td><td>Coefficient of neural stress</td></tr><tr><td>C</td><td>L<sup>n</sup> (mol)<sup>-n</sup></td><td>Coefficient of neural stress function</td></tr><tr><td>n</td><td>-</td><td>Power coefficient of neural stress function</td></tr><tr><td>α<sub>1</sub></td><td>L<sup>2</sup> (mol)<sup>-2</sup>(months)<sup>-1</sup></td><td>Proliferation rate coefficient of microglial cells stimulated by oligomer concentration</td></tr><tr><td>α<sub>2</sub></td><td>L (mol)<sup>-1</sup></td><td>Saturation constant moderating the effect of oligomers on microglial proliferation</td></tr><tr><td>λ<sub>M</sub></td><td>mol L<sup>-1</sup>(months)<sup>-1</sup></td><td>Source term of microglial cells</td></tr><tr><td>M̃</td><td>mol L<sup>-1</sup></td><td>Carrying capacity of microglial cells</td></tr><tr><td>σ</td><td>(months)<sup>-1</sup></td><td>Degradation rate of microglial cells</td></tr><tr><td>D<sub>1</sub></td><td>m<sup>2</sup> (months)<sup>-1</sup></td><td>Diffusion coefficient of microglial cells</td></tr><tr><td>D<sub>l</sub></td><td>m<sup>2</sup> (months)<sup>-1</sup></td><td>Diffusion coefficient of interleukins</td></tr><tr><td>ν<sub>1</sub></td><td>m<sup>2</sup> (months)<sup>-1</sup></td><td>Diffusion coefficient of monomers</td></tr><tr><td>ν<sub>2</sub></td><td>m<sup>2</sup> (months)<sup>-1</sup></td><td>Diffusion coefficient of oligomers</td></tr><tr><td>α</td><td>m<sup>2</sup> L (mol)<sup>-1</sup></td><td>Coefficient of chemotaxis</td></tr></tbody></table>

In the article [27], the authors studied the spatially homogeneous steady states of the model (1). In particular, the point $V_0 := (0,0,0,\lambda_M/\sigma,0)$ is always an equilibrium point for the model (1), corresponding to the disease-free equilibrium, which we denote as $V_0$. Moreover, when the following condition is satisfied

$$ \sigma\gamma_0\tau_3 < \tau_1\tau_S\lambda_M, \qquad (2) $$

there is a critical threshold for the monomer degradation, denoted $d_m$, so that for $d < d_m$ the model (1) admits two other positive spatially-homogeneous steady states, denoted $V_u$ and $V_s$. In the absence of chemotaxis and diffusion, the first positive $V_u$ is unstable and the second one $V_s$ is stable.

## Objectives and results

In this article, we focus on the emergence of spatially periodic patterns corresponding biologically to the emergence of amyloid plaques. To achieve this goal, we study the possible instability of the positive equilibrium $V_s$ when the chemotaxis parameter $\alpha$ is strictly positive. First, we use a linear analysis to study the impact of the chemotaxis parameter on the stability of the equilibria in our problem. In particular, we prove in section 2 that the disease-free equilibrium is linearly asymptotically stable. Conversely, the positive equilibrium $V_s$ becomes unstable when the chemotaxis parameter is large enough. Secondly, we use numerical simulations to study the spatial patterns of the model (1). These patterns correspond to successions of zones of high and low oligomer concentration in the plaques.

We focus on the bistable regime of the model (1). Even if the analysis of the patterns of bistable systems is less studied, this type of numerical analysis has already been carried out for several bistable reaction-diffusion systems, notably for predatory prey systems or for systems describing chemical processes (see for example [28, 29, 30, 31]).

In the bistable regime, stable pattern branches compete with the stable homogeneous state. Therefore, it is not surprising to observe a convergence towards the stable equilibrium point. In the article by Krause *et al.* [32], the authors prove, through numerous examples, that mul-
tistability can significantly alter the predictions made by linear analysis. Indeed, they observe that the solution may converge towards a stable homogeneous state, even though the parameter values within the Turing space suggest the emergence of patterns. This phenomenon is also illustrated in [33], where Al-Karkhi et al. study the patterns of a predator-prey model with multiple stable homogeneous equilibria without diffusion. By conducting a bifurcation diagram, they show that the pattern branches can all lose their stability, thus explaining the absence of convergence towards a spatially heterogeneous pattern. In this article, we also observe and numerically analyze this phenomenon.

In the one-dimensional case, we find several possible scenarios: low-amplitude or high-amplitude stationary patterns, time-oscillating patterns, and convergence to the disease-free equilibrium. In the two-dimensional case, numerical simulations revealed the following phenomena: the emergence of stripe patterns, dots, complex stationary patterns between stripes and dots, and convergence to the disease-free equilibrium.

To study the existence and stability of these different possible patterns, we perform several bifurcation diagrams for the chemotaxis parameter $\alpha$, using the package pde2path [34, 35]. We notice a difference between these bifurcation diagrams and the typical ones, which are generally obtained for a monostable reaction-diffusion system. Generally, close enough to the Turing bifurcation, the only stable pattern is the one associated with the admissible critical frequency. Then moving away from this bifurcation, several other frequencies gradually become admissible and thus several stable branches emerge and coexist. In general, linear analysis shows that close to the Turing bifurcation, only a single mode can be destabilized, thereby predicting a stationary pattern with a spatial period corresponding to the critical frequency. Moving further away from this bifurcation, the study of the dispersion relation shows that several other modes progressively become admissible, and thus, several stable branches associated with these modes emerge and coexist. For our model (1), we notice an inverse trend: initially, several branches of stable patterns coexist, and the further we move away from the critical Turing threshold, the more the branches lose their stability, to the point where they all lose their stability. The result predicted by our numerical simulations is presented in the following conjecture:

**Conjecture 1.1 (Selection of disease-free equilibrium for $\alpha$ large enough)** We assume that the condition (2) is satisfied.

Then there exists $\alpha_0 > 0$ such that for all $\alpha > \alpha_0$,

$$ \lim_{t \to \infty} \mathbf{V}(x,t) \to \mathbf{V}_0, \text{ uniformly for all } x \in \Omega. $$

From a modeling point of view, this means counterintuitively that when the chemotaxis coefficient is large enough the system converges to the disease-free equilibrium without any amyloid plaques or inflammation.

In the case where $d > d_m$, the disease-free equilibrium is the only stable steady state and it is a global attractor for all initial data. Therefore, the value of $\alpha$ does not influence the selection of the final stationary state and we have $\alpha_0 = 0$. In the case where $d \le d_m$, the spatially homogeneous states have bistable dynamics and there is a critical threshold $\alpha_c$ for which the equilibrium $\mathbf{V}_s$ loses its stability. We observe numerically that there exists a $\alpha_0 \ge \alpha_c$ such that for $\alpha > \alpha_0$ the solution converges, independently of the slight noise, to the disease-free equilibrium. This threshold $\alpha_0$ depends on the values of the rest of the parameters in the model but also depends on the size of the domain $\Omega$. In particular, we also note that $\alpha_0 \to \alpha_c$ when $d \to d_m$.

Therefore, in the context of modeling, the parameter $d$, representing the natural degradation rate of monomers, is a reliable indicator of disease progression. When this parameter is sufficiently
large ($d > d_m$), natural degradation is sufficient to prevent the formation of amyloid plaques. In this case, the individual can be considered healthy or in a pre-symptomatic stage where the disease has not yet manifested. Conversely, when natural degradation is insufficient ($d < d_m$), plaque formation becomes possible, with a probability that depends on its value: the lower the value of $d$, the wider the range of parameter values for which plaques form and persist.

Our numerical simulations reveal that the stability range of pattern branches seems to be correlated with the frequencies of these patterns. For the examples given, we notice that the greater the number of peaks in a pattern, the faster its associated branch tends to lose stability as $\alpha$ increases. Based on these numerical results, we propose the following conjecture:

**Conjecture 1.2 (Decrease in the average number of amyloid plaques when $\alpha$ increases)**

We assume that the condition (2) is satisfied. Let $V(\cdot, t = 0)$ be an initial data corresponding to a slight Gaussian noise from the positive equilibrium $V_s$.

We assume that $\alpha > \alpha_c$, we denote $\#\max(\alpha)$ the number of local maxima of the final stationary solution of the system (1), with $\#\max$ defined as 0 in the spatially homogeneous case. The average of $\#\max(\alpha)$ over several independent realizations for different noisy initial conditions is denoted with some abuse of notation $E\#\max(\alpha)$. In this setting, we find the following property:

the function $E\#\max(\alpha)$ is non-increasing with respect to $\alpha$.

From a modeling perspective, this implies that as the chemotaxis coefficient increases, the number of amyloid plaques decreases.

The plan of the paper is the following: in Section 3, we begin with a linear analysis of the model (1) to obtain the conditions for pattern emergence, and in Section 4 we numerically study the stationary patterns generated by the model (1) in the one-dimensional case and the two-dimensional case.

# 2 Stability Analysis

## 2.1 Linear stability analysis of the disease-free equilibrium

We start by studying the stability of the disease-free equilibrium. The linearization around this steady state can be rewritten in a system of the form

$$ \frac{\partial \tilde{\omega}}{\partial t} = D \Delta \tilde{\omega} + J \tilde{\omega}, \qquad (3) $$

where $\tilde{\omega} = (\tilde{u}, \tilde{u}_p, \tilde{m}, \tilde{M}, \tilde{I})^\top$ is a vector for a small perturbation of the steady state $(0, 0, 0, \frac{\lambda}{\sigma}, 0)$ with Neumann boundary conditions and

$$ D = \begin{pmatrix} \nu_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \nu_1 & 0 & 0 \\ -\alpha \frac{\lambda}{\sigma} & 0 & 0 & D_1 & 0 \\ 0 & 0 & 0 & 0 & D_I \end{pmatrix}, \quad J = \begin{pmatrix} -\gamma_0 & 0 & 0 & 0 & 0 \\ \gamma_0 & -\tau_p & 0 & 0 & 0 \\ 0 & 0 & -d & 0 & \tau_S \\ (\hat{M} - \frac{\lambda}{\sigma}) \frac{\lambda}{\sigma} & 0 & 0 & -\sigma & 0 \\ \tau_1 \frac{\lambda}{\sigma} & 0 & 0 & 0 & -\tau_3 \end{pmatrix}. $$

In this setting, the corresponding eigenvalue problem is given by

$$ \lambda \tilde{\omega} = D \Delta \tilde{\omega} + J \tilde{\omega}. \qquad (4) $$
Concerning the stability, we have the following result.

**Proposition 2.1 (Disease-free equilibrium is asymptotically-stable)** *The eigenvalue problem (4) has a countable set of eigenvalues that are all strictly negative. Therefore the disease-free equilibrium of the system (1) is linearly asymptotically stable.*

**Proof.** Consider the following eigenvalue problem

$$ 
\begin{array}{l}
\Delta\varphi = \lambda\varphi \quad \text{in } \Omega, \\
\nabla\varphi \cdot \hat{n} = 0 \quad \text{on } \partial\Omega,
\end{array}
\qquad (5) $$

which has a set of eigenvalues {$\lambda_k$}$_{k \in \mathbb{N}}$ with $\lambda_k \le 0$ and $\lambda_k \to -\infty$ when $k \to \infty$. Let $\varphi_k$ be the corresponding eigenfunction.

We prove that the eigenvalues of problem (4) are strictly negative.

From Equation (5), we deduce that the solution of the first equation of the eigenvalue problem (4) given by

$$ 
\left\{
\begin{array}{ll}
\nu_2 \Delta \tilde{u} = (\lambda + \gamma_0) \tilde{u} & \text{in } \Omega, \\
\nabla \tilde{u} \cdot \hat{n} = 0 & \text{on } \partial \Omega,
\end{array}
\right.
\qquad (6) $$

is of the form $\tilde{u} = \varphi_k$ with $\lambda = \nu_2\lambda_k - \gamma_0 < 0$ and thus the system (4) has a solution for the corresponding eigenvector (take for instance $(\varphi_k, 0, 0, 0, 0)^\top$).

Assume now that $\lambda \ne \nu_2\lambda_k - \gamma_0$. From the theory of elliptic operators [38], we get that $\tilde{u} \equiv 0$ and the system (4) is reduced as follows.

$$ 
\left\{
\begin{array}{l}
(\lambda + \tau_p) \tilde{u}_p = 0, \\
\nu_1 \Delta \tilde{m} = (\lambda + d) \tilde{m} - \tau_S \tilde{I}, \\
D_1 \Delta \tilde{M} = (\lambda + \sigma) \tilde{M}, \\
D_I \Delta \tilde{I} = (\lambda + \tau_3) \tilde{I}.
\end{array}
\right.
\qquad (7) $$

By applying the same argument used for $\tilde{u}$ to rest of components of $\tilde{\omega}$ (first for $\tilde{u}_p$, $\tilde{M}$, $\tilde{I}$ and then for $\tilde{m}$), we conclude that set of eigenvalues $\mathcal{E}$ is given by

$$ \mathcal{E} = \{-\tau_p\} \cup \bigcup_{k \in \mathbb{N}} \{\nu_1\lambda_k - d, \nu_2\lambda_k - \gamma_0, D_1\lambda_k - \sigma, D_I\lambda_k - \tau_3\} \subset (-\infty, 0). $$

This proves the desired result. From the theory of elliptic equations [38], we remark that the spectrum of the linear operator determined by the system (4) consists only of eigenvalues.

□

## 2.2 Linear stability analysis of the positive equilibrium $V_s$

The corresponding linearization around the point $V_s$ gives the following equation

$$ \frac{\partial \tilde{\omega}}{\partial t} = D \Delta \tilde{\omega} + J \tilde{\omega}, \qquad (8) $$

where $\tilde{\omega} = (\tilde{u}, \tilde{u}_p, \tilde{m}, \tilde{M}, \tilde{I})^\top$ is a vector that corresponds to a small perturbation of the equilibrium $V_s$ with

$$ D = \begin{pmatrix}
\nu_2 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & \nu_1 & 0 & 0 \\
-\alpha M_0 & 0 & 0 & D_1 & 0 \\
0 & 0 & 0 & 0 & D_I
\end{pmatrix}, $$
$$ \text{and, } J = \begin{pmatrix} -\gamma_0 & 0 & 2r_1 m_0 & 0 & 0 \\ \gamma_0 & -\tau_p & 0 & 0 & 0 \\ -r_2 m_0 - I_0 \tau_S \frac{n C u_0^{n-1}}{(1+C u_0^n)^2} & 0 & -d - r_2 u_0 - 2r_1 m_0 & 0 & \frac{\tau_S}{1+C u_0^n} \\ (\hat{M} - M_0) M_0 \frac{\alpha_1}{(1+\alpha_2 u_0)^2} & 0 & 0 & \frac{\alpha_1 u_0}{1+\alpha_2 u_0} (\hat{M} - 2M_0) - \sigma & 0 \\ M_0 \frac{\tau_1}{(1+\tau_2 u_0)^2} & 0 & 0 & \frac{\tau_1 u_0}{1+\tau_2 u_0} & -\tau_3 \end{pmatrix}. $$

The solution of (8) with the form $\tilde{\omega} = e^{ikx+\lambda t}$ leads to the following equation,

$$ \lambda \tilde{\omega} = (J - k^2 D) \tilde{\omega}, \qquad (9) $$

where $\lambda \in \mathbb{C}$ is the eigenvalue and $k \in \mathbb{R}$ corresponds to the wavenumber. We denote $P_{J-k^2D}$, the characteristic polynomial associated with the problem (9). We obtain the following factorization, $P_{J-k^2D}(\lambda) = (-\tau_p - \lambda)Q_{J-k^2D}(\lambda)$. From this, we derive the first eigenvalue of the problem, equals to $\lambda_0 = -\tau_p < 0$ and the polynomial $Q_{J-k^2D}$ is given by the following relation

$$ Q_{J-k^2D}(\lambda) = \lambda^4 + a_1 \lambda^3 + a_2 \lambda^2 + a_3 \lambda + a_4, $$

with

$$ \begin{cases} a_1 = A + F + I + P, \\ a_2 = -BE + A(F + I + P) + F(I + P) + IP, \\ a_3 = -BGW + BE(I + P) - AFI - AFP - AIP - FIP, \\ a_4 = AFIP - BEIP + BGIW - BGKH(\alpha), \end{cases} $$

and where the constants are given by

$$ A = -\gamma_0 - k^2 \nu_2, \quad B = 2r_1 m_0, \quad E = -r_2 m_0 - I_0 \tau_S \frac{n C u_0^{n-1}}{(1 + C u_0^n)^2}, \quad F = -d - r_2 u_0 - 2r_1 m_0 - k^2 \nu_1, $$

$$ G = \frac{\tau_S}{1 + C u_0^n}, \quad H(\alpha) = (\hat{M} - M_0) M_0 \frac{\alpha_1}{(1 + \alpha_2 u_0)^2} + \alpha M_0 k^2, \quad I = \frac{\alpha_1 u_0}{1 + \alpha_2 u_0} (\hat{M} - 2M_0) - \sigma - k^2 D_1, $$

$$ K = \frac{\tau_1 u_0}{1 + \tau_2 u_0}, \quad P = -\tau_3 - k^2 D_1, \quad W = M_0 \frac{\tau_1}{(1 + \tau_2 u_0)^2}. $$

To study Turing bifurcation we first prove the following lemma.

**Proposition 2.2** Assume that there is a Turing bifurcation. Then the positive equilibrium $V_s$ becomes unstable when $\det(J - k^2D) > 0$ holds for at least one $k > 0$ and it remains stable when $\det(J - k^2D) < 0$ holds for all $k \ge 0$.

**Proof.** Using the Routh-Hurwitz criterion, we state that all eigenvalues have strictly negative real parts if and only if these five conditions hold

$$ a_i > 0 \text{ for all } i \in \{1, 2, 3, 4\} \text{ and } a_1 a_2 a_3 > a_3^2 + a_1^2 a_4. $$

We assume that $(\lambda_i)_{i \in \{1,2,3,4\}}$ represent the roots of the characteristic equation and accordingly to the characteristic equation, we obtain

$$ \begin{aligned} -a_1 &= \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4, & a_2 &= \lambda_1 \lambda_2 + \lambda_1 \lambda_3 + \lambda_1 \lambda_4 + \lambda_2 \lambda_3 + \lambda_2 \lambda_4 + \lambda_3 \lambda_4, \\ -a_3 &= \lambda_1 \lambda_2 \lambda_3 + \lambda_1 \lambda_2 \lambda_4 + \lambda_1 \lambda_3 \lambda_4 + \lambda_2 \lambda_3 \lambda_4, & a_4 &= \lambda_1 \lambda_2 \lambda_3 \lambda_4. \end{aligned} $$

Turing instability emerges when the real part of one eigenvalue passes through zero and the remaining eigenvalues have negative real parts. We denote $k_c$ as the critical wavenumber. Thus,
at this critical wavenumber $k = k_c$ we can assume without any loss of generality that $\lambda_1(k_c) = 0$ and $\lambda_2, \lambda_3, \lambda_4 < 0$.

Under this assumption, it follows that $a_1(k_c) > 0$, $a_2(k_c) > 0$, $a_3(k_c) > 0$ and $a_1(k_c)a_2(k_c) > a_3(k_c)$. However, for $k = k_c$, we obtain $a_4(k_c) = 0$. Thus, if there is a Turing bifurcation, that comes from the nullity of the term $a_4$ which also corresponds to the nullity of the determinant of the matrix $J - k^2D$. $\square$

The condition $a_4 < 0$ is equivalent to the following inequality

$$ \alpha > \frac{AFIP - BEIP + BGIW}{BGKM_0k^2} - \frac{(\hat{M} - M_0) \frac{\alpha_1}{(1+\alpha_2u_0)^2}}{k^2}, \quad (10) $$

by positivity of terms $B$, $G$, $K$ and $M_0$. Therefore, a spatial pattern appears on a one-dimensional domain that is sufficiently large, and for a sufficiently large chemotaxis parameter.

The term $a_4$ can be rewritten as a polynomial function of degree 4 in terms of $k^2$,

$$ a_4(k^2) = b_1(k^2)^4 + b_2(k^2)^3 + b_3(k^2)^2 + b_4(\alpha)k^2 + b_5, $$

with functions $b_i$, $i = 1, 2, 3, 4$, which have cumbersome expressions not reported in this work.

To determine whether a Turing bifurcation occurs, we analyze the minimal value of $a_4(k^2)$. For a given value of $\alpha$, we define $k_{\min}^2(\alpha)$ as the point where $a_4(k^2; \alpha)$ reaches its minimum. Turing instability arises when this minimum becomes zero, that is, when

$$ a_4(k_{\min}^2(\alpha_c); \alpha_c) = 0. $$

We then define $k_c^2 := k_{\min}^2(\alpha_c)$ as the critical wavenumber at bifurcation. Thus, a Turing bifurcation can occur at the critical value $\alpha_c$, which is given by a root of the function

$$ \alpha \mapsto b_1(k_{\min}^2(\alpha))^4 + b_2(k_{\min}^2(\alpha))^3 + b_3(k_{\min}^2(\alpha))^2 + b_4(\alpha)k_{\min}^2(\alpha) + b_5. \quad (11) $$

We compute this root using numerical tools. The results of these calculations are consistent with the numerical simulations and bifurcation diagrams presented in the next section.

## 3 Pattern analysis using numerical simulations

In this section, we use numerical tools to study the patterns induced by the loss of stability of the equilibrium $V_s$. The parameter values used in these simulations are similar to those presented in the article [27], whose orders of magnitude respect biological laws. These parameter values are given in Table 2. In this study, we only investigate the influence of the following variable parameters: chemotaxis coefficient, $\alpha$, and monomer degradation, $d$. First, we work on the domain $\Omega = [0, L]$ and then we perform a numerical study of the patterns for the two-dimensional square domain, $\Omega = [0, L] \times [0, L]$, $L > 0$.

The initial conditions chosen for our numerical simulations play a crucial role in pattern selection. Moreover, adding noise in these initial conditions is necessary for consistent modeling. Throughout this article, except for the Section 3.4, we choose a weak perturbation of the positive equilibrium state $V_s$, described in the remark below.

### Remark 3.1 (Choice of the noisy initial condition)

* For a one-dimensional domain, we discretize the space into the points $(i\Delta x)_i$, and for each grid point we add a slight Gaussian noise,

$$ V(i\Delta x, t = 0) = V_s + \epsilon N^i, $$
where $N^i$ corresponds to a Gaussian vector of dimension 5 and $\epsilon$ corresponds to a positive constant (e.g. $\epsilon = 10^{-5}$) small enough to preserve positivity.

* For the two-dimensional domain, we use the same initial conditions, but this time on the $(i\Delta x, j\Delta y)_{i,j}$ grid.

<table><caption>Table 2: Parameter values for the numerical simulations of the system (1).</caption><thead><tr><th>Parameter</th><th>Value</th><th>Units</th><th>Parameter</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>r<sub>1</sub></td><td>10<sup>-1</sup></td><td>L (mol)<sup>-1</sup>(months)<sup>-1</sup></td><td>α<sub>1</sub></td><td>1</td><td>L<sup>2</sup> (mol)<sup>-2</sup>(months)<sup>-1</sup></td></tr><tr><td>r<sub>2</sub></td><td>10<sup>-1</sup></td><td>L (mol)<sup>-1</sup>(months)<sup>-1</sup></td><td>α<sub>2</sub></td><td>1</td><td>L (mol)<sup>-1</sup></td></tr><tr><td>d</td><td>Variable</td><td>(month)<sup>-1</sup></td><td>λ<sub>M</sub></td><td>10<sup>-3</sup></td><td>mol L<sup>-1</sup>(months)<sup>-1</sup></td></tr><tr><td>γ<sub>0</sub></td><td>5 × 10<sup>-2</sup></td><td>(month)<sup>-1</sup></td><td>M&#770;</td><td>1</td><td>mol L<sup>-1</sup></td></tr><tr><td>τ<sub>1</sub></td><td>1</td><td>L (mol)<sup>-1</sup>(months)<sup>-1</sup></td><td>σ</td><td>10<sup>-3</sup></td><td>(months)<sup>-1</sup></td></tr><tr><td>τ<sub>2</sub></td><td>1</td><td>L (mol)<sup>-1</sup></td><td>D<sub>1</sub></td><td>1</td><td>m<sup>2</sup> (months)<sup>-1</sup></td></tr><tr><td>τ<sub>3</sub></td><td>1</td><td>(months)<sup>-1</sup></td><td>D<sub>I</sub></td><td>1</td><td>m<sup>2</sup> (months)<sup>-1</sup></td></tr><tr><td>τ<sub>p</sub></td><td>3 × 10<sup>-2</sup></td><td>(months)<sup>-1</sup></td><td>ν<sub>1</sub></td><td>1</td><td>m<sup>2</sup> (months)<sup>-1</sup></td></tr><tr><td>τ<sub>S</sub></td><td>1</td><td>(months)<sup>-1</sup></td><td>ν<sub>2</sub></td><td>1</td><td>m<sup>2</sup> (months)<sup>-1</sup></td></tr><tr><td>C</td><td>1</td><td>L<sup>n</sup> (mol)<sup>-n</sup></td><td>α</td><td>Variable</td><td>m<sup>2</sup> L (mol)<sup>-1</sup></td></tr><tr><td>n</td><td>2</td><td>-</td><td></td><td></td><td></td></tr></tbody></table>

## 3.1 Numerical simulation in the one-dimensional case

### Low and high amplitude patterns

For the one-dimensional domain, we numerically observe the existence of two different stationary pattern profiles. When the Turing bifurcation is supercritical, we notice the presence of low-amplitude patterns. In this case, the pattern corresponds to a sinusoidal perturbation of the positive equilibrium state with the critical frequency $k_c$. The amplitude of this perturbation increases with distance from the bifurcation. However, the stable branch associated with this pattern rapidly loses its stability. An example of a low-amplitude pattern is shown in figure 3 (left). When the bifurcation is subcritical, we notice the emergence of an unstable branch as expected. For our system and the parameter values tested, this branch remains unstable. Note that a weakly nonlinear analysis could be applied to approximate the stationary pattern in the supercritical case. This method can also determine the parameter values for which the bifurcation is sub- or supercritical (see, for example, [39, 40, 41]).

The nonlinear interactions of the system (1) can cause the emergence of large amplitude patterns. These correspond to higher peaks, distributed periodically in space, according to a certain frequency which is not the expected one and it is not even necessarily an admissible frequency according to the linear study. In figure 3 (right) we represent an example of high amplitude patterns. These high amplitude patterns are similar to biological images of amyloid plaques in the human brain, with the sequence of areas of high concentration and areas with the absence of beta-amyloid proteins. Note that these two types of low and high amplitude patterns can mathematically coexist.

### Time-oscillating patterns

We also notice the presence of heterogeneous patterns in space oscillating in time. Time-oscillating patterns have already been studied for classical reaction-diffusion systems [42, 43], for a prey-predator model [44], with a cross-diffusion term [46, 47, 48] and as well as for bistable reaction-diffusion models [29, 30].
<figure><img src="image_4.png" /><figcaption>Figure 3: Examples of final stationary patterns of the system (1). **Left.** Example of a low-amplitude pattern. **Right** Example of a high-amplitude pattern. The concentration of oligomers is shown in red, the concentration of monomers in magenta, the concentration of microglial cells in blue, the concentration of interleukins in green, and the concentration of oligomers in amyloid plaques in black. Parameter values are those given in Table 2, with $d = 0.05$ and $\alpha = 21.7$ for the figure on the left and $\alpha = 28$ for the figure on the right.</figcaption></figure>

An example of a time-periodic oscillating pattern is shown in figure 4A. This figure illustrates the evolution of oligomer concentration in amyloid plaques in space and time. The solution becomes periodic and oscillates between low-amplitude and high-amplitude patterns. In figure 4B, we illustrate the phase portrait of the concentration of oligomers at $x = 20$ as a function of the concentration of oligomers in the amyloid plaques. At the same time, we notice the emergence of a limit cycle in agreement with the periodicity in time of the solution. The dots shown in black and red correspond to the phase portrait at time $t = t_1$ and $t = t_2$ respectively and these two-time points are accordingly represented by two vertical lines of the same color in figure 4A. In figure 4C, we illustrate an example of the high amplitude in the oscillating pattern. This corresponds to the solution shown in figure 4A at time $t = t_1$. The solid curve shows the spatial distribution of oligomer concentration in amyloid plaques, while the dotted curve shows the evolution of microglial cell concentration. Finally, we also illustrate a low-amplitude pattern in figure 4D, taking a solution at time $t = t_2$. We compare it with the high-amplitude pattern, shown here in red.

Like in other articles (e.g. [46, 42]), we also notice the presence of quasi-periodicity and chaotic oscillations. For example, for a slightly larger value of $\alpha$ we notice the emergence of a torus in the phase portrait shown in figure 5.

A further slight increase in the magnitude of $\alpha$ drives the system towards more complex dynamics, as the resulting patterns become chaotic in both space and time. The time evolution of spatial averages for all components is quite irregular over time and exhibits specific chaotic characteristics. The chaotic time series are not shown here for brevity. However, the chaotic behavior is confirmed using two numerical tools: sensitivity to initial conditions and Lyapunov exponents. Figure 18 in the appendix demonstrates the divergence in the time evolution of
<figure><img src="image_5.png" /><figcaption>Figure 4: Spatially heterogeneous patterns that also exhibit regular oscillations in time. **A.** Concentration of oligomers in amyloid plaques over time and space. **B.** Evolution over time of the concentration of oligomers at $x = 20$ compared to the concentration of oligomers in the amyloid plaques also at $x = 20$. The black dot corresponds to time $t = t_1$, while the red dot corresponds to time $t = t_2$. **C.** Concentration of microglial cells and oligomers in amyloid plaques at time $t = t_1$. (dotted black line in sub-figure A.). **D.** Concentration of microglial cells (dotted line) and oligomers in amyloid plaques (solid line) at time $t = t_1$ in black and at time $t = t_2$ in red (vertical dotted red line in sub-figure A.).</figcaption></figure>

spatial averages for $u$ and $u_p$ when numerical simulations are performed with two slightly different initial conditions. Specifically, a small perturbation is introduced only at the center of the domain $x = 46$. The time evolution of the two simulation results is presented in red and black in Figure 18 in the appendix, clearly indicating the divergence of two nearby trajectories. Furthermore, the largest Lyapunov exponent for the spatiotemporal simulation results was calculated using the approach introduced in [49] and further employed in the literature [45, 50]. Positive largest Lyapunov exponents for the variables $u, u_p$, and $m$ are shown in Figure 19 in the appendix. The irregular oscillatory behavior of these variables at a fixed spatial point, as well as their spatial averages (presented in Figure 6), reinforces the chaotic oscillation of the constituent variables. The significant differences seen in the phase portraits in the two panels of Figure 6 confirm the irregularity of oscillations in both space and time.

Biologically speaking, to our knowledge, there is no evidence of temporal oscillation in amyloid plaque formation. Thus, we focus our work on time-stationary solutions. According to our numerical simulations, time-oscillating and low-amplitude exist only very close to the Turing bifurcation, which occurs at $\alpha = \alpha_c$. Therefore we mainly study the patterns for $\alpha$ assumed to be large enough to obtain only high amplitude. In the following section, we study in particular the influence of the parameter $d$ on stationary patterns.
<figure><img src="image_6.png" /><figcaption>Figure 5: Heterogeneous patterns in space with quasi-periodic oscillations in time. **Top.** Concentration of amyloid plaques over time. **Bottom-left.** Evolution over time of the concentration of oligomers at $x = 46$ as a function of the concentration of amyloid plaques also at $x = 46$. **Bottom-right.** Concentration of oligomers at $x = 46$ over time.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Illustration of chaotic phenomena. **Top.** Evolution over time of the spatial average of the concentration of oligomers as a function of the spatial average of the concentration of amyloid plaques. **Bottom.** Evolution over time of the concentration of oligomers at $x_0 = 46$ as a function of the concentration of amyloid plaques also at $x_0 = 46$.</figcaption></figure>

## Influence of monomer degradation on high-amplitude patterns

We successively study the formation of high-amplitude for the following parameter values: $d = 0.15$, $d = 0.35$ and $d = 0.45$.

* **For** $d = 0.15$. To study the stationary patterns of the system (1), we perform a bifurcation shown in figure 7 diagram by using the `pde2path` package [34, 35] and more specifically the code presented in the article [37] for a reaction-diffusion system including a chemotaxis term. Stable branches are represented by a thick line, while unstable branches are represented by a thin line and the black horizontal branch corresponds to the positive steady state $V_s$. As predicted in the section above, this branch loses its stability for $\alpha \ge \alpha_c$. At this Turing bifurcation, a branch is created which is associated with a periodic pattern in space. For $d = 0.15$, we notice that the Turing bifurcation is supercritical, indeed the branch is initially stable and associated with low amplitude, nearly sinusoidal pattern. The branch appears almost vertical near $\alpha_c$ in the bifurcation diagram, this is due to a rapid increase in amplitude over a narrow range of $\alpha$. As $\alpha$ increases, the pattern quickly evolves into a larger-amplitude solution with a more complex spatial structure, and the branch loses stability. Nevertheless, following several other bifurcations of this branch, we obtain new stable branches associated with high amplitude. First of all, we note the stability, over a short interval, of the branch associated with the pattern *One peak in the middle and one peak on the edge* represented by the blue line in figure 7. Then we notice two similar branches for the $L^1$ norm, associated respectively with the patterns *Middle peak* and *Two peaks at the edges*. The pattern with a middle peak can be seen as a reorganization of the pattern with two half-peaks at the edges, where each half-peak is shifted towards the center of the domain to form a complete, centered peak. Finally, we also notice the stability of the branch in magenta, corresponding to the pattern *One peak at the edge*. This diagram shows that for a wide range of $\alpha$, the system admits multi-stability of different branches, which is expected for models with a Turing-type bifurcation.
We notice that the size of the stability domain of the different branches of high amplitude decreases when the number of peaks or half-peaks (peaks on the edge) increases. Moreover, we notice that these domains of stability are finite and that for a large enough $\alpha$, none of the branches that we have detected is stable. For this bifurcation diagram, we have assumed a relatively small domain for the spatial variable, $\Omega = [0, 40]$, to have fewer possible stable patterns. For a larger domain, we notice again this negative correlation between the number of peaks and the length of the stability domain. Note that in addition to the stable branches shown in figure 7, the disease-free equilibrium, $V_0$, is also stable and can be an attractor. Therefore, when $\alpha$ is sufficiently large, this disease-free equilibrium becomes the unique stable steady state, as all other branches lose their stability. For smaller values of $\alpha$, $V_0$ coexists with other stable patterned states, contributing to the system's multistability.

<figure><img src="image_8.png" /><figcaption>Figure 7: Bifurcation diagram for $d = 0.15$. This bifurcation diagram is obtained using the package pde2path [34, 35]. Thick lines indicate stable branches, while thin lines indicate unstable branches. The horizontal branch in black corresponds to the stationary state $V_s$, which loses its stability in $\alpha_c$. For $d = 0.15$, this bifurcation point is supercritical, giving rise to a stable branch associated with a low-amplitude pattern. Following several bifurcations, we obtain the presence of other stable branches, associated with high-amplitude patterns described in the legend on the right. The term $M_s$ corresponds to the microglial cell concentration component of the positive equilibrium $V_s$. Parameter values are given in Table 2.</figcaption></figure>

In addition to the bifurcation diagram, which exhibits a multi-stability of different branches of possible patterns, we illustrate the probability of selection of these different stable structures in figure 8 by performing 200 independent simulations for different values of parameter $\alpha$. Pattern selection is highly dependent on the energy of each branch and the choice of initial condition. In our case, we choose an initial condition presented in Remark 3.1 with a very slight noise with an order of magnitude $10^{-5}$. Despite this relatively low noise level, we find at least once, a wide disparity. Indeed for several values of parameter $\alpha$, we notice three patterns out of the four determined using the bifurcation diagram. When $\alpha$ is relatively close to the critical threshold $\alpha_c$ the pattern *Two peaks at the edges* seems to be strongly selected. Then when $\alpha$ is between $\alpha = 40$ and $\alpha = 90$ the probabilities of selection of the *Two peaks at the edges*, *One peak at the edge* and *Middle peak* are in the same order of magnitude. Finally for $\alpha > 120$, in agreement with the Conjecture 1.1, the only possible stationary solution is the disease-free equilibrium.
In the top panel of figure 8, we also illustrate the domains where stationary solutions are stable. These domains corroborate the results shown at the bottom panel, obtained independently. Indeed, no pattern was selected outside its stability zone determined by the bifurcation diagram.

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

Figure 8: Pattern Selection for $d = 0.15$. **Top.** Stability of the different branches of stationary patterns obtained using the package *pde2path* for different values of parameter $\alpha$ (See figure 7), for $d = 0.15$. The black branch corresponds to the disease-free equilibrium, this branch is stable for all $\alpha \in \mathbb{R}$. **Bottom.** Percentage of obtaining the different stationary patterns following 200 numerical simulations independent of the (1) system. The noisy initial conditions used are described in Remark 3.1. This process was repeated for different values of parameter $\alpha$. The vertical black dashed line indicates the critical value $\alpha_c$ at which the Turing bifurcation occurs. The bifurcation diagram shown in figure 7 is consistent with our independent numerical simulations. Other parameter values are given in the Table 2.

In figure 9, we plot the average number of peaks in the final pattern associated with the numerical simulations in figure 8. We note that this average number of peaks tends to decrease as the chemotaxis coefficient, $\alpha$, increases. In particular, we note a strong discontinuity linked to the loss of stability of the *Two peaks at the edges* and *Middle peak* branches. These results are in agreement with the Conjecture 1.2.

* For $d = 0.35$. We now assume that the degradation coefficient of the monomers is equal to $d = 0.35$ and we perform the same analysis as before. The bifurcation diagram is shown in figure 10. This diagram has strong similarities with the previous diagram in figure 7. Again, close to the Turing-type bifurcation, we notice a strong multi-stability of different branches associated with different spatially heterogeneous patterns. The red branch corresponds to a pattern with two peaks on the edges, the green branch with one peak in the middle, and finally the blue branch corresponds to a single peak on the edge. As before, we notice that the length of the domains of stability is correlated with the number of half-peaks.

The main difference between the bifurcation diagram in figure 10 and the previous bifurcation
<figure><img src="image_10.png" /><figcaption>Figure 9: Average Number of peaks for $d = 0.15$. Average number of amyloid plaques for the final stationary pattern for 200 independent numerical simulations of the system (1), on the domain $\Omega = [0.40]$, for different values of the parameter $\alpha$. The noisy initial conditions used are described in Remark 3.1. The other parameter values are given in the Table 2.</figcaption></figure>

diagram shown in figure 7, concerns the length of the stable domains of the branches. In particular, we notice that the last branch to lose its stability, loses its stability close to $\alpha = 22$, much less than in the previous case.

To determine the probability of selection of these different possible patterns, we perform 50 simulations with a slightly noisy initial condition (see Remark 3.1), and the results are shown in figure 11. When the branch *Two peaks at the edges* is stable, we always select it even if other branches are also stable. As soon as the red branch loses stability, the blue branch *One peak at the edge*, becomes predominant in the selection of the stationary pattern. But in this range of value of $\alpha$ the solution can also, with low probability, converge towards the disease-free equilibrium. According to our bifurcation diagram, once the blue branch loses its stability, no branch associated with a spatially heterogeneous pattern is stable and this result agrees with our numerical simulations of figure 11. For parameter values $\alpha = 22$, 22.5, and 23, we systematically select the disease-free equilibrium.

The *Middle peak* structure was never selected in our simulations, despite its stability. Note that the noise is again chosen to be very low. For initial values with a larger variance, the results could differ greatly.

In the top panel of figure 11, we illustrate the average of the peaks of the selected stationary patterns presented in figure 10. As expected, we obtain a decreasing function as a function of $\alpha$, which is close to a step function, where the discontinuity comes from the loss of stability of a stable branch. For a larger domain, for example, $\Omega = [0, 100]$, more patterns are theoretically eligible. In this case, we notice that the decrease in the average of the peaks is less discontinuous than before (see the bottom panel of the figure 11).

* For $d = 0.45$. We finally study the case $d = 0.45$. For this parameter value, we are still in the bistable domain of the model 1, so there is a Turing type bifurcation such that for $\alpha > \alpha_c$ the stable equilibrium without chemotaxis $V_s$, becomes unstable. This bifurcation gives rise to a branch of spatially heterogeneous patterns. However, according to our bifurcation diagram and
<figure><img src="image_11.png" /></figure>

Figure 10: Bifurcation diagram for $d = 0.35$. This bifurcation diagram is obtained using the package `pde2path` [34, 35]. Thick lines indicate stable branches, while thin lines indicate unstable branches. The horizontal branch in black corresponds to the stationary state $V_s$, which loses its stability in $\alpha_c$. For $d = 0.35$, this bifurcation point is subcritical. Following several bifurcations, we obtain the presence of other branches, which are stable, and associated with high-amplitude patterns described in the legend on the right. The term $M_s$ corresponds to the microglial cell concentration component of the positive equilibrium $V_s$. Parameter values are given in Table 2.

independent numerical simulations, this branch is always unstable. So when the equilibrium $V_s$ loses its stability, the problem becomes monostable, and the solution of the model (1) converges towards the disease-free equilibrium.

## 3.2 Numerical simulation in the two-dimensional case

### Stripes and dots

We now study spatially heterogeneous patterns in a two-dimensional square domain, $\Omega = [0, L] \times [0, L]$. In this case, we can obtain many families such as stripes, dots, and hexagons, or even more complex ones, such as mazes.

Again we notice the presence of low-amplitude and high-amplitude. These low-amplitude, correspond to weak perturbations of the positive steady state $V_s$ and exist only for $\alpha$ very close to the critical threshold $\alpha_c$. Note that the weakly nonlinear analysis is a tool to theoretically determine the family of low-amplitude [51]. However, as in the previous section, we focus only on high amplitude. Our numerical simulations reveal several possible high amplitudes, including stripes and dots. These two types are illustrated in figure 13 and can coexist for certain ranges of parameter values. Note also that we obtain complex stationary patterns corresponding to a structure between dots and stripes. For a large enough $\alpha$, we only notice the presence of a dot. Subsequently, we focus on those which are more reasonable in terms of modeling. Indeed, biologically, amyloid plaques form spots distributed heterogeneously in space (see figure 1).
<figure><img src="image_12.png" /><figcaption>Figure 11: Pattern Selection for $d = 0.35$. **Top.** Stability of the different branches of stationary patterns obtained using the package *pde2path* for different values of the parameter $\alpha$ (See figure 10), for $d = 0.35$. The black branch corresponds to the disease-free equilibrium, this branch is stable for all $\alpha \in \mathbb{R}$. **Bottom.** Percentage of obtaining the different stationary patterns following 50 independent numerical simulations of the system (1). The noisy initial conditions used are described in Remark 3.1. This process was repeated for different values of parameter $\alpha$. The bifurcation diagram shown in figure 7 is consistent with our independent numerical simulations. Other parameter values are given in the Table 2.</figcaption></figure>

Pattern selection

In the two-dimensional case, there are many more possible structures. Therefore, it is difficult to detect all stable branches using a continuation method. In particular, there is a good chance of forgetting branches that are difficult to select but theoretically stable. Thus for the square domain $\Omega = [0, L] \times [0, L]$ we only work with successive independent numerical simulations.

In figure 14 we plot the evolution of the average number of peaks for 50 independent simulations of the model (1). The choice of the initial conditions is detailed in Remark 3.1. Since it is impossible to quantify the number of peaks for stripe patterns, we study the average only for $\alpha > 25$. For these values of $\alpha$ we always obtained stationary dot patterns. As for the study in the one-dimensional case, presented above, the average number of peaks tends to decrease when the chemotaxis parameter increases. From a certain stage, this average number of peaks is equal to 0. These results are in agreement with the Conjecture 1.2.
<figure><img src="image_13.png" /><figcaption>Figure 12: Average Number of peaks for $d = 0.35$. Average number of amyloid plaques for the final stationary pattern for 50 independent numerical simulations of the system (1), for different values of the parameter $\alpha$ and for the following domain: **Top.** $\Omega = [0, 40]$ and **Bottom.** $\Omega = [0, 100]$. The noisy initial conditions used are described in Remark 3.1. The other parameter values are given in the Table 2.</figcaption></figure>

### 3.3 Parameter value ranges without amyloid plaque

For modeling reasons, it is interesting to determine the zones of parameters for which it is *a priori* impossible to observe the emergence and persistence of amyloid plaques. In other words, for the initial conditions presented in Remark 3.1, we want to determine a threshold $\alpha_0$ such that for any $\alpha > \alpha_0$, the solution of the system (1) converges to the disease-free equilibrium independently of the noise of the initial condition.

We approach this threshold $\alpha_0$ using numerical simulations. We discretize the possible values of $\alpha$ and look for the smallest possible point such that for any $\alpha$ greater than this point the n independent simulations have led to convergence to the disease-free equilibrium. Our algorithm first explores a possible candidate and then refines this candidate to obtain a precise approximation of the threshold $\alpha_0$. We apply this algorithm for different values of parameter $d$ and
<figure><img src="image_14.png" /><figcaption>Figure 13: Examples of high-amplitude final stationary patterns for the system (1). **Top.** Stripe pattern. **Bottom.** Dot pattern. The parameter values are those shown in the Table 2, with $d = 0.15$ and $\alpha = 24$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Average Number of peaks for $d = 0.15$ and $\Omega = [0, 40] \times [0, 40]$. Average number of amyloid plaques for the final stationary pattern for 50 independent numerical simulations of the system (1), for different values of the parameter $\alpha$. For $\alpha > 25$, we did not obtain any stripes pattern for our various numerical simulations. The noisy initial conditions used are described in Remark 3.1. The other parameter values are given in the Table 2.</figcaption></figure>

obtain an approximation of the threshold $\alpha_0^{\exp}(d)$. For $d < d_m$, and if we assume that $d$ is large enough, the threshold $\alpha_0$ is equal to the critical threshold $\alpha_c$. Conversely, for $d$ small enough, we observe a strong increase in the $\alpha_0$ threshold.

In the one-dimensional case, we represent the approximation obtained using our algorithm, $\alpha_0^{\exp}(d)$, by the dotted red line marked by the red arrow in figure 15. In the bistable domain, we note that this curve separates two distinct parts. The first part corresponds to parameter values for which the solutions of our $n$ consecutive independent simulations have converged toward the disease-free equilibrium. We therefore assume that for these parameter values, we always converge to the disease-free equilibrium. Note that in theory, this assumption could turn out
to be wrong, especially for small parameter values *d* where the threshold seems very high. We denote $\mathbf{P}_d(\alpha)$ the probability that the solution does not converge to the equilibrium without disease. Even if this probability is positive for any $\alpha > 0$, we theoretically get the existence of $\alpha_0^{\exp}$. Nevertheless, the bifurcation diagrams presented in section 3.1 show that pattern branches lose their stability for a large enough $\alpha$; this motivates and justifies the choice of such an assumption. From a modeling point of view, this zone of parameter values corresponds to a zone where, after a certain time $t$, no amyloid plaques are present.

In the second part, shown in red, there are several possible phenomena. In particular for $\alpha < \alpha_c$, for very slightly noisy initial conditions around $\mathbf{V}_s$, the solution converges to the positive homogeneous stationary state $\mathbf{V}_s$. When $\alpha$ is between $\alpha_c$ and $\alpha_0^{\exp}$, there are plenty of possibilities, including convergence to a spatially heterogeneous stationary pattern, convergence to a time oscillating pattern and convergence to the disease-free equilibrium. Nevertheless, according to our numerical simulations, we obtained at least once a solution that is not the disease-free equilibrium.

In the monostable domain ($d > d_m$), the positive equilibrium $\mathbf{V}_s$ does not exist. For this range of parameter values, only the disease-free equilibrium is stable. Consequently, the solution converges, regardless of the chemotaxis parameter value and independently of the noise of the initial condition, to this disease-free equilibrium. In the two-dimensional case, i.e. for $\Omega = [0, L] \times [0, L]$ we obtain a similar result, which is represented in the figure 20 in Appendix A.

<figure><img src="image_16.png" /><figcaption>Figure 15: Asymptotic behavior of the solutions of system (1) as a function of the degradation of monomers *d* and the chemotaxis parameter $\alpha$, for the domain $\Omega = [0.40]$. For $d > d_m$, there is only one equilibrium which is the disease-free equilibrium. From the stability of this equilibrium, the solutions of the system converge towards it. For $d < d_m$, it is possible to converge to the disease-free equilibrium, the positive equilibrium $\mathbf{V}_s$, and a spatially heterogeneous (possibly time-oscillating) pattern. The red zone represents an area of parameter values for which we have noticed, at least once, convergence towards a stationary solution other than the disease-free equilibrium. The green zone corresponds to parameter values for which we systematically converge towards the disease-free equilibrium. This result is highly dependent on the chosen initial conditions. The choice of the noisy initial condition is described in Remark 3.1.</figcaption></figure>
## 3.4 Influence of the initial data

We assume that the initial concentrations of monomers, oligomers, and microglial cells are initially uniformly distributed in space and there is a complete absence of amyloid plaques. In this section, we investigate the influence of initial inflammation on the final pattern. The initial data is provided in Table 3.

<table><caption>Table 3: Initial data used in Section 3.4</caption><thead><tr><td>Parameter</td><td>Value</td><td>Units</td><td>Description</td></tr></thead><tbody><tr><td>u<sub>0</sub></td><td>10<sup>-4</sup></td><td>mol L<sup>-1</sup></td><td>Concentration of free oligomers</td></tr><tr><td>u<sub>p</sub></td><td>0</td><td>mol L<sup>-1</sup></td><td>Concentration of oligomers in the amyloid plaques</td></tr><tr><td>m<sub>0</sub></td><td>10<sup>-3</sup></td><td>mol L<sup>-1</sup></td><td>Concentration of monomers</td></tr><tr><td>M<sub>0</sub></td><td>1</td><td>mol L<sup>-1</sup></td><td>Concentration of microglial cells</td></tr><tr><td>I<sub>0</sub></td><td>Variable</td><td>mol L<sup>-1</sup></td><td>Concentration of interleukins</td></tr></tbody></table>

### 3.4.1 Initial homogeneous inflammation in space

First, we assume that the inflammation, corresponding to the interleukin concentration, is initially uniformly distributed in space with a certain quantity $I_0$ that we vary.

In figure 16, we follow a similar approach as previously done to identify parameter regions where the persistence of amyloid plaques is not observed. Here, we examine the influence of monomer degradation $d$ and the initial value of interleukin concentration $I_0$. The algorithm used to obtain this figure is similar to the one presented in Section 3.3. In the monostable case ($d > d_m$), where the only equilibrium is the disease-free equilibrium, all the solutions of the system (1) converge to this equilibrium, while in the bistable case ($d < d_m$), the steady-state solution can be either homogeneous or spatially heterogeneous. The green region corresponds to parameter values for which the steady-state solution of the model (1) coincides with the disease-free equilibrium. In contrast, the red region represents other possibilities, namely a homogeneous solution equal to the other positive equilibrium or a spatially heterogeneous solution. The dichotomy observed in figure 16 is similar to the results presented in the article [27].

### 3.4.2 Initial heterogeneous inflammation in space

We now assume that the initial concentration of interleukins is heterogeneous in space. More precisely, we assume that the interleukin concentration is initially distributed in three Gaussians evenly spaced across the domain $[0, 40]$. Mathematically, the initial data is defined by the following family of functions,

$$I_0^\eta(x) = \frac{1}{\eta\sqrt{2\pi}}e^{-\frac{(x-10)^2}{2\eta^2}} + \frac{1}{\eta\sqrt{2\pi}}e^{-\frac{(x-20)^2}{2\eta^2}} + \frac{1}{\eta\sqrt{2\pi}}e^{-\frac{(x-30)^2}{2\eta^2}}, \quad (12)$$

dependent on the variable $\eta > 0$, which corresponds to the variance of the Gaussians. Note that the initial mass of the interleukin concentration, $\int I_0$, is independent of the variable $\eta$.

In figure 17, we investigate the impact of the variance $\eta$ on the potential persistence of inflammation. In the left panel, we illustrate four different initial data of interleukin concentration, each associated with different variances, according to the formula (12). In the right panel, we depict the final steady-state solution associated with the corresponding initial data on the left. When the variance $\eta$ is small, the initial data corresponds to very narrow Gaussian curves, resulting in highly concentrated inflammation in small spatial regions. In this case, we observe that the solutions of the model (1) converge towards the disease-free equilibrium. Conversely,
<figure><img src="image_17.png" /></figure>

Figure 16: Asymptotic behavior of the solutions of system (1) as a function of the degradation of the monomers $d$ and the initial concentration of interleukin $I_0$. For $d > d_m$, there is only one equilibrium which is the disease-free equilibrium. By stability of this equilibrium, the solutions of the system converge towards this steady state. For $d < d_m$ it is possible to converge to the disease-free equilibrium, the positive equilibrium $V_s$, and a spatially heterogeneous (possibly time-oscillating) pattern. The red zone represents an area of parameter values for which we have noticed at least one convergence towards a stationary solution other than the disease-free equilibrium. The green zone corresponds to parameter values for which we systematically converge towards the disease-free equilibrium. This result is very similar to the one obtained in the article [27]. The other parameter values are given in Table 2 and the initial data is provided in Table 3.

when the variance $\eta$ is larger, the initial inflammation is more dispersed than before but with less amplitude. We notice that in this case the steady-state solutions of the model (1) correspond to spatially positive heterogeneous solutions and the initial inflammation leads to a persistent inflammation over time.

The result presented in figure 17 is in some sense quite similar to the proposition outlined in Conjecture 1.1. As $\alpha$ increases, heterogeneous steady-state solutions are more prone to converge to asymptotic Dirac masses. In other words, as $\alpha$ grows, the heterogeneous spatial concentration of interleukins becomes increasingly strong but also more concentrated in small regions of space, leaving a significant portion of space with an almost negligible interleukin density. In the previous sections, through numerical simulations, we proved that these spatially heterogeneous solutions lose their stability when $\alpha$ becomes excessively large ($\alpha > \alpha_0$), meaning when inflammation is too concentrated in a small fraction of space.

In terms of modeling, it is, therefore, simpler to eliminate inflammation if it is highly concentrated in a small part of the space, preserving a large part of the healthy tissue, rather than to treat milder inflammation spread over a wider area.
<figure><img src="image_18.png" /><figcaption>Figure 17: Steady-state solutions of the model (1) for various initial data corresponding to sums of Gaussians. **Left.** Different initial data for interleukin concentration given by the equation (12), associated with different variances from lower ones (Top) to higher ones (Bottom). **Right.** Steady-state solutions of the model (1) associated with the initial data from the left panel. The parameter values are listed in Table 2, and the remaining initial data is given in Table 3.</figcaption></figure>

# 4 Conclusion and perspectives

In this work, we have extended the research presented in the article by Ciuperca *et al.* [27] by studying the model (1) with its spatial components. Our linear analysis shows that the equilibrium point $V_s$ is unstable to spatially heterogeneous perturbations when $\alpha > \alpha_c$, with $\alpha_c$ being a critical threshold corresponding to a root of the function defined by (11). To study pattern formation, we plotted numerous bifurcation diagrams using the `pde2path` package, revealing a wide variety of patterns. We also performed numerous independent numerical simulations to ensure that all branches were found. With this complete bifurcation diagram, we can predict interesting modeling results. For instance, according to our model, the higher the coefficient $\alpha$, corresponding to the intensity of the attraction of microglial cells to oligomers, the fewer plaques are expected. Additionally, when $\alpha$ is sufficiently large, neither plaques nor inflammation are expected. This phenomenon also depends on other parameters; we specifically study the influence of the parameter $d$, representing the natural degradation rate of monomers.

However, to confirm or refute these predictions, it would be necessary to obtain brain images at different stages of plaque formation dynamics, as well as methods to modify parameters present in our system, such as the administration of anti-inflammatory drugs.

In our diagrams, we have highlighted that for sufficiently large values of $\alpha$, all pattern branches become unstable. This unusual behavior appears to be a consequence of the bistability of our model [32]. Similar results, showing a total loss of branch stability, were presented in the article by Al-Karkhi *et al.* [33]. Additionally, our bifurcation diagrams reveal that the loss of branch stability strongly depends on their associated pattern: the higher the number of peaks, the faster the loss of stability. It would be particularly interesting to theoretically study how bistability can be responsible for such phenomena.

The model (1) introduced in the article [27] is a simplified subsystem of a more complex model,
also introduced in the same article, which also models proto-oligomer structures of different sizes. Therefore, it would be interesting to study whether our results extend to the more general model. Additionally, the model was initially introduced for a study on a bounded open domain $\Omega \subset \mathbb{R}^3$. In this work, we performed analyses only in one and two dimensions, as these dimensions are simpler to handle. A three-dimensional analysis, although more complex, would be particularly interesting and is also made possible by the Matlab package `pde2path` (see, for example, [36]). In our work, amyloid plaque formation is modeled using a reaction-diffusion-chemotaxis system. However, other modeling frameworks can also produce spatial patterning relevant to Alzheimer's disease. For instance, phase-separation models based on Cahn-Hilliard-type equations may also lead to plaque-like structures with different long-term dynamics, such as merging or coarsening behaviors [52].

A direct application of our model is to study the effects of different therapeutic strategies. In this work, we have highlighted that microglial cells play a key role in plaque aggregation, particularly through the influence of the coefficient $\alpha$. Therefore, it would be interesting to identify factors, such as molecules, able to influence plaque structuring via microglial cells.

Anti-inflammatory treatments can also be studied using our model. As proven in the article [27], we show that initial inflammation may or may not lead to plaque formation. We also highlight that inflammation concentrated in a specific tissue area does not have the same impact as more diffuse inflammation. Given the strong multistability observed in this model, it is obvious that the administration of anti-inflammatory drugs induces transitions between different stable branches. Anti-inflammatory treatment has already been tested for Alzheimer's disease. In the article by Rivers-Auty et al. [53], the authors indicated that diclofenac could slow cognitive decline and offer promising prospects for the progression of Alzheimer's disease. However, rigorous clinical trials have failed to confirm the efficacy of these anti-inflammatory treatments [55, 54, 56]. In contrast, other studies suggest that the effectiveness of these treatments could be significant if administered well before the onset of clinical symptoms [57, 58].

# A Appendix: Supplementary Figures
<figure><img src="image_19.png" /><figcaption>Figure 18: Divergence of the spatial average from two simulations with very small change in initial condition. Divergence of the spatial average of $u_p$ (upper panel) and $u$ (lower panel) obtained from two simulations with very small change in initial condition for prey density at the spatial location $x = 46$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Determination of Lyapunov exponents for the variables $u$ (red), $u_p$ (black) and $m$ (blue). The largest Lyapunov exponent for the space-time simulation results was calculated using the approach introduced in [49].</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20: Asymptotic behavior of the solutions of system (1) as a function of the degradation of the monomers $d$ and the chemotaxis parameter $\alpha$, for the domain $\Omega = [0, 30] \times [0, 30]$. For $d > d_m$, there is only one equilibrium which is the disease-free equilibrium. By stability of this equilibrium, the solutions of the system converge towards it. For $d < d_m$, it is possible to converge to the disease-free equilibrium, the positive equilibrium $V_s$, and a spatially heterogeneous (possibly time-oscillating) pattern. The red zone represents an area of parameter values for which we have noticed, at least once, convergence towards a stationary solution that is not the disease-free equilibrium. The green zone corresponds to parameter values for which we systematically converge towards the disease-free equilibrium. This result is highly dependent on the chosen initial conditions. The choice of the noisy initial condition is described in Remark 3.1.</figcaption></figure>
References

[1] Haass, C. & Selkoe, D. Soluble protein oligomers in neurodegeneration: lessons from the Alzheimer's amyloid β-peptide. *Nature Reviews Molecular Cell Biology*. **8**, 101-112 (2007)

[2] Sakono, M. & Zako, T. Amyloid oligomers: formation and toxicity of Aβ oligomers. *The FEBS Journal*. **277**, 1348-1358 (2010)

[3] Sengupta, U., Nilson, A. & Kayed, R. The Role of Amyloid-β Oligomers in Toxicity, Propagation, and Immunotherapy. *EBioMedicine*. **6** pp. 42-49 (2016)

[4] Soto, C. Unfolding the role of protein misfolding in neurodegenerative diseases. *Nature Reviews Neuroscience*. **4**, 49-60 (2003)

[5] Cohen, S., Linse, S., Luheshi, L., Hellstrand, E., White, D., Rajah, L., Otzen, D., Ven- druscolo, M., Dobson, C. & Knowles, T. Proliferation of amyloid-β42 aggregates occurs through a secondary nucleation mechanism. *Proceedings Of The National Academy Of Sciences*. **110**, 9758-9763 (2013)

[6] Murphy, R. & Pallitto, M. Probing the Kinetics of β-Amyloid Self-Association. *Journal Of Structural Biology*. **130**, 109-122 (2000)

[7] Nag, S., Sarkar, B., Bandyopadhyay, A., Sahoo, B., Sreenivasan, V., Kombrabail, M., Muralidharan, C. & Maiti, S. Nature of the Amyloid-β Monomer and the Monomer-Oligomer Equilibrium. *Journal Of Biological Chemistry*. **286**, 13827-13833 (2011)

[8] Forloni, G. & Balducci, C. Alzheimer's Disease, Oligomers, and Inflammation. *Journal Of Alzheimer's Disease*. **62**, 1261-1276 (2018)

[9] Kinney, J., Bemiller, S., Murtishaw, A., Leisgang, A., Salazar, A. & Lamb, B. Inflamma- tion as a central mechanism in Alzheimer's disease. *Alzheimer's & Dementia: Translational Research & Clinical Interventions*. **4**, 575-590 (2018)

[10] Al-Ghraiybah, N., Wang, J., Alkhalifa, A., Roberts, A., Raj, R., Yang, E. & Kaddoumi, A. Glial Cell-Mediated Neuroinflammation in Alzheimer's Disease. *International Journal Of Molecular Sciences*. **23**, 10572 (2022)

[11] Lopategui Cabezas, I., Herrera Batista, A. & Pentón Rol, G. The role of glial cells in Alzheimer disease: potential therapeutic implications. *Neurología (English Edition)*. **29**, 305-309 (2014)

[12] Kreutzberg, G. Microglia: a sensor for pathological events in the CNS. *Trends In Neuro- sciences*. **19**, 312-318 (1996)

[13] Akiyama, H. Inflammation and Alzheimer's disease. *Neurobiology Of Aging*. **21**, 383-421 (2000,6)

[14] Kettenmann, H., Hanisch, U., Noda, M. & Verkhratsky, A. Physiology of Microglia. *Phys- iological Reviews*. **91**, 461-553 (2011)

[15] Andrade-Restrepo, M., Lemanre, P., Pujo-Menjouet, L., Tine, L. & Ciuperca, S. Modeling the spatial propagation of Aβ oligomers in Alzheimer's Disease. *ESAIM: Proceedings And Surveys*. **67** pp. 30-45 (2020)

[16] Ciuperca, I., Dumont, M., Lakmeche, A., Mazzocco, P., Pujo-Menjouet, L., Rezaei, H., Tine, L., Alzheimer's disease and prion: An *in vitro* mathematical model. *Discrete & Continuous Dynamical Systems - B*. **22**, 1-36 (2017)
[17] Hao, W. & Friedman, A. Mathematical model on Alzheimer's disease. *BMC Systems Biology*. **10**, 108 (2016)

[18] Matthäus, F. Diffusion versus network models as descriptions for the spread of prion diseases in the brain. *Journal Of Theoretical Biology*. **240**, 104-113 (2006)

[19] Matthäus, F. The spread of prion diseases in the brain — models of reaction and transport on networks. *Journal Of Biological Systems*. **17**, 623-641 (2009)

[20] Bertsch, M., Franchi, B., Marcello, N., Tesi, M. & Tosin, A. Alzheimer's disease: a mathematical model for onset and progression. *Mathematical Medicine And Biology: A Journal Of The IMA*. (2016)

[21] M. Bertsch, B. Franchi, M. C. Tesi, and V. Tora, The role of Aβ and Tau proteins in Alzheimer's disease: a mathematical model on graphs, *Journal of Mathematical Biology*, 87(3):49, 2023. doi:10.1007/s00285-023-01985-7.

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

[23] I. A. Kuznetsov & A. V. Kuznetsov, How the formation of amyloid plaques and neurofibrillary tangles may be related: a mathematical modelling study, *Proc. R. Soc. A*, 474(2210):20170777, 2018.

[24] I. A. Kuznetsov & A. V. Kuznetsov, Simulating the effect of formation of amyloid plaques on aggregation of tau protein, *Proc. R. Soc. A*, 474(2220):20180511, 2018.

[25] A. V. Kuznetsov, Effect of diffusivity of amyloid beta monomers on the formation of senile plaques, *Math. Med. Biol.*, 41(4):346-362, 2024.

[26] C. Chu, Y. L. C. Low, L. Ma, et al., How can we use mathematical modeling of amyloid-β in Alzheimer's disease research and clinical practices?, *Journal of Alzheimer's Disease*, 97(1):89-100, 2023. doi:10.3233/JAD-230938.

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

[28] Pearson, J. Complex Patterns in a Simple System. *Science*. **261**, 189-192 (1993)

[29] Mazin, W., Rasmussen, K., Mosekilde, E., Borckmans, P. & Dewel, G. Pattern formation in the bistable Gray-Scott model. *Mathematics And Computers In Simulation*. **40**, 371-396 (1996)

[30] Vanag, V. & Epstein, I. Out-of-phase oscillatory Turing patterns in a bistable reaction-diffusion system. *Physical Review E*. **71**, 066212 (2005)

[31] Weide Rodrigues, V., Cristina Mistro, D. & Díaz Rodrigues, L. Pattern Formation and Bistability in a Generalist Predator-Prey Model. *Mathematics*. **8**, 20 (2019)

[32] Krause, A., Gaffney, E., Jewell, T., Klika, V. & Walker, B. Turing instabilities are not enough to ensure pattern formation. (2023)

[33] Al-Karkhi, T., Kusdiantara, R., Susanto, H. & Codling, E. Bloom Formation and Turing Patterns in an Infochemical Mediated Multitrophic Plankton Model. *International Journal Of Bifurcation And Chaos*. **30**, 2030028 (2020)
[34] Uecker, H., Wetzel, D. & Rademacher, J. pde2path - A Matlab package for continuation and bifurcation in 2D elliptic systems. (2012)

[35] Uecker, H. Numerical continuation and bifurcation in nonlinear PDEs. Book. (2021)

[36] Uecker, H. & Wetzel, D. Snaking branches of planar BCC fronts in the 3D Brusselator. *Physica D: Nonlinear Phenomena.* **406** pp. 132383 (2020)

[37] Uecker, H. Pattern formation with pde2path – a tutorial. (2019)

[38] Gilbarg, D. & Trudinger, N. Elliptic Partial Differential Equations of Second Order. (Springer Berlin Heidelberg, 2001)

[39] Cross, M. & Hohenberg, P. Pattern formation outside of equilibrium. *Reviews Of Modern Physics.* **65**, 851-1112 (1993)

[40] Wollkind, D., Manoranjan, V. & Zhang, L. Weakly Nonlinear Stability Analyses of Prototype Reaction-Diffusion Model Equations. *SIAM Review.* **36**, 176-214 (1994)

[41] Gambino, G., Lombardo, M. & Sammartino, M. Turing instability and traveling fronts for a nonlinear reaction–diffusion system with cross-diffusion. *Mathematics And Computers In Simulation.* **82**, 1112-1132 (2012)

[42] Aragón, J., Barrio, R., Woolley, T., Baker, R. & Maini, P. Nonlinear effects on Turing patterns: Time oscillations and chaos. *Physical Review E.* **86**, 026201 (2012)

[43] Liu, R., Liaw, S. & Maini, P. Oscillatory Turing patterns in a simple reaction-diffusion system. (2007)

[44] Manna, K. & Banerjee, M. Spatiotemporal pattern formation in a prey–predator model with generalist predator. *Mathematical Modelling Of Natural Phenomena.* **17** pp. 6 (2022)

[45] Manna, K., Banerjee, M. & Petrovskii, S. On the structural sensitivity of some diffusion–reaction models of population dynamics. *Physica D: Nonlinear Phenomena.* **467** pp. 134220 (2024)

[46] Gambino, G., Lombardo, M. & Sammartino, M. Cross-diffusion-induced subharmonic spatial resonances in a predator-prey system. *Physical Review E.* **97**, 012220 (2018)

[47] Giunta, V., Lombardo, M. & Sammartino, M. Pattern formation and transition to chaos in a chemotaxis model of acute inflammation. *SIAM Journal On Applied Dynamical Systems.* **20**, 1844-1881 (2021)

[48] Mukherjee, N., Ghorai, S. & Banerjee, M. Cross-Diffusion Induced Turing and non-Turing Patterns in Rosenzweig–MacArthur Model. *Letters In Biomathematics.* **6**, 1-22 (2019)

[49] Morozov, A., Petrovskii, S. & Li, B. Spatiotemporal complexity of patchy invasion in a predator-prey system with the Allee effect. *Journal Of Theoretical Biology.* **238**, 18-35 (2006)

[50] Banerjee, M., Ghosh, S., Manfredi, P. & D'Onofrio, A. Spatio-temporal chaos and clustering induced by nonlocal information and vaccine hesitancy in the SIR epidemic model. *Chaos, Solitons & Fractals.* **170** pp. 113339 (2023)

[51] Ermentrout, B. Stripes or spots? Nonlinear effects in bifurcation of reaction—diffusion equations on the square. *Proceedings Of The Royal Society Of London. Series A: Mathematical And Physical Sciences.* **434**, 413-417 (1991)
[52] S. Sampaoli, A. Agosti, G. Pozzi, and P. Ciarletta, A toy model of misfolded protein aggregation and neural damage propagation in neurodegenerative diseases, *International Journal of Non-Linear Mechanics*, 144:104083, 2022.

[53] Rivers-Auty, J., Mather, A., Peters, R., Lawrence, C. & Brough, D. Anti-inflammatories in Alzheimer's disease—potential therapy or spurious correlate?. *Brain Communications*. **2**, fcaa109 (2020)

[54] Group, A. Cognitive function over time in the Alzheimer's disease anti-inflammatory prevention trial (ADAPT): results of a randomized, controlled trial of naproxen and celecoxib. *Archives Of Neurology*. **65**, 896 (2008)

[55] ADAPT Research Group Naproxen and celecoxib do not prevent AD in early results from a randomized controlled trial. *Neurology*. **68**, 1800-1808 (2007)

[56] Ozben, T. & Ozben, S. Neuro-inflammation and anti-inflammatory treatment options for Alzheimer's disease. *Clinical Biochemistry*. **72** pp. 87-89 (2019)

[57] Imbimbo Are NSAIDs useful to treat Alzheimer's disease or mild cognitive impairment?. *Frontiers In Aging Neuroscience*. (2010)

[58] Ali, M., Ghouri, R., Ans, A., Akbar, A. & Toheed, A. Recommendations for Anti-inflammatory Treatments in Alzheimer's Disease: A Comprehensive Review of the Literature. *Cureus*. (2019)