Microscopic and macroscopic models for the onset and progression of Alzheimer's disease

TL;DR

This paper reviews and connects microscopic and macroscopic mathematical models for Alzheimer's disease progression, demonstrating how particle-based models relate to continuum descriptions and showing qualitative agreement with clinical data.

Contribution

It establishes a link between microscopic neuron-level models and macroscopic population models, deriving the kinetic transport equation from particle dynamics.

Findings

Models qualitatively match clinical data

Microscopic and macroscopic models are consistent

Kinetic transport equation derived from particle models

Abstract

In the first part of this paper we review a mathematical model for the onset and progression of Alzheimer's disease (AD) that was developed in subsequent steps over several years. The model is meant to describe the evolution of AD in vivo. In [Y. Achdou et al., 2013] we treated the problem at a microscopic scale, where the typical length scale is a multiple of the size of the soma of a single neuron. Subsequently, in [M. Bertsch at al., 2016] we concentrated on the macroscopic scale, where brain neurons are regarded as a continuous medium, structured by their degree of malfunctioning. In the second part of the paper we consider the relation between the microscopic and the macroscopic models. In particular we show under which assumptions the kinetic transport equation, which in the macroscopic model governs the evolution of the probability measure for the degree of malfunctioning of neurons, can be derived from a particle-based setting. In the microscopic model we consider basically mechanism i), modelling it by a system of Smoluchowski equations for the amyloid concentration (describing the agglomeration phenomenon), with the addition of a diffusion term as well as of a source term on the neuronal membrane. At the macroscopic level instead we model processes i) and ii) by a system of Smoluchowski equations for the amyloid concentration, coupled to a kinetic-type transport equation for the distribution function of the degree of malfunctioning of the neurons. The second equation contains an integral term describing the random onset of the disease as a jump process localized in particularly sensitive areas of the brain. Even though we deliberately neglected many aspects of the complexity of the brain and the disease, numerical simulations are in both cases (microscopic and macroscopic) in good qualitative agreement with clinical data.

Figures

Equations

$$\Omega :={\Omega }_0\setminus \bigcup _{j=1}^M{\overline{\Omega }}_i.$$

$$\left\{\begin{array}[]{c}\partial_{\tau}u_{1}=d_{1}\Delta_{x}u_{1}-u_{1}\sum_{j=1}^{N}a_{1,j}u_{j}\\ \partial_{\nu}u_{1}=\psi_{0}\equiv 0\quad\mbox{on }\partial\Omega_{0}\\ \partial_{\nu}u_{1}=\psi_{j}\quad\mbox{on }\partial\Omega_{j},j=1,\dots,M\\ u_{1}(x,0)=U_{1}(x)\geq 0,\end{array}\right.$$

$$\left\{\begin{array}[]{c}\partial_{\tau}u_{m}=d_{m}\Delta_{x}u_{m}-u_{m}\sum_{j=1}^{N}a_{m,j}u_{j}+\frac{1}{2}\sum_{j=1}^{m-1}a_{j,m-j}u_{j}u_{m-j}\\ \partial_{\nu}u_{m}=0\quad\mbox{on }\partial\Omega_{0}\\ \partial_{\nu}u_{m}=0\quad\mbox{on }\partial\Omega_{j},j=1,\dots,M\\ u_{m}(x,0)=0,\end{array}\right.$$

$$\left\{\begin{array}[]{c}\partial_{\tau}u_{N}=d_{N}\Delta_{x}u_{N}+\frac{1}{2}\sum_{j+k\geq N,k<N,j<N}a_{j,k}u_{j}u_{k}\\ \partial_{\nu}u_{N}=0\quad\mbox{on }\partial\Omega_{0}\\ \partial_{\nu}u_{N}=0\quad\mbox{on }\partial\Omega_{j},j=1,\dots,M\\ u_{N}(x,0)=0,\end{array}\right.$$

$${\partial }_{\tau }{u}_1=d_1{\Delta }_x{u}_1-u_1\sum _{j=1}^N{a}_{1,j}{u}_j+\mathcal{F}$$

$$\mathcal{F}=\mathcal{F}[f]=C_{\mathcal{F}}{\int }_0^1({\mu }_0+a)(1-a)df(x,a,t).$$

$${\partial }_{t}f+{\partial }_a(fv[f])=0.$$

$$v[f]={\int }_0^1\mathcal{G}(x,a,b)df(x,b,t)+\mathcal{S}(x,a,u_1(x,\tau ),\dots ,u_{N-1}(x,\tau )).$$

$$\mathcal{G}(x,a,b)=C_{\mathcal{G}}(b-a{)}^{+},$$

$$\mathcal{S}=C_{\mathcal{S}}(1-a){\left(\sum _{m=1}^{N-1}m{u}_m(x,\tau )-\overline{U}\right)}^{+}.$$

$${\partial }_{t}f+{\partial }_a\left(fv[f]\right)=J[f],$$

$$J[f]=\eta \left({\int }_0^{1}P(t,a_{\ast }\to a)f(x,a_{\ast },t)d{a}_{\ast }-f(x,a,t)\right)\chi (x,t).$$

$$P(t,a_{\ast }\to a)\equiv P(a_{\ast }\to a)=\{\begin{array}{ll}\frac{2}{1-a_{\ast }}& \text{if }{a}_{\ast }\le a\le \frac{1+a_{\ast }}{2}\\ 0& \text{otherwise},\end{array}$$

$$\{\begin{array}{ll}{\partial }_{t}f+{\partial }_a\left(fv[f]\right)=J[f]& \\ \epsilon {\partial }_t{u}_1=d_1{\Delta }_x{u}_1-u_1\sum _{j=1}^N{a}_{1,j}{u}_j+\mathcal{F}[f]-{\sigma }_1{u}_1& \\ \epsilon {\partial }_t{u}_m=d_m{\Delta }_x{u}_m+\frac{1}{2}\sum _{j=1}^{m-1}{a}_{j,m-j}{u}_j{u}_{m-j}& \\ -u_m\sum _{j=1}^N{a}_{m,j}{u}_j-{\sigma }_m{u}_m& (2\le m<N)\\ \epsilon {\partial }_t{u}_N=\frac{1}{2}\sum _{\begin{array}{c}j+k\ge N\\ k,j<N\end{array}}{a}_{j,k}{u}_j{u}_k,& \end{array}$$

$$\{\begin{array}{ll}\text{vskip}3.0ptplus1.0ptminus1.0ptf(x,a,0)=f_0(x,a)& \text{if }x\in \Omega ,0\le a\le 1\\ \text{vskip}3.0ptplus1.0ptminus1.0pt{u}_i(x,0)=u_{0,i}(x)& \text{if }x\in \Omega ,1\le i\le N,\end{array}$$

$$-\frac{d_m}{\epsilon }{\Delta }_x{u}_m\cdot \mathbf{n}=0\text{on }\partial {\Omega }_{\text{out}},m=1,\dots ,N-1,$$

$$-\frac{d_m}{\epsilon }{\Delta }_x{u}_m\cdot \mathbf{n}=\gamma u_m\text{on }\partial {\Omega }_{\text{in}},m=1,\dots ,N-1,$$

$$(f,u_1,\cdots ,u_N)\in L^{\infty }(\Omega ;C([0,T];X_{[0,1]}))\times C([0,T];C^1(\overline{\Omega }){)}^N,$$

$${\int }_0^t\left(\int ({\partial }_s\varphi +v{\partial }_a\varphi )df(x,\cdot ,s)+\int \varphi dJ(x,\cdot ,s)\right)ds$$

$$=\int \varphi (x,\cdot ,t)df(x,\cdot ,t)-\int \varphi (x,\cdot ,0)d{f}_0(x,\cdot ).$$

$$H_{X,Y}\mathcal{G}(X,A_{\tau },B_{\tau }),$$

$$\mathbb{P}(H_{X,Y}=1)=h(X,Y),\mathbb{P}(H_{X,Y}=0)=1-h(X,Y)$$

$$\mathcal{S}=\mathcal{S}(X,A_{\tau },u(X,\tau )),$$

$$0\le J_{\tau }\le 1-A_{\tau },$$

$$\mathbb{P}(T_{\kappa }=1)=\kappa \Delta \tau ,\mathbb{P}(T_{\kappa }=0)=1-\kappa \Delta \tau ,\kappa =\nu ,\mu ,\eta ,$$

$$A_{\tau +\Delta \tau }=A_{\tau }+T_{\nu }{H}_{X,Y}\mathcal{G}(X,A_{\tau },B_{\tau })+T_{\mu }\mathcal{S}(X,A_{\tau },u(X,\tau ))+T_{\eta }{J}_{\tau }.$$

$$\phi (X,A_{\tau +\Delta \tau })=\phi (X,A_{\tau }+T_{\nu }{H}_{X,Y}\mathcal{G}(X,A_{\tau },B_{\tau })+T_{\mu }\mathcal{S}(X,A_{\tau },u(X,\tau ))+T_{\eta }{J}_{\tau }),$$

$$\displaystyle\begin{aligned} \left\langle\varphi(X,\,A_{\tau+\Delta{\tau}})\right\rangle=\left\langle\varphi(X,\,A_{\tau})\right\rangle+\Delta{\tau}\Bigl{[}&\nu\left\langle\varphi(X,\,A_{\tau}+H_{X,Y}\mathcal{G}(X,\,A_{\tau},\,B_{\tau}))\right\rangle\\ &+\mu\left\langle\varphi(X,\,A_{\tau}+\mathcal{S}(X,\,A_{\tau},\,u(X,\,\tau)))\right\rangle\\ &+\eta\left\langle\varphi(X,\,A_{\tau}+J_{\tau})\right\rangle\\ &-(\nu+\mu+\eta)\left\langle\varphi(X,\,A_{\tau})\right\rangle\Bigr{]}+o(\Delta{\tau}),\end{aligned}$$

$$⟨\phi (X,A_{\tau }+H_{X,Y}$$

$$=⟨\phi (X,A_{\tau }+\mathcal{G}(X,A_{\tau },B_{\tau }))h(X,Y)⟩$$

$$+⟨\phi (X,A_{\tau })(1-h(X,Y))⟩,$$

$$\frac{⟨\phi (X,A_{\tau +\Delta \tau })⟩-⟨\phi (X,A_{\tau })⟩}{\Delta \tau }$$

Full text

In the second part of the paper we consider the relation between the microscopic and the macroscopic models. In particular we show under which assumptions the kinetic transport equation, which in the macroscopic model governs the evolution of the probability measure for the degree of malfunctioning of neurons, can be derived from a particle-based setting.

The models are based on aggregation and diffusion equations for $\beta$ Amyloid, a protein fragment that healthy brains regularly produce and eliminate. In case of dementia A$\beta$ monomers are no longer properly washed out and begin to coalesce forming eventually plaques. Two different mechanisms are assumed to be relevant for the temporal evolution of the disease:

i) diffusion and agglomeration of soluble polymers of amyloid, produced by damaged neurons;

ii) neuron-to-neuron prion-like transmission.

In the microscopic model we consider basically mechanism i), modelling it by a system of Smoluchowski equations for the amyloid concentration (describing the agglomeration phenomenon), with the addition of a diffusion term as well as of a source term on the neuronal membrane. At the macroscopic level instead we model processes i) and ii) by a system of Smoluchowski equations for the amyloid concentration, coupled to a kinetic-type transport equation for the distribution function of the degree of malfunctioning of the neurons. The second equation contains an integral term describing the random onset of the disease as a jump process localized in particularly sensitive areas of the brain.

Even though we deliberately neglected many aspects of the complexity of the brain and the disease, numerical simulations are in both cases (microscopic and macroscopic) in good qualitative agreement with clinical data.

1. Introduction

The aim of the present paper is twofold: to provide an overview of the research carried on in the last few years by several authors in different and variated collaborations on both microscopic and macroscopic mathematical models for Alzheimer’s disease (AD) in the human brain [1, 12, 2, 10, 11], and to present a new result about the consistency of the microscopic and the macroscopic model. AD has a huge social and economic impact  [20, 4, 23]. Until 2040 its global prevalence, estimated as high as 44 millions in 2015, is expected to double every 20 years [35]. Not by chance AD-related issues belong to the cutting edge of scientific research. Apart from the classical in vivo and in vitro approaches, there is increasing interest in in silico research, based on mathematical modelling and computer simulations [40, 6, 13, 28, 8, 18, 17].

To cover the diverse facets of the AD in a single model, different spatial and temporal scales must be taken into account: microscopic spatial scales to describe the role of the neurons, macroscopic spatial and short temporal (minutes, hours) scales for the description of the relevant diffusion processes in the brain, and large temporal scales (years, decades) for the description of the global development of AD. In [1, 12] the authors began by attacking the problem at a microscopic scale, that is by considering as size scale of the model a multiple of the size of the soma of a single neuron (from 4 to 100 $\mu$m). Subsequently, see [2], the authors concentrated on a macroscopic scale, where they treat brain neurons as a continuous medium, and structure them by their degree of malfunctioning. Mathematically, the bridge between the two models is provided using two quite different techniques: through the so-called homogenisation technique in [10, 11]; and by adapting some arguments of the modern Boltzmann-type kinetic theory for multi-agent systems in Section  4 of this paper.

Following closely the biomedical literature on the AD, we briefly describe the processes (both microscopic and macroscopic) which we include in our models.

In the neurons and their interconnections several microscopic phenomena take place. It is largely accepted that beta amyloid (A$\beta$), especially its highly toxic oligomeric isoforms A$\beta_{40}$ and A$\beta_{42}$, play an important role in the process of the cerebral damage (the so-called amyloid cascade hypothesis [22]). In our papers we focus on the role of A$\beta_{42}$ in its soluble form, which recently has been suggested to be the principal cause of neuronal death and eventually dementia [41]. Indeed nowadays there are several evidences, such as enzyme-linked-immunosorbent assays (ELISAs) and mass spectrometry analysis, suggesting that the presence of plaques is not related to the severity of the AD. On the other hand, high levels of soluble A$\beta$ correlate much better with the presence and degree of cognitive deficits than plaque statistics. As a matter of fact some authors (see for instance [16]) overturn the traditional perspective, claiming that large aggregates of A$\beta$ can actually be inert or even protective to healthy neurons.

At the level of the neuronal membrane, monomeric A$\beta$ peptides originate from the proteolytic cleavage of a transmembrane glycoprotein, the amyloid precursor protein (APP). By unknown and partially genetic reasons, some neurons present an unbalance between produced and cleared A$\beta$ (we refer to such neurons as damaged neurons). In addition, it has been proposed that neuronal damage spreads in the neuronal net through a neuron-to-neuron prion-like propagation mechanism [5, 34].

On the other hand, macroscopic phenomena take place at the level of the cerebral tissue. The monomeric A$\beta$ produced by damaged neurons diffuses through the microscopic tortuosity of the brain tissue and undergoes a process of agglomeration, leading eventually to the formation of long, insoluble amyloid fibrils, which accumulate in spherical deposits known as senile plaques. In addition, soluble A$\beta$ shows a multiple neurotoxic effect: it induces a general inflammation that activates the microglia which in turn secretes proinflammatory innate cytokines [15] and, at the same time, increases intracellular calcium levels [13] yielding ultimately apoptosis and neuronal death.

The mathematical models which we derive in Sections 2 and 3 do not describe all the above-mentioned phenomena involved in the pathological process of the AD. They also neglect as well other additional phenomena, that we do not even mention. For example, we do not enter the details of the tortuosity of the brain tissue, we neglect the action of the $\tau$-protein, we simplify the role of the microglia, and neglect its multifaceted mechanism. In fact, we simply assume that high levels of soluble amyloid are toxic for neurons at all scales. Our primary goal was to overcome the fundamental mathematical difficulties and set the basis for a highly flexible model, which can be easily fine-tuned to include other issues. On the other hand, when we work at macroscopic scale we take into account also a neuron-to-neuron prion-like propagation mechanism ([34, 5]),

The models are minimal but effective: the numerical simulations produce a posteriori images and graphs which are in good qualitative agreement with clinical findings and confirm the validity of our assumptions. They also capture, at different scales, the cerebral damage in the early stage of the Mild Cognitive Impairment (MCI [33]).

In Section 4 we derive the equation for the progression of AD in the macroscopic model from a microscopic description of three main biophysical processes, among those recalled above. Namely, a prion-like spread of the disease over the neural network, the poisoning effect of soluble A$\beta$ polymers diffusing in the brain tissue and stochastic jumps in the level of neuron malfunctioning due to uncontrolled causes, such as e.g. external agents or genetic factors. We use mathematical techniques coming from the modern Boltzmann-type kinetic theory for multi-agent systems [32], such as microscopic binary interaction schemes and mean-field asymptotic limits.

In Section 5 we highlight some shortage of the present approach and we discuss possible extensions of the models, inspired by future research directions.

2. Mathematical model at the microscopic scale

When aiming at producing mathematical models of biological phenomena we have to fix preliminarily a spatial scale, as well as a time scale. Thus, we consider a portion of the cerebral cortex comparable in size to the size of a neuron, and we omit both the description of intracellular phenomena and clinical manifestations of the disease at a macroscopic scale, which will be considered instead in the model at the macroscopic scale. On the other hand, the experimental data of [25, 24] show that the process of plaques formation takes few days and therefore our temporal scale is chosen of the order of hours. In particular, no anatomical alteration of the neurons and of the surrounding cerebral tissue is taken into account.

The portion of cerebral tissue we consider is represented by a bounded smooth region $\Omega_{0}\subset\mathbb{R}^{3}$ (or $\Omega_{0}\subset\mathbb{R}^{2}$ in numerical simulations to reduce the computational complexity). To fix our ideas, we can think that the diameter of $\Omega_{0}$ is of the order of $10$ $\mu$m. The neurons are represented by a family of regular regions $\Omega_{j}$ such that

1. (1)

$\overline{\Omega}_{j}\subset\Omega_{0}$ if $j=1,\dots,M$; 2. (2)

$\overline{\Omega}_{i}\cap\overline{\Omega}_{j}=\emptyset$ if $i\neq j$.

We set

[TABLE]

To describe the evolution of the amyloid in $\Omega$, we consider a vector-valued function $u=(u_{1},\dots,u_{N})$, where $N\in\mathbb{N}$, $u_{m}=u_{m}(\tau,x)$, $m=1,\dots,N$, $x\in\Omega$ is the space variable and $\tau\geq 0$ is the time variable. If $1\leq m<N-1$ then $u_{m}(\tau,x)$ denotes the (molar) concentration at time $\tau\geq 0$ and point $x\in\Omega$ of A$\beta$ assemblies of polymers of length $m$. In addition, $u_{N}$ takes into account aggregations of more than $N-1$ monomers. Although $u_{N}$ has a different meaning from the other $u_{m}$’s, we keep the same letter $u$ in order to avoid cumbersome notations. Clusters of oligomers of length $\geq N$ (fibrils) may be thought of as a medical parameter (the plaques), clinically observable through PIB-PET (Pittsburgh compound B - PET [30]).

Coherently with this choice of the scales, it is coherent to assume that the diffusion of A$\beta$ in $\Omega$ is uniform, and therefore employ the usual Fourier linear diffusion equation (see, for instance, [29]).

In addition, we describe the agglomeration phenomena by means of the so-called finite Smoluchowski system of equations with diffusion. Classical references are [37, 7]; applications of Smoluchowski system to the description of the agglomeration of A$\beta$ amyloid appeared for the first time in [28].

The production of A$\beta$ in the monomeric form at the level of neuron membranes is modelled by a an inhomogeneous Neumann condition on $\partial\Omega_{j}$, the boundary of $\Omega_{j}$, for $j=1,\dots,M$. Finally, a homogeneous Neumann condition on $\partial\Omega_{0}$ is meant to neglect the neighbouring cerebral regions.

Thus, we are led to the following Cauchy-Neumann problem

[TABLE]

where $0\leq\psi_{j}\leq 1$ is a smooth function for $j=1,\dots,M$ describing the production of the amyloid near the membrane of the neuron.

We only take into account neurons affected by the disease, i.e. we assume $\psi_{j}\not\equiv 0$ for $j=1,\dots,M$. Moreover, to avoid technicalities, we assume that $U_{1}$ is smooth, more precisely $U_{1}\in\mathbf{C}^{2+\alpha}(\bar{\Omega})$ for some $\alpha\in(0,1)$, and that $\partial_{\nu}U_{1}=\psi_{j}$ on $\partial{\Omega_{j}}$, $j=0,\dots,M$.

In addition, if $1<m<N$,

[TABLE]

and

[TABLE]

where $d_{j}>0$, $j=1,\dots,N$ and $a_{i,j}=a_{j,i}>0$, $i,j=1,\dots,N$ (but $a_{N,N}=0$).

We assume that the diffusion coefficients $d_{j}$ are small when $j$ is large, since big assemblies do not move. In fact, the diffusion coefficient of a soluble peptide scales approximately as the reciprocal of the cube root of its molecular weight (see [14, 29]).

The form of the coefficients $a_{i,j}$ (the coagulation rates) considered in [28, formula (13)] rely on sophisticated statistical mechanics considerations (see also [19, 38]). In our numerical simulations, we use a slightly approximate form of these coefficients, taking $a_{i,j}=\dfrac{\alpha}{ij}$ where $\alpha>0$. In fact, this approximation basically consists in neglecting logarithmic terms in front of linear ones for large $i,j$. Concerning the coefficient $a_{N,N}$, it is clearly consistent with experimental data to assume $a_{N,N}=0$ for large $N$, which is equivalent to say that large oligomers do not aggregate with each other.

In our simulations we identify senile plaques with the sets $\{x\in\Omega\,:\,u_{N}(\tau,x)>c>0\}$. The following picture is provided by a numerical implementation of the model  (1), (2), (3). As in clinical observations, plaques grow near a neuron (the circle in the picture). The picture has been obtained by taking appropriate level sets of $u_{N}(\tau,\cdot)$.

We notice that our model leads to a smooth shape of the senile plaques (because of standard regularity properties of diffusion equations), in disagreement with evidences found in vivo. This may be explained by Figure 3 in [8] and related comments on the role of the microglia.

Besides numerical simulations, the main result obtained in [1] for this model is the following existence theorem:

Theorem 2.1**.**

For all $T>0$ the Neumann-Cauchy problem (1), (2), (3) has a unique classical positive solution $u\in\mathbf{C}^{1+\alpha/2,2+\alpha}([0,T]\times\bar{\Omega})$.

3. Mathematical model at the macroscopic scale

We identify now a large portion of the cerebral tissue with a 3-dimensional region $\Omega$, with diameter of $\Omega$ of the order of $10$ cm. As for the time scale, a new phenomenon occurs: two temporal scales are needed to simulate the evolution of the disease over a period of years, i.e. besides the short (i.e., rapid) $\tau$-scale (whose unit time coincides with hours) for the diffusion and agglomeration of A$\beta$ [24] that we used for the microscopic model, we need a long (i.e., slow) $t$-scale (whose unit time coincides with several months) to take into account the progression of AD. We can write the relation between the two scales as $\Delta{t}=\varepsilon\Delta\tau$ for a small parameter $\varepsilon\ll 1$.

At the macroscopic scale, the boundary vale problem for monomers (1) must have a different form. Indeed, the information given on the microscale by the non-homogeneous Neumann boundary condition is transferred into a source term $\mathcal{F}$ appearing in the macroscopic equation. This is due to the fact that at this scale neurons are reduced to points. Therefore, we have the following macroscopic equation for monomers:

[TABLE]

while the equations in (2) and (3) remain unchanged.

Mathematically, the transition from system (1) to equation (4) has been obtained by a two-scale homogenisation procedure described in [10] and [11].

The source term $\mathcal{F}$ in (4) will depend on the health state of the neurons. We model the degree of malfunctioning of a neuron with a parameter $a$ ranging from [math] to $1$: $a$ close to [math] stands for “the neuron is healthy” whereas $a$ close to $1$ for “the neuron is dead”. This parameter, although introduced for the sake of mathematical modelling (see also [34]), can be compared with medical images from Fluorodeoxyglucose PET (FDG-PET [27]).

For fixed $x\in\Omega$ and $t\geq 0$, let $f(x,\,a,\,t)$ be a probability measure, supported in $[0,1]$, that indicates the fraction of neurons close to $x$ with degree of malfunctioning between $a$ and $a+da$ at time $t$. From now on, we denote by $X_{[0,1]}$ the space of probability measures on $\mathbb{R}$ that are supported in $[0,1]$.

Since A$\beta$ monomers are produced by neurons and the production increases if neurons are damaged, we choose in (4)

[TABLE]

The small constant $\mu_{0}>0$ accounts for the physiologic A$\beta$ production by healthy neurons, and the factor $1-a$ for the fact that dead neurons do not produce amyloid.

The progression of AD occurs in the slow time scale $t$, over decades, and is determined by the deterioration rate $v=v(x,\,a,\,t)$ of the health state of the neurons through the continuity equation:

[TABLE]

Here $v[f]$ indicates that the deterioration rate depends on $f$ itself.

We assume that

[TABLE]

The notation $\mathcal{G}$ takes into account the spreading of the disease by proximity, while $\mathcal{S}$ models the action of toxic A$\beta$ oligomers, ultimately leading to apoptosis. For instance, we can choose

[TABLE]

and

[TABLE]

The threshold $\overline{U}>0$ indicates the minimal amount of toxic A$\beta$ needed to damage neurons, assuming that the toxicity of soluble $A\beta$-polymers does not depend on $m$. In reality length dependence has been observed [31], but, to our best knowledge, quantitative data are only available for very short molecules (see [31, Table 2]). For long molecules any analytic expression would be arbitrary.

At this point, we stress that equation (6), by its own nature, fails to describe the onset of the disease. To describe the onset of AD we assume that in small, randomly chosen parts of the cerebral tissue, concentrated for instance in the hippocampus, the degree of malfunctioning of neurons randomly jumps to higher values due to external agents or genetic factors. This leads to an additional term in the equation for $f$,

[TABLE]

where

[TABLE]

Here, $P(t,\,a_{\ast}\to a)$ is the probability to jump from the state $a_{\ast}$ to a state $a\in[0,\,1]$ (obviously, $P(t,\,a_{\ast}\to a)=0$ if $a<a_{\ast}$), $\chi(x,\,t)$ describes the random jump distribution, and $\eta$ is the jump frequency. For instance we can choose

[TABLE]

and neglect randomness, taking $\chi(x,\,t)\equiv\chi(x)$, concentrated in the hippocampus.

Finally, to model the phagocytic activity of the microglia as well as other bulk clearance processes [21], we add a term $-\sigma_{m}u_{m}$ in equations (2), (3) and (4), where $\sigma_{m}>0$.

We consider a transversal section (i.e. a horizontal planar section) of the brain that can be compared with radiological imaging (see, e.g., [2, Fig. 1]). For the sake of simplicity we schematise the section of the brain as a bounded connected region $\Omega\subset\mathbb{R}^{2}$, with two inner disjoint “holes” representing the sections of the cerebral ventricles. Consistently, we assume that the boundary of $\Omega$ consists of two disjoint parts: an outer boundary $\partial\Omega_{\text{out}}$ and an inner boundary $\partial\Omega_{\text{in}}$, i.e. the boundary of the cerebral ventricles consisting of two disjoint closed simple curves.

Eventually, we are led to the system

[TABLE]

with $\tau$ replaced by $\varepsilon^{-1}t$. Since we are interested in longitudinal modelling, we assume that initially, at $t=0$, there is a small uniform distribution of soluble amyloid $u_{0}=(u_{0,1},\dots,u_{0,N})$.

Thus system (12) has to be coupled with Cauchy initial data

[TABLE]

where the $u_{0,i}\in C^{1}(\overline{\Omega})$ are nonnegative functions for $i=1,\dots,N$, and $f_{0}\in L^{\infty}(\Omega;X_{[0,1]})$ describes the distribution of the disease at time $t=0$.

On the outer boundary $\partial\Omega_{\text{out}}$ we assume vanishing normal polymer flow. Therefore we impose a homogeneous Neumann condition for the diffusing amyloid polymers:

[TABLE]

$\mathbf{n}$ being the outward normal unit vector to $\partial\Omega_{\text{out}}$. Notice that no boundary condition is required for the concentration $u_{N}$ of the fibrillar amyloid, since its equation does not feature space dynamics (cf. the last equation in (12)). On the inner boundary $\partial\Omega_{\text{in}}$, that is the boundaries of the cerebral ventricles, we model the removal of A$\beta$ from cerebrospinal fluid (CSF) through the choroid plexus (cf. [21, 36]) by an outward polymer flow proportional to the concentration of the amyloid. For this, we impose a Robin boundary condition of the form:

[TABLE]

with $\gamma>0$ a constant.

An existence and uniqueness theorem for system (12) with Cauchy initial data (13) and boundary conditions (14) and (15) is proved in [3]. With our choices of $P$, $\mathcal{G}$ and $\mathcal{S}$ in (11), (8) and (9), it reads as follows:

Theorem 3.1**.**

For all $T>0$ there exist a unique $(N+1)$-ple

[TABLE]

$u_{i}\geq 0$* for $i=1,\dots,N$, solving (12) in a weak sense in $[0,T]$, with Cauchy data (13) and boundary data (14) and (15).*

In particular, the first equation in (12) is satisfied in the following weak sense: for a.e. $x\in\Omega$, for $\phi=\phi(x,\cdot,\cdot)\in\mathcal{D}(\mathbb{R}\times[0,T])$ and for all $t\in[0,T]$

[TABLE]

Concerning the outputs of the numerical simulation of (12) with Cauchy initial data (13) and boundary conditions (14) and (15), it is instructive to compare plots of $f$, at different times, with FDG-PET images (see e.g. [9]): we create a schematic image of a transverse section of the brain and attribute different colors (varying from red to blue) to those parts of the brain where probabilistically the level of malfunctioning lies in different ranges. As in the FDG-PET, the red corresponds to a healthy tissue. Here AD originates only from the hippocampus and propagates, at the beginning, along privileged directions (such as those corresponding to denser neural bundles) mimicked by two triangles. Obviously the details of the numerical output depend on the choice of the constants used in the mathematical model. Performing a considerable amount of numerical runs with different values of the constants in the model, we have reached the conclusion that, at least qualitatively, the behaviour of the solutions does not depend on the precise choice of those constants, as long as their variation is restricted to reasonable ranges.

The constant $\gamma$ enters the model through condition (15) at the boundary of the cerebral ventricles. Smaller values of $\gamma$ mean that less A$\beta$ is removed from the CSF through the choroid plexus. The comparison of the cases $\gamma=1$ and $\gamma=0.01$ becomes quite clear when we create spatial plots of $f$ (taking into account the average degree of malfunctioning of the brain in every point) at fixed computational times $t=T$. In Figures 2 and 3, where we compare plots of $f$ at, respectively, times $T=30$ and $T=40$ for the two different values of $\gamma=0.01$ and $\gamma=1$, AD originates only from the hippocampus and propagates, at the beginning, along privileged directions (such as those corresponding to denser neural bundles) mimicked by two triangles. Two remarks are now in order. First of all, though our images represent a mean value of brain activity instead of a single patient’s brain activity, still they show a good agreement with clinical neuroimaging (obviously representing the specific situation of an individual patient). See, e.g. [26], reproduced also with permission in [2], Fig. 6. The specificities (both anatomic and physiologic) of the single patient might account for the discrepancies between the outputs of our simulations and clinical neuroimaging. Secondly, we notice that if $\gamma$ becomes smaller (corresponding to a lower rate of clearance of the amyloid), we observe a temporal acceleration of the development of the illness; this could indicate the potential importance of the removal of $A\beta$ through the choroid plexus to slow down the temporal development of AD.

Looking for even more realistic images, we take now into account the randomness of the spatial distribution of the sources of the disease. Therefore we perform some runs where the AD does not only originate from the hippocampus, but also from several sources of A$\beta$ randomly distributed in the occipital part of the brain. We report the outputs of such runs, for $\gamma=1$ and two different values of time $T$, in Figure 4. The randomly distributed sources appear as the small blue spots.

Comparing the random sources case with the one in which AD originates only in the hippocampus, for the same values of $\gamma=1$ and $T=30$, it is clear that the brain sickness is more advanced when the number of sources is increased.

4. Derivation of the macroscopic equation for the progression of the disease

While the macroscopic system of Smoluchowski equations in (12) has been obtained from a smaller neuron-size scale through the homogenisation technique described in [10, 11], the mesoscopic equation for the distribution $f$ of the disease has been so far postulated on a mainly heuristic basis. In this section we provide its derivation from more fundamental particle-based dynamics, by adapting some arguments of the modern Boltzmann-type kinetic theory for multi-agent systems [32].

4.1. Particle-based neuron dynamics

Let $\tau\geq 0$ be the short (i.e. rapid) time variable, like in Section 3, and $\Omega\subset\mathbb{R}^{n}$ ($n=2,\,3$) a bounded subset of the physical space representing the brain, or possibly a two-dimensional section of it. We denote by $X\in\Omega$ the position of a neuron in the brain and by $A_{\tau}\in[0,\,1]$ its degree of malfunctioning at time $\tau$, which we assume to evolve according to the main biophysical mechanisms mentioned in the Introduction:

• •

a neuron-to-neuron prion-like transmission of the disease regarded as a binary interaction with a neighbouring neuron in the point $Y\in\Omega$ with degree of malfunctioning $B_{\tau}\in[0,\,1]$ at time $\tau$. We model the effect of such a binary interaction by a term of the form

[TABLE]

where $\mathcal{G}:\Omega\times[0,\,1]\times[0,\,1]\to[0,\,1]$ is a prescribed function which accounts for the prionic transmission of the disease and $H_{X,Y}\in\{0,\,1\}$ is a variable which describes the structure of the neural network. Specifically, $H_{X,Y}=1$ if the neurons in $X$ and $Y$ are connected by a synapse while $H_{X,Y}=0$ if they are not.

Due to the extremely complicated structure of the neural network, we consider a simple probabilistic description of it by assuming that $H_{X,Y}$ is a Bernoulli random variable parameterised by $X$, $Y$, in the sense that its law is

[TABLE]

for a given $h:\Omega\times\Omega\to[0,\,1]$ such that $h(x,\,y)=h(y,\,x)$ for all $x,\,y\in\Omega$.

On the other hand, as an approximation we explicitly disregard the variability in time of the connections among the neurons;

• •

the poisoning effect of soluble A$\beta$ polymers diffusing in the brain tissue, which we model by a function

[TABLE]

where we have set $u(x,\,\tau):=(u_{1}(x,\,\tau),\,\dots,\,u_{N-1}(x,\,\tau))$, $x\in\Omega$, for brevity;

• •

stochastic jumps in the degree of malfunctioning due to uncontrolled causes, such as external agents or genetic factors, which we model by means of a random variable $J_{\tau}$ such that

[TABLE]

as in no case the new degree of malfunctioning after the jump, i.e. $A_{\tau}+J_{\tau}$, can be greater than $1$.

In order to introduce a rule for the time variation of $A_{\tau}$ we assume that in a short time interval $\Delta{\tau}>0$ there is a probability proportional to $\Delta{\tau}$ that the neuron undergoes any of the mechanisms mentioned above. Furthermore, we assume that each mechanism is independent of the others. A simple way to formalise this is to introduce three independent Bernoulli random variables $T_{\nu},\,T_{\mu},\,T_{\eta}\in\{0,\,1\}$ such that

[TABLE]

where $\nu,\,\mu,\,\eta>0$ are the frequencies associated to each of the mechanisms above while the time interval has to be chosen in such a way that $\Delta{\tau}<1/\max\{\nu,\,\mu,\,\eta\}$. Under this assumption we set

[TABLE]

4.2. Boltzmann-type kinetic description

Let $\varphi=\varphi(x,\,a):\Omega\times[0,\,1]\to\mathbb{R}$ be a test function representing any observable quantity that can be computed out of the microscopic state $(X,\,A_{\tau})$ of a neuron. From (17) we get:

[TABLE]

whence, averaging both sides and computing first the mean with respect to the variables $T_{\nu},\,T_{\mu},\,T_{\eta}$,

[TABLE]

where $\left\langle\cdot\right\rangle$ denotes the average. Furthermore, using (16) to compute the mean with respect to $H_{X,Y}$ in the first term in brackets at the right-hand side, we obtain

[TABLE]

hence (18) specialises in

[TABLE]

In the limit $\Delta{\tau}\to 0^{+}$ this produces the continuous-in-time master equation

[TABLE]

Let us now introduce the probability density function

[TABLE]

of the microscopic state $(X,\,A_{\tau})$, i.e. $g(x,\,a,\,\tau)\,dx\,da$ is the fraction of neurons which at time $\tau$ are in the infinitesimal volume $dx$ centred at $x\in\Omega$ with a degree of malfunctioning in $[a,\,a+da]$. In the spirit of a Boltzmann-type ansatz, we assume that the processes $(X,\,A_{\tau})$ and $(Y,\,B_{\tau})$ are independent, so that their joint law is $g(x,\,a,\,\tau)g(y,\,b,\,\tau)$, cf. the next Remark 4.1. Moreover, we denote by $p(\tau,\,j|x,\,a)$, $0\leq j\leq 1-a$, the law of $J_{\tau}$ conditioned to $(X,\,A_{\tau})$, which is such that

[TABLE]

In view of these positions, we compute explicitly the remaining averages in (19) as:

[TABLE]

where the starred variables denote the state of the neuron after one of the three types of interactions according to (19):

[TABLE]

Equations (21), (22) provide the Boltzmann-type kinetic description of the microscopic model formulated in Section 4.1.

Remark 4.1*.*

Inspired by the discussion set forth in [32, Chapter 1], we observe that the assumption of stochastic independence of the states $(X,\,A_{\tau})$, $(Y,\,B_{\tau})$ is not fully justified from the biological point of view, being mostly dictated by the wish to obtain a closed equation in terms of the sole distribution function $g$. However, as it often happens in this type of problems, one needs to mediate between the high complexity of the biological phenomenon and the possibility to construct a usable, though necessarily approximated, mathematical model. In this respect, the aforesaid assumption should be regarded as a reasonable compromise, which permits a quite complete description and analysis of the evolution of the system.

4.3. The quasi-invariant degree of malfunctioning limit

As recalled in Section 3, the progression of AD occurs in a much slower time scale than that of the diffusion and agglomeration of A$\beta$ polymers. This implies that the actual time scale where the macroscopic effects of the progression of AD are observable is much longer than the $\tau$-scale at which the particle dynamics discussed in Sections 4.1, 4.2 take place. For this reason, as anticipated in Section 3 and inspired by the quasi-invariant interaction limits introduced in [32, 39], we now define the new time variable

[TABLE]

where $\varepsilon$ is a dimensionless parameter. Under this scaling, the typical time of a single particle transition (22), which is $O(1)$ in the $\tau$-scale, becomes much shorter in the $t$-scale, precisely $O(\varepsilon)$. Simultaneously, we scale by $\varepsilon$ also the interactions (22), considering that the effect of a single transition is attenuated in the longer $t$-scale. In particular, we set

[TABLE]

As far as the stochastic jumps are concerned, we assume instead that the strength of a single jump is independent of the time scale, hence we still have $a^{\ast\ast\ast}=a+j$ also in the $t$-scale, but the frequency $\eta$ of the jumps scales as $\varepsilon\eta$, i.e. single jumps are rarer, thus less probable, in the longer time scale.

In order to get from (21) an evolution equation in the $t$-scale, which avoids the detailed computation of the particle transitions in the unobservable $\tau$-scale, we introduce the scaled distribution function

[TABLE]

which satisfies the relations $\int_{0}^{1}\int_{\Omega}f(x,\,a,\,t)\,dx\,da=1$ and $\partial_{t}f=\frac{1}{\varepsilon}\partial_{\tau}g$, and, by (21), the equation

[TABLE]

where we have defined $P(t,\,j|x,\,a):=p(t/\varepsilon,\,j|x,\,a)$.

Taking $\varphi\in C^{\infty}(\Omega\times[0,\,1])$, with $\varphi(\cdot,\,0)=\varphi(\cdot,\,1)=0$, we expand

[TABLE]

with $\tilde{a},\,\bar{a}\in(0,\,1)$, then we plug into (21) to get

[TABLE]

where the remainder $R$ is

[TABLE]

Using that $0\leq h(x,\,y)\leq 1$ and $\int_{0}^{1}\int_{\Omega}f(x,\,a,\,t)\,dx\,da=1$, we see that

[TABLE]

Therefore, if $\mathcal{G}$ and $\mathcal{S}$ are bounded, $R(\varepsilon)\to 0$ as $\varepsilon\to 0^{+}$, and we obtain the equation

[TABLE]

for

[TABLE]

Let us further inspect the last term at the right-hand side of (26). Substituting $j$ with $a^{\ast\ast\ast}$ according to (22) yields:

[TABLE]

On the whole (26) rewrites as

[TABLE]

which can be finally recognised as a weak form of

[TABLE]

This equation is a balance law with non-local transport velocity given by (27). The term $J[f]$ at the right-hand side is the jump operator

[TABLE]

which, owing to (20), is such that $\int_{0}^{1}\int_{\Omega}J[f](t,\,x,\,a)\,dx\,da=0$. Denoting the probability law of the jumps as $P(t,\,a_{\ast}\to a|x)$ we see that it is the term modelled in (10), with the dependence on the spatial distribution of the jumps hidden in the dependence of $P$ on $x$.

Remark 4.2*.*

If we assume that the neural network is composed by an extremely large number of neurons, which are connected mainly with other close neurons, we can take the probability $h$ as $h(x,\,y)=\chi_{B_{1/N}(x)}(y)$, where $N$ is the total number of neurons and $B_{1/N}(x)\subset\Omega\subset\mathbb{R}^{n}$ is the ball centred in the point $x$ with radius $1/N$. If we suppose also $\nu\propto N^{d}$ (i.e. the more the neurons the more frequent the prionic transmission of the disease among them) then in the limit $N\to\infty$ we obtain that $\nu h(x,\,y)\to\delta_{0}(y-x)$ and

[TABLE]

which is indeed the form of $v$ postulated in (7).

Notice that, in general, $v$ depends on contributions from multiple time scales, in fact the prionic transmission of the disease (first term at the right-hand side of the formula above) takes place in the slower $t$-scale while the poisoning of the neurons by A$\beta$ polymers (second term at the right-hand side of the formula above) takes place in the faster $\tau$-scale.

5. Discussion

There are two main features of our macroscopic model that deserve to be highlighted: first of all, as we showed above, it derives mathematically from a microscopic model we carefully developed relying on an accurately selected set of phenomena described in biomedical literature. This has been possible also thanks to the interdisciplinary character of our team, including mathematicians and medical doctors. Secondly, the model is characterised by a high level of flexibility, which potentially allows one to simulate different working hypotheses and compare them with clinical data. In fact, the model provides a flexible tool to test alternative conjectures on the evolution of the disease. Up to now we have chosen some specific aspects of the illness, such as aggregation, diffusion and removal of $A\beta$, possible spread of neuronal damage in the neural pathway, and, to describe the onset of AD, a random neural deterioration mechanism. Some mathematical results have been obtained, and numerical simulations have been compared with clinical data. Although we have restricted ourselves to a 2-dimensional rectangular section of the brain, the results are in good qualitative agreement with the spread of the illness in the brain at various stages of its evolution. In particular, the model captures the cerebral damage in the early stage of MCI.

There are several research issues that we would like to address in the future. As mentioned in the Introduction, the main shortage of the present model is the complete omission of the action of the $\tau$ protein and the microglia (in this context we mention a mathematical model proposed in [18]). The mechanisms related to the presence of the two substances should eventually be considered in a subsequent evolution of this model, both to obtain optimal quantitative agreement with clinical data and also to investigate the possible formation of patterns in the distribution of the level of malfunctioning of the brain.

From a numerical point of view, simulations should become more realistic, in a 3-dimensional domain which matches the geometric characteristics of the brain. Moreover, a certain sensitivity of the numerical output to the value of the constant $\gamma$ in (15), which models the removal of $A\beta$ through the choroid plexus, spontaneously leads to the question whether dialysis-mechanisms can be introduced to enhance $A\beta$-removal artificially. Most probably, a serious answer to this question requires, in addition to a detailed comparison with clinical data, a more refined modelling of the clearance of soluble A$\beta$ by the cerebral fluid.

In [1] and [2] the reader can find an exhaustive discussion of the literature on mathematical models for AD. Besides that, in the recent paper [17] a large system of reaction-diffusion equations is proposed as a macroscopic model which takes into account many of the processes which possibly play a role in the development of AD. In addition the paper contains some simulations of medical treatments. The authors do not consider the onset of AD, and model the action of the $\beta$-amyloid in a way which is quite different from our approach. In a forthcoming paper we shall compare both approaches.

Acknowledgments

This paper is dedicated to Stu with deep affection and gratitude.

The authors would like to express their gratitude to MD Norina Marcello for many stimulating and fruitful discussions over several years.

B. F. and M. C. T. are supported by the University of Bologna, funds for selected research topics, and by MAnET Marie Curie Initial Training Network. B. F. is supported by GNAMPA of INdAM (Istituto Nazionale di Alta Matematica “F. Severi”), Italy, and by PRIN of the MIUR, Italy.

A. T. is member of GNFM (Gruppo Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di Alta Matematica “F. Severi”), Italy.