A Network Based Study of the Dynamics of Aβ and τ Proteins in Alzheimer’s Disease

1. Introduction

Alzheimer’s disease (AD) is a neurodegenerative disorder characterized by a progressive decline in memory and other cognitive functions, leading inevitably to death. It is an incurable disease affecting more than 50 million people, with figures set to increase significantly in the coming years (World Alzheimer Report 2023). The aetiology of the disease is unclear to date, but two proteins are universally recognised to play a crucial role in the development of the disease: amyloid-beta $\left(A\beta \right)$ and tau $\left(\tau \right)$. Both proteins are physiologically present in the brain, but in the presence of the disease they form abnormal aggregates, in a progressive and irreversible way. Indeed $A\beta$ plaques and neurofibrillary tangles (NFT) of $\tau$ are pathological hallmarks of AD. Experimental studies suggest non-uniform distributions of pathological proteins in the brain [1]. Both $A\beta$ and $\tau$ exhibit characteristic spatiotemporal deposition patterns. NFT tangles appear first in the entorhinal cortex and then spread to the amygdala, temporal areas and lately throughout the cortex [2]. On the other hand, $A\beta$ plaques form first in the temporal and frontal areas and then spread to other areas of the brain [2,3]. Recent literature suggests that the interplay between the two proteins should be crucial in the development of the disease and must be taken into account for the development of new therapies [4,5,6]. Clearly when one is interested in modelling complex pathologies such as AD, comparison with clinical and experimental data is crucial to the reliability of the study. Clinical multimodal neuroimages, if properly processed in terms of harmonization, combination and quantitative approach can give information of local character, providing for instance the concentration of proteins in different brain regions [7].

This work has two main objectives. One is to capture such types of non-uniform distributions through appropriate mathematical models. The other purpose is to compare the results from the mathematical models with clinical data in the context of Alzheimer’s disease. Comparison with clinical data is a fundamental step in the development of a trustable mathematical model but it is rarely addressed in the mathematical modelling literature. To this aim, we have processed a considerable amount of clinical neuroimages to extract quantitative features which can be compared with the $\tau$ concentration values provided by the mathematical models.

Networks are mathematical tools which provide an ideal setting for comparisons between clinical data and numerical results from models. Complex network properties have been identified with some consistency in all modalities of neuroimaging data and over a range of spatial and time scales. Thus, brain networks allow to describe the information process in neurons and its characteristic of being both locally specialized or segregated and globally distributed or integrated [8]. In particular, small-worldness characterizes an intermediate regime of networks between the two extreme regimes of a regular lattice and a random network, efficiently representing the brain structure [9]. Mathematical modeling of AD can be very useful for integrating clinical and experimental data in a theoretical setting which could allow, once reliability has been verified, to test various hypotheses otherwise difficult to verify. For an exhaustive overview of existing mathematical models of proteins spreading on networks and of various related challenging questions we refer to [10]. Here we wish to emphasise that the dynamics of the two proteins we are interested in differ in several respects, both spatial ($A\beta$ spreads over small distances, $\tau$ over large distances) and temporal (the dynamics of $A\beta$ is much faster than that of $\tau$). Differences in spatial dynamics should in some way be reflected by the choice of network used for the evolution of the dynamics. This paper is inspired by the paper [11] as far as the idea of using distinct networks for different protein dynamics is concerned; however its purpose is quite different. Indeed, in [11] the authors were mainly interested in investigating the synergisitic interactions of the two proteins $A\beta$ and $\tau$, whose relevance was confirmed there by testing various modelling hypotheses. However the authors were not interested to any comparison with medical data. Here we strongly rely on the comparison with medical data to test different modelling hypothesis concerning diffusion of $\tau$ on properly chosen networks, to see if there is an optimal way to deal with the spreading of proteins when it comes to modelling AD. Therefore the main target of this paper is twofold: 1. to consider different models for the spreading of $A\beta$ and especially of $\tau$, developed on networks with a precise biomedical meaning; 2. to have medical data for a comparison as meaningful as possible with the models, and use that comparison to test the trustability of the model.

The plan of the work is the following: in Section 2 we describe the networks we use for our models. We start from data publicly available, and from that we construct several different weighted graphs on which we evolve proteins. In Section 3 we describe the models we study: in particular, we present four different models for various possible ways of the spreading of $\tau$. Section 4 is concerned with the procedures we have adopted to obtain reliable medical data. We emphasize that we used a multimodal approach, based on MRI and PET, to determine the concentrations of $\tau$ protein in the brain. Finally in Section 5 we present the results obtained by numerical simulations and we make comparisons with medical data.

2. Brain networks via weigheted graphs

The central concept of a brain network constructed to segregate and integrate information processing comes from the advent of "disconnection syndromes" hypothesis [12] mainly on the basis of clinico-pathological correlations. The conceptualization of the relationship between structure and function of the brain has lead to understand that brain regions participate to many functions introducing the structure-function framework described by a network of brain regions. From Mesulam and over [13], it has been introduced the network approach to understand the localization of complex functions. The connection matrix of the human brain, known as the human ”connectome”, today represents an indispensable tool for mapping brain structure to functional processes and has a valuable impact on understanding brain diseases.

To model the spreading of different proteins in the brain, we consider several connectomes, each corresponding to a different undirected weighted graph. Data used for the construction of these connectomes have been downloaded at the website https://braingraph.org and consist of an averaged graph of 477 healthy subjects each with 1015 nodes (the weight calculation mode set as median and the number of fibers launched set to 20) [14,15]. We stress that the vertices of such a graph do not have associated coordinates that determine their position in space, since it is precisely an averaged graph. This will have consequences when we talk about “distances", which must be understood in an intrinsic sense and not in the classical euclidean metric sense.

In the following, we first recall the definition of weighted graph and how to construct an associated Laplacian, mimicking diffusion on it. Then, we describe the weighted graphs we will use to model the spreading of $A\beta$ and $\tau$ proteins. A graph is a pair $\mathcal{G}=(V,E)$, where V is a set of vertices and $E\subset V\times V$ is a set of the edges; $\mathcal{G}$ is said to be undirected if $(i,j)\in E$ implies that $(j,i)\in E$. The graph $\mathcal{G}$ is called weighted if a measure $\omega :E\to {\mathbb{R}}^{+}$ exists, assigning a unique positive number to each edge; the value $\omega (i,j)$ is referred to as weight of the edge $(i,j)$. A weighted graph can be represented through the adjacency matrix A whose entries $A_{i,j}$ represent the weights of the edge $(i,j)$. Let $N=\left|V\right|$ be the number of vertices in $\mathcal{G}$, then $A\in {\mathbb{R}}^{N\times N}$ is defined as

$$A_{i,j}=\left\{\begin{array}{cc}w(i,j)\hfill & \mathrm{if}(i,j)\in E;\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

There are several possible definitions for the Laplacian associated to a given graph. Given the adjacency matrix $A\in {\mathbb{R}}^{N\times N}$ of $\mathcal{G}$, following [16] we define the graph Laplacian L as $L=D-A$ where D is the weighted degree matrix whose jth diagonal element $D_{j,j}$ is given by

$$D_{j,j}=\sum _{i=1}^N{A}_{i,j}.$$

2.1. Structural Connectome

We call structural connectome a weighted graph $\mathcal{G}=(V,E)$ extracted from tractography of diffusion tensor images of 477 healthy subjects of the Human Connectome Project [17] using the Budapest Reference Connectome v3.0 [15]. In this graph, vertices correspond to parcellated regions of grey matter formed by neurons which share similarity in cytoarchitecture, functional activity, and structural connections to other regions, and edges represent the connectivity between the regions. We use a high-resolution connectome with $N=1015$, which can be downloaded at the website https://pitgroup.org/connectome/. It is possible to choose for the edges of this graph different weights, based on the mean number of fibers connecting two regions, and on their mean lenght. Accordingly, we will denote by ${\mathcal{G}}_{\mathrm{NL}}$ the graph with weights given by

$$w_{\mathrm{NL}}(i,j)={\displaystyle \frac{n_{ij}}{{\ell }_{ij}}},$$

(1)

where $n_{ij}$ is the mean number of fibers connecting verices i and j, and ${\ell }_{ij}$ is the mean lenght of such fibers. Analogously, we will denote by ${\mathcal{G}}_{\mathrm{L}}$ the graph with weights given by

$$w_{\mathrm{L}}(i,j)={\ell }_{ij}.$$

(2)

Therefore, the two weighetd graphs ${\mathcal{G}}_{\mathrm{NL}}$ and ${\mathcal{G}}_{\mathrm{L}}$ have the same vertices and edges, but different weights.

2.2. Intrinsic Proximity Connectome

Let $R_{\mathrm{P}}\in {\mathbb{R}}^{+}$ be a given positive value, the intrinsic proximity connectome is the weighted graph ${\mathcal{G}}_{\mathrm{P}}=(V,E)$ whose set of vertices V is the same as ${\mathcal{G}}_{\mathrm{L}}$’s and whose set of edges $E_{\mathrm{P}}\subset E_{\mathrm{L}}$ is the subset of edges of ${\mathcal{G}}_{\mathrm{L}}$ connected by a fiber with length ℓ less than $R_{\mathrm{P}}$. The weights are given by

$$w_{\mathrm{P}}(i,j)=\left\{\begin{array}{cc}{e}^{-{\ell }_{ij}^2/\sigma }\hfill & \mathrm{if}{\ell }_{ij}\le R_{\mathrm{P}}\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

where $\sigma \in {\mathbb{R}}^{+}$ is a fixed parameter. In this way, we connect two vertices only if they are close in an intrinsic sense, as opposed to a “geometric" vicinity measured with some kind of euclidean distance. Indeed, as already mentioned the vertices of the graph do not have associated spatial coordinates. The weights we assign are stronger for intrinsically close vertices. For this reason, ${\mathcal{G}}_{\mathrm{P}}$ is referred to as intrinsic proximity connectome. It is appropriate to emphasize here that this proximity connectome is quite different from the one used in [11]. Indeed here we use only intrisic connections, without adding any new connection, and we give them a weight proportional to the intrinsic distance.

2.3. Cumulative Connectome

The cumulative connectome is a weighted graph ${\mathcal{G}}_{\mathrm{C}}=(V,E)$ with the same V and E as ${\mathcal{G}}_{\mathrm{NL}}$, but with weights defined as follows. Let ${\mathcal{S}}_{ij}$ be the set of all paths in the graph ${\mathcal{G}}_{\mathrm{NL}}$ starting at vertex i and ending at vertex j, i.e:

$${\mathcal{S}}_{ij}=\{p_{ij}^1,p_{ij}^2,\dots ,p_{ij}^{m_{ij}}\},$$

(3)

where $m_{ij}=\left|{\mathcal{S}}_{ij}\right|$ and each path $p_{ij}^k\in {\mathcal{S}}_{ij}$ consists of a sequence of edges:

$$p_{ij}^k=\{e_{ij}^{k_1},e_{ij}^{k_2},\dots ,e_{ij}^{k_r}\},k=1,\dots ,m_{ij}.$$

Recall that the graphs ${\mathcal{G}}_{\mathrm{NL}}$ and ${\mathcal{G}}_{\mathrm{L}}$ differ only in the weights assigned to the edges, thus we can define the length of the path $p_{ij}^k$ considering it living in the graph ${\mathcal{G}}_{\mathrm{L}}$ or in the graph ${\mathcal{G}}_{\mathrm{NL}}$. Since the length of a path is usually defined as the sum of the weights of the edges forming the path, in ${\mathcal{G}}_{\mathrm{NL}}$ we obtain the length $l_{\mathrm{NL}}\left(p_{ij}^k\right)$ defined as

$$l_{\mathrm{NL}}\left(p_{ij}^k\right)=\sum _{s\in \{k_1,\dots ,k_r\}}{w}_{\mathrm{NL}}\left(e_{ij}^s\right),k=1,\dots ,m_{ij},$$

and, in ${\mathcal{G}}_{\mathrm{L}}$ we obtain the length $l_{\mathrm{L}}\left(p_{ij}^k\right)$ given by

$$l_{\mathrm{L}}\left(p_{ij}^k\right)=\sum _{s\in \{k_1,\dots ,k_r\}}{w}_{\mathrm{L}}\left(e_{ij}^s\right),k=1,\dots ,m_{ij}.$$

Finally let ${\mathcal{S}}_{ij}^{\delta }$ be the subset of ${\mathcal{S}}_{ij}$ such that

$${\mathcal{S}}_{ij}^{\delta }=\{p_{ij}^k\in {\mathcal{S}}_{ij}|l_{\mathrm{L}}\left(p_{ij}^k\right)\le \delta \}$$

(4)

where $\delta \in {\mathbb{R}}^{+}$ is a fixed parameter. Two vertices i and j are connected in the cumulative connectome ${\mathcal{G}}_{\mathrm{C}}$ if there is at least one path of length in ${\mathcal{G}}_{\mathrm{L}}$ less or equal $\delta$ joining them (we call such a path an admissible path.). In the affirmative case, the weight of the edge connecting vertices i and j is given by the sum of the lengths in ${\mathcal{G}}_{\mathrm{NL}}$ of all admissible paths (hence the name cumulative connectome):

$${\displaystyle w_{\mathrm{C}}(i,j)=\sum _{p_{ij}^k\in {\mathcal{S}}_{ij}^{\delta }}{l}_{\mathrm{NL}}\left(p_{ij}^k\right).}$$

Summarising, in the cumulative connectome, two vertices i and j are connected if there is at least one path starting at i and ending at j whose length is smaller than a fixed value $\delta$. Vertices that are not connected in ${\mathcal{G}}_{\mathrm{NL}}$ will be connected in ${\mathcal{G}}_{\mathrm{C}}$ if there is an admissible path in ${\mathcal{G}}_{\mathrm{NL}}$ joining them. In this way, we encode in ${\mathcal{G}}_{\mathrm{C}}$ connections between brain regions over long distances. In other words, we connect two vertices if the corresponding brain regions are joined by axonal paths with length determined by a $\delta$ parameter that we choose. The weights of the connections depend on how many and how long the axonal fibers in the connections are.

3. The models

By choosing appropriate graphs among the ones introduced in the previous section, we set up different models to describe the dynamics of $A\beta$ and $\tau$ proteins. The choice of the graphs is guided by the main biological features characterizing the two proteins and their dynamics.

Concerning $A\beta$, it is well known that monomeric $A\beta$ peptides originate physiologically from the cleavage of the transmembrane protein APP (amyloid precursor protein), and are considered not toxic [18]. Then, for reason not yet clarified, an imbalance between production and clearence of the protein can occur, giving rise to a process of agglomerations. This leads to the formation of toxic amyloid fibrils [19,20], often referred to as olygomers. These olygomers eventually accumulate in insoluble agglomerates known as senile plaques, nowadays considered not toxic [19]. To describe this process we adopt a so-called compartmental model, that is we consider $A\beta$ proteins assuming only three types of conformations: monomers, olygomers (which include all soluble conformations that are not monomers) and plaques (insoluble conformations). Moreover, we want to take into account that $A\beta$ protein diffuses on short distances and its dynamics (production, aggregation, clearence) is fast [21]. The molar concentration of a protein on the vertices V of a graph will be denoted by an N-dimensional vector, where $N=\left|V\right|$ is the number of vertices of the graph. Accordingly, the vector valued functions ${\mathbf{u}}_j:\mathbb{R}\to {\mathbb{R}}^N$ describe the molar concentration of $A\beta$ monomers ($j=1$), olygomers ($j=2$) and plaques ($j=3$).

We model $A\beta$ dynamics with the following equations [11]:

	(5a)
	(5b)
	(5c)

where $.\ast$ denotes the element-wise product between vectors. System (5) is endowed with initial conditions at $t=0$ as follows (at the initial time the brain is healthy, therefore only $A\beta$ monomers are present):

$$\left\{\begin{array}{c}{\mathbf{u}}_1\left(0\right)=C_{{\mathbf{u}}_1}\mathbf{1}\hfill \\ {\mathbf{u}}_2\left(0\right)=\mathbf{0}\hfill \\ {\mathbf{u}}_3\left(0\right)=\mathbf{0}\hfill \end{array}\right.$$

(6)

where $\mathbf{0}$ and $\mathbf{1}$ denote the N-dimensional vectors whose components are all equal to zero and one, respectively. Since t denotes a slow time variable, the $ϵ$ in front of the equations takes into account the fact that the processes described are fast [21]. The first term in equations (5a) and (5b) models diffusion of $A\beta$ monomers and olygomers along the network. Since this protein spreads along short distances, we choose the intrinsic proximity connectome to model the network. Indeed on the intrinsic proximity connectome ${\mathcal{G}}_{\mathrm{P}}$ two vertices i and j are connected only if the mean length of the axonal fibers connecting the corresponding regions of the brain is sufficiently short; in addition the shorter the fibre length, the greater the weight of the corresponding edge. Accordingly, here diffusion is driven by the graph Laplacian $L_{\mathrm{P}}$. In the first equation, the term $C_{{\mathbf{u}}_1}$ represents a source term (gain) due to the physiological production of $A\beta$ monomers, the term - $\alpha {\mathbf{u}}_1\left(t\right).\ast {\sum }_{j=1}^3{\mathbf{u}}_j\left(t\right)$ a loss due to the aggregtion of monomers with other $A\beta$ proteins and the term $-{\sigma }_1{\mathbf{u}}_1\left(t\right)$ a loss due to clearence phenomena. In the second equation the term ${\displaystyle \frac{\alpha }{2}}{\mathbf{u}}_1\left(t\right).\ast {\mathbf{u}}_1\left(t\right)$ is a gain in olygomers due to the aggregation of two monomers, while the next two term are losses analogous to those in the first equation. Finally, in the third equation the diffusion term is absent since plaques are insoluble (they are too heavy to diffuse), there is a gain term due to aggregation (for a detailed explanation of this term and in general on the use of Smoluchowki’s equations in this context see [6]) and a loss term due to clearence. We stress that in the equations governing the dynamics of $A\beta$ proteins there is no coupling with the $\tau$ protein. Indeed, such a coupling will be present in the equation governing the dynamics of $\tau$.

Protein $\tau$ is a physiological microtubule-associated protein: its main function is to stabilise microtubules and to regulate axonal transport. In brains affected by AD it has been found misfolded $\tau$ in form of aggregates called neurofribrillary tangles, which together with $A\beta$ plaques are an hallmark of AD [22]. Misfolded $\tau$ is toxic for the neuron: indeed it compromises stabilization, transport and in general interferes with neuronal functions. As already mentioned in the Introduction, in recent years there has been a consensus in the scientific community to consider a synergistic effect of $A\beta$ and $\tau$ when it comes to AD [4,23,24]. More precisely, the “trigger and bullet" hypothesis identifies in toxic $A\beta$ peptides the trigger for $\tau$ misfolding: the bullet is the misfolded $\tau$ which in turn causes neural damage until eventually neuron’s death [25,26,27]. Concerning the spreading mechanism of misfolded $\tau$, it has been proposed that neuronal damage spreads in the neuronal net through a neuron-to-neuron prion-like propagation mechanism [28,29,30,31]. Based on the considerations made so far, the dynamics of misfolded $\tau$, whose concentration is given by the vector valued function $\mathbf{w}:\mathbb{R}\to {\mathbb{R}}^N$, is governed by the following equation [11]:

$${\displaystyle \frac{d\mathbf{w}\left(t\right)}{dt}}={\gamma }_{3}K\left[\mathbf{w}\right]+C_{\mathbf{w}}·{\left({\mathbf{u}}_2\left(t\right)-U_{\mathbf{w}}\right)}^{+}+s_{\mathbf{w}}\left(t\right)-{\sigma }_4\mathbf{w}\left(t\right)$$

(7)

endowed with the initial condition (at the initial time the brain is healthy, therefore there is no misfolded $\tau$):

$$\mathbf{w}\left(0\right)=\mathbf{0}.$$

(8)

Here we see the coupling between $A\beta$ and $\tau$. Indeed the second term $C_{\mathbf{w}}·{\left({\mathbf{u}}_2\left(t\right)-U_{\mathbf{w}}\right)}^{+}$ in the equation rules the interaction of the two proteins: if ${\mathbf{u}}_2$, the toxic form of $A\beta$, is above a given threshold, then it triggers the misfolding of $\tau$. The term $s_{\mathbf{w}}\left(t\right)$ represents a source of misfolded $\tau$, tipically located in the enthorinal region of the brain [2], while the last term is the loss due to clearence.

The main purpose of this work is to investigate different possible forms for the operator K, appearing as first term in (7), modelling the spreding of $\tau$. As already mentioned, the spreading of such protein can take place over long distances [32], possibly following a prion-like type of process [30], on a slow time scale. Several models have appeared in literature concerning various possible mechanisms of $\tau$ spreading on networks [11,33,34,35,36,37,38,39,40]. Here we consider four different possibilities (corresponding to four different mathematical models) for the operator K, and we make a comparison between the results obtained for each model and clinical data. In this way we try to determine if there is one form for K better describing the dynamics we are interested in reproducing. The full system of equations we are going to study is therefore given by:

	(9a)
	(9b)
	(9c)
	(9d)

3.1. Model 1: diffusion of $\tau$ along the structural connectome ${\mathcal{G}}_{\mathrm{NL}}$

Following previous works in literature [6,11,33,38], a quite standard choice to model the spreading of $\tau$ protein along edges of a connectome consists in choosing the Laplacian $L_{\mathrm{NL}}$ associated to the structural connectome ${\mathcal{G}}_{\mathrm{NL}}$. In this way spreading corresponds to diffusion and we have

$$K\left[\mathbf{w}\right]=-L_{\mathrm{NL}}\mathbf{w}.$$

(10)

3.2. Model 2: diffusion of $\tau$ along the cumulative connectome ${\mathcal{G}}_c$

In the cumulative connectome, two vertices are connected if there exists at least one admissible path between them. Such a path can be formed by several consecutive edges, each with its own weight contributing to the final weight. The corresponding Laplacian $L_{\mathrm{C}}$ therefore encodes the structure of the brain with connections spanning long distances. For this reason, we choose

$$K\left[\mathbf{w}\right]=-L_{\mathrm{C}}\mathbf{w}.$$

(11)

We stress again that all distances mentioned in the paper are to be understood not in a geometric but in an intrinsic sense, i.e. referring to fibre lengths, that is edges, in the graph.

3.3. Model 3: spreading of $\tau$ via convolution on ${\mathcal{G}}_{\mathrm{L}}$

A common way of dealing with dynamics occurring over long distances is to make use of integral operators, tipically convolutions with an appropriate kernel taking into account distances. With notation similar to (4), let ${\mathcal{S}}_{ij}^{\tilde{\delta }}$ be the subset of ${\mathcal{S}}_{ij}$ containing all the paths in ${\mathcal{G}}_{\mathrm{L}}$ from node i to node j with length less than $\tilde{\delta }$:

$${\mathcal{S}}_{ij}^{\tilde{\delta }}=\{p_{ij}^k\in {\mathcal{S}}_{ij}|l_{\mathrm{L}}\left(p_{ij}^k\right)\le \tilde{\delta }\}$$

(12)

where $\tilde{\delta }\in {\mathbb{R}}^{+}$ is a fixed parameter. We define the signal $k_i$ on each node i of the graph ${\mathcal{G}}_{\mathrm{L}}$ as follows:

$${\displaystyle k_i=\sum _{p_{ij}^k\in {\mathcal{S}}_{ij}^{\tilde{\delta }}}{l}_{\mathrm{L}}\left(p_{ij}^k\right).}$$

(13)

Intuitively, $k_i$ accounts for all the paths between i and each vertex j in ${\mathcal{G}}_{\mathrm{L}}$ that is not too far from i and plays the role of the convolution kernel. The convolution product on ${\mathcal{G}}_{\mathrm{L}}$ between $\mathbf{k}$ and $\mathbf{w}$ is defined by using the graph Laplacian eigenvectors [41]. Let ${\left\{{\mathbf{v}}_{\ell }\right\}}_{\ell =1,\dots ,N}$ be a complete set of orthonormal eigenvectors of the Laplacian $L_{\mathrm{L}}$ of ${\mathcal{G}}_{\mathrm{L}}$ with eigenvalues ${\left\{{\lambda }_{\ell }\right\}}_{\ell =1,\dots ,N}$ ordered in nondecreasing order and let $U\in {\mathbb{R}}^{N\times N}$ be the matrix whose columns correspond to the eigenvectors ${\mathbf{v}}_{\ell }$, $\ell =1,\dots ,N$. The graph Fourier transform $\widehat{\mathbf{k}}$ of the signal $\mathbf{k}$ on ${\mathcal{G}}_{\mathrm{L}}$ is defined as the expansion of $\mathbf{k}$ in terms of the eigenvectors of $L_{\mathrm{L}}$, i.e.:

$$\widehat{\mathbf{k}}=U^{\ast }\mathbf{k}$$

(14)

where ^∗ denotes the conjugate transpose of a matrix. The inverse graph Fourier transform is given by

$$\mathbf{k}=U\widehat{\mathbf{k}}.$$

(15)

The convolution product on ${\mathcal{G}}_{\mathrm{L}}$ between $\mathbf{k}$ and $\mathbf{w}$ is the vector $\mathbf{k}\ast \mathbf{w}$ such that

$$\mathbf{k}\ast \mathbf{w}=U(\widehat{\mathbf{k}}.\ast \widehat{\mathbf{w}})$$

(16)

where $\widehat{\mathbf{w}}$ is the graph Fourier transform of $\mathbf{w}$ defined as in (14). We choose $K\left[\mathbf{w}\right]$ as the convolution product of $\mathbf{w}$ with $\mathbf{k}$:

$$K\left[\mathbf{w}\right]=\mathbf{k}\ast \mathbf{w}.$$

(17)

3.4. Model 4: spreading of $\tau$ via convolution on ${\mathcal{G}}_{\mathrm{NL}}$

In [36,37], the authors introduce a graph convolution operator as a nonlocal model for the conversion from an healthy protein to a toxic one. Inspired by these works, we use their nonlocal operator to model the prion-like spreading of misfolded $\tau$ protein. Let $M\in {\mathbb{R}}^{N\times N}$ be the matrix whose elements $m_{i,j}$, $i,j=1,\dots ,N$, are the length of the shortest path from vertex i to vertex j in the graph ${\mathcal{G}}_{\mathrm{NL}}$. If two vertices are not connected, $m_{i,j}=0$. Moreover, let $\tilde{M}\in {\mathbb{R}}^{N\times N}$ be the matrix with elements

$${\tilde{m}}_{i,j}=e^{-{\eta }^2{m}_{i,j}^2},i,j=1,\dots ,N,$$

(18)

where $\eta$ is a positive parameter. Finally, let $C\in {\mathbb{R}}^{N\times N}$ be the matrix defined as

$$c_{i,j}={\displaystyle \frac{{\tilde{m}}_{i,j}}{{\sum }_{k=1}^N{\tilde{m}}_{i,k}}},i,j=1,\dots ,N.$$

(19)

Then, the operator $K\left[\mathbf{w}\right]$ is chosen as

$$K\left[\mathbf{w}\right]=C\mathbf{w}.$$

(20)

4. Methods

4.1. Data Availability

Data used in the preparation of this article were obtained from the ADNI Initiative (ADNI) database (https://adni.loni.usc.edu/). As described on the ADNI website, ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial MRI, positron emission tomography (PET), other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of mild cognitive impairment (MCI) and early Alzheimer’s disease. For up-to-date information, see www.adni-info.org.

4.2. Subjects

The study included a total of 261 participants: 238 Cognitive Normal (CN) subjects (mean age $\pm st.dev.=72.7y\pm 8.7y$; $M:F=86:152$), and 23 AD (mean age $\pm st.dev.=73.1y\pm 9.6y$; $M:F=11:12$) from the Alzheimer’s Disease Neuroimaging Initiative (ADNI). All data were downloaded from ADNI3. To quantitatively measure the concentration of $\tau$ in the brain regions, we considered subjects undergone a PET imaging with the radiotracer [¹⁸F]-AV1451 and a corresponding Accelerated sagittal T1-weighted anatomical image acquired using a 3D magnetization prepared rapid acquisition gradient echo (MPRAGE) sequence from Siemens 3T MRI scanners to minimize effects of inter-scanner variability. PET and MRI imaging was chosen to be acquired at no more than 3 months of temporal distance one from the other.

4.3. PET acquisition and processing

[¹⁸F]-AV1451 imaging acquisition and pre-processing in ADNI3 consists of the baseline [¹⁸F]-AV1451 PET scan (first scan acquired) for each subject, injection of 370 MBq (10 mCi) of tracer and 30 min dynamic scan consisting of six 5-min frames 75 min post-injection. Pre-processing consists of co-registration of separated frames to minimize patient motion effect; average of the 5 frames to obtain a single 30 min PET scan; re-orientering to a standardized space and normalization of intensity and finally the application of a filter specific for the type of scanner.

4.4. MRI structural imaging acquisition and processing

The MRI structural images chosen within 3 months to the date of [¹⁸F]-AV151 PET, were pre-processed, segmented and parcellated with FreeSurfer 6.0 (http://surfer.nmr.mgh.harvard.edu/), in order to subdivide the brain volumes into a set of 83 anatomical cortical and sub-cortical region of interest (ROIs), separated between left and right hemisphere (with the addition of the brain stem) and belonging to six main networks within the brain: the frontal, the parietal, the occipital, the temporal and the limbic lobes, plus the basal ganglia [42,43]. The subdivision just described was chosen also based on the prior hypothesis that the limbic system is specially affected by the disease. A complete list of the regions, the abbreviation of their names and the network they belong to can be seen in the Appendix ??. Nodes of the brain graphs were put in correspondence with the 83 ROIs, considering the anatomical information of each vertex, to compare $\tau$ concentrations relative to the 83 ROIs in the PET images and the concentrations estimated by the mathematical models.

4.5. PET and MRI co-registration and quantification of $\tau$

To calculate the differences in [¹⁸F]-AV1451 uptake between AD and CN we co-registered the PET images of each subject to the corresponding MRI structural image using PETSurfer (https://surfer.nmr.mgh.harvard.edu/fswiki/PetSurfer). Once the PET images are aligned in the same space of the structural image, since the ROIs parcellization is in the latter space, we can calculate the $\tau$ concentration in each of the 83 ROIs. To have an absolute quantification of each subject we normalized the concentration values to the cerebellum value, since it is considered a reference region for [¹⁸F]-AV1451 metabolism.

4.6. Statistical analysis

Firstly, the normality and homogeneity of variances of the distributions of $\tau$ concentrations for each group (AD and CN) were tested using the Shapiro-Wilks and Levene’s tests. Then the distribution of [¹⁸F]-AV1451 PET uptake between AD and CN was calculated by a Mann-Whitney-Wilcoxon Test. The null hypothesis was rejected when p-values were below 0.05. In these cases, the post hoc test, Tukey’s HSD test, was employed. Post hoc tests control the family-wise error rate by the Benjamini- Hochberg step-up procedure. The ROIs for which the null hypothesis was rejected were considered significantly different between AD and CN. Significant ROIs were ordered from the most significantly different between AD and CN to the less significant until the non-significant ROIs and the corresponding value of [¹⁸F]-AV1451 concentration was considered. In Figure 1 we depict the 83 cortical and subcortical nodes of the brain with colors corresponding to different $\tau$ concentration.

• For the statistical analysis, we used Python libraries.

5. Results

This section presents the results of the numerical simulations of the four mathematical models corresponding to the different choices for the operator K governing $\tau$ spreading. All simulations were conducted using Matlab R2021a on an Intel Core i5 processor with 2.50GHZ and Windows operating system. The codes used for the current experiments can be made available upon request to the authors.

5.1. Experimental setting

The differential system (9) has been numerically solved by using the Matlab function ode45 implementing the Dormand-Prince 5(4) Runge-Kutta method with a variable time step for efficient computation [44]. Since misfolded $\tau$ protein has been shown to originate in the entorhinal region of the brain, the source term $s_{\mathbf{w}}\left(t\right)$ has been chosen non-null and equal to one only on those vertices belonging to the entorhinal region, i.e. for all $t>0$ we have:

$${\left(s_{\mathbf{w}}\left(t\right)\right)}_j=\left\{\begin{array}{cc}j\hfill & j\in \mathrm{entorhinal}\mathrm{region}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(21)

where ${\left(s_{\mathbf{w}}\left(t\right)\right)}_j$ denotes the jth component of $s_{\mathbf{w}}\left(t\right)$. The seeding region for misfolded $\tau$ is shown in Figure 2.

The differential system has been integrated from the initial time $t_0=0$ to the final time $t_f=150$. The final time was chosen to correspond to the stabilisation of the solutions. In the numerical simulations of the four models, we have fixed the values of all parameters of system (9), except ${\gamma }_3$, as reported in Table 1. These values have been chosen, on the basis of different tests, as the ones giving better results compared with clinical data.

Moreover we have fixed for the four models, with the same criterion as above, the parameters reported in Table 2. The procedure used for identifying, for each model, a good ${\gamma }_3$ value will be described in the following subsection.

5.2. Identification of the parameter ${\gamma }_3$

The statistical analysis described in Section 4.6 shows that the distribution of ¹⁸F- AV1451 uptake is not normal for all ROIs. Subsequently, the non-parametric Mann-Whitney-Wilcoxon Test with a Bonferroni correction for multiple comparisons showed that 29 ROIs significantly differ between AD and CN. Table 3 shows the six ROIs which are mostly different (p-value< ${10}^{-5}$) and their respective avaraged $\tau$ concentrations.

The most significant ROIs for the accumulation of $\tau$ in AD with respect to CN are, in decreasing order, the inferior temporal region, the fusiform region, the middle temporal region, the lateral orbitofrontal region, the amygdala and the lingual region. These six regions belong to three different networks (see Appendix ??): the inferior temporal and middle temporal regions belong to the temporal network, the fusiform and the lingual regions are part of the occipital network, and the lateral orbitofrontal region and the amygdala are included in the limbic network. Since the brain is organized in several functional discrete networks, as introduced in Section 4, it is convenient to refer to the networks to which these regions belong rather than to the single ROIs. For this reason, we average the values of $\tau$ protein in the significant regions belonging to the same network. Let $w_t^{(\ast )}$ be the average of the $\tau$ values obtained from clinical data in the inferior temporal and middle temporal regions. Analogously, let $w_o^{(\ast )}$ be the mean $\tau$ value in the fusiform and the lingual regions and let $w_{\mathrm{L}}^{(\ast )}$ be the mean $\tau$ value in the later orbital and the amygdala regions.

To evaluate the results obtained from four models, we also consider the value of $\tau$ protein in the sensorimotor network that is known not to be characterized by pathological accumulation of $\tau$ in AD. A good model for AD progression should, in fact, describe the evolution of $\tau$ both in the brain areas affected by the disease and in those that are not influenced. Therefore, let $w_s^{(\ast )}$ be the average $\tau$ value in the paracentral, postcentral, precentral and superior frontal regions forming the sensorimotor network. The $\tau$ mean values obtained from clinical data constitute our benchmark; they are ordered in decreasing order as follows

$$w_t^{(\ast )}>w_l^{(\ast )}>w_o^{(\ast )}>w_s^{(\ast )}$$

resulting in the string

$$s^{(\ast )}=“tlo{s}^{″}$$

which therefore represents the decreasing order of $\tau$ values in clinical data. Brain regions with a higher concentration of misfolded $\tau$ are more deteriorated; therefore, $s^{(\ast )}$ represents a clinical deterioration pattern of a brain with AD. We stress that this deterioration pattern considers only the brain regions that are significative for AD, according to our statistical analysis. Mathematical models should reproduce in the best possible way this deterioration pattern to be reliable tools for the analysis of AD progression. To evaluate the ability of the models to reproduce the deterioration pattern, we consider for each model j, $j=1,\dots ,4$, the values of $\tau$ protein at the final time $t_f$. We average the values on the vertices belonging to the inferior temporal and middle temporal regions to get $w_t^{\left(j\right)}$, $j=1,\dots ,4$. Proceeding similarly, we obtain $w_l^{\left(j\right)}$, $w_o^{\left(j\right)}$ and $w_s^{\left(j\right)}$ as mean $\tau$ values on the vertices belonging to the later orbital and the amygdala regions, the fusiform and the lingual regions and the paracentral, postcentral, precentral and superior frontal regions, respectively. For each model j, we define the string $s^{\left(j\right)}$ formed using the letters t, l, o, s sorted in decreasing order of $\tau$ values $w_t^{\left(j\right)}$, $w_l^{\left(j\right)}$, $w_o^{\left(j\right)}$ and $w_s^{\left(j\right)}$ . We use the Hamming Distance (HD) between the strings $s^{(\ast )}$ and $s^{\left(j\right)}$, $j=1,\dots ,4$ to measure the distance between the clinical and the models’ deterioration patterns. We recall that the Hamming distance between two strings is the number of positions at which the two strings differ [45]. The value of the parameter ${\gamma }_3$ has been experimentally determined by trial and error to identify, for each model, the best outcome in terms of Hamming distance. Since, for each model, intervals of values exist for ${\gamma }_3$ giving the same optimal Hamming distance value, we identify the best ${\gamma }_3$ value in such intervals as follows. We define the vector

$${\mathbf{w}}^{(\ast )}={({\tilde{w}}_{\mathrm{T}}^{(\ast )},{\tilde{w}}_{\mathrm{L}}^{(\ast )},{\tilde{w}}_{\mathrm{O}}^{(\ast )},{\tilde{w}}_{\mathrm{S}}^{(\ast )})}^T$$

whose components, that are normalized $\tau$ values ranging from 0 to 1, are defined as follows:

$${\tilde{w}}_{\mathrm{T}}^{(\ast )}={\displaystyle \frac{w_{\mathrm{T}}^{(\ast )}}{w_{\mathrm{T}}^{(\ast )}}},{\tilde{w}}_{\mathrm{L}}^{(\ast )}={\displaystyle \frac{w_{\mathrm{L}}^{(\ast )}}{w_{\mathrm{T}}^{(\ast )}}},{\tilde{w}}_{\mathrm{O}}^{(\ast )}={\displaystyle \frac{w_{\mathrm{O}}^{(\ast )}}{w_{\mathrm{T}}^{(\ast )}}},{\tilde{w}}_{\mathrm{S}}^{(\ast )}={\displaystyle \frac{w_{\mathrm{S}}^{(\ast )}}{w_{\mathrm{T}}^{(\ast )}}}.$$

Then, we define the vectors

$${\mathbf{w}}^{\left(j\right)}={({\tilde{w}}_{\mathrm{T}}^{\left(j\right)},{\tilde{w}}_{\mathrm{L}}^{\left(j\right)},{\tilde{w}}_{\mathrm{O}}^{\left(j\right)},{\tilde{w}}_{\mathrm{S}}^{\left(j\right)})}^T,$$

where

$${\tilde{w}}_{\mathrm{T}}^{\left(j\right)}={\displaystyle \frac{w_{\mathrm{T}}^{\left(j\right)}}{M^{\left(j\right)}}},{\tilde{w}}_{\mathrm{L}}^{\left(j\right)}={\displaystyle \frac{w_{\mathrm{L}}^{\left(j\right)}}{M^{\left(j\right)}}},{\tilde{w}}_{\mathrm{O}}^{\left(j\right)}={\displaystyle \frac{w_{\mathrm{O}}^{\left(j\right)}}{M^{\left(j\right)}}},{\tilde{w}}_{\mathrm{S}}^{\left(j\right)}={\displaystyle \frac{w_{\mathrm{S}}^{\left(j\right)}}{M^{\left(j\right)}}},$$

and $M^{\left(j\right)}=max\{w_{\mathrm{T}}^{\left(j\right)},w_{\mathrm{L}}^{\left(j\right)},w_{\mathrm{O}}^{\left(j\right)},w_{\mathrm{S}}^{\left(j\right)}\}$ for $j=1,\dots ,4$. As a measure of the distance between ${\mathbf{w}}^{(\ast )}$ and ${\mathbf{w}}^{\left(j\right)}$ we can use the Root Mean Square Error (RMSE) and the Mean Absolute Error (MAE). We recall that, given a computed vector $\mathbf{w}\in {\mathbb{R}}^4$ of $\tau$ values, these error measures are respectively defined as

$$\begin{array}{cc}\hfill \mathrm{RMSE}& =\sqrt{{\displaystyle \frac{\parallel {\mathbf{w}}^{(\ast )}{-\mathbf{w}\parallel }_2^2}{4}}},\hfill \\ \hfill \mathrm{MAE}& ={\displaystyle \frac{\parallel {\mathbf{w}}^{(\ast )}{-\mathbf{w}\parallel }_1}{4}}\hfill \end{array}$$

where ${\parallel ·\parallel }_2$ and ${\parallel ·\parallel }_2$ denote the $L_2$ and $L_1$ norms. Among all the ${\gamma }_3$ values corresponding to the same optimal HD value, we heuristically determine the value which achieves the best RMSE reduction. This leads to four values for ${\gamma }_3$, one for each model, shown in Table 4. Observe that the identified optimal ${\gamma }_3$ value is the same for Model 1 and Model 2.

5.3. Numerical results

Table 5 shows the clinical and computed deterioration patterns; Table 6 reports the corresponding HD, RMSE and MAE values.

In Figure 3, we show the deterioration pattern in the significative brain regions from clinical data and from the four different models.

From the results shown above, it is clear that the only model capable of reproducing the clinical deterioration pattern is Model 2. A diffusion operator along the cumulative connectome seems to be the most appropriate operator for modelling the spreading of $\tau$ protein in a brain with AD. Model 1 uses a diffusion operator too, but comparing the deterioration patterns of Model 1 and Model 2, it is evident that the cumulative connectome, encoding information on brain regions connected along long distances, is a key tool. Even if Model 1 does not exactly reproduce the clinical deterioration pattern, it is able to distinguish between deteriorated and non-deteriorated regions, since the letter S is in the last position in the string $w^{\left(1\right)}$, as it should be. We observe from Table 6 that models 3 and 4 give the best results in terms of RMSE and MAE reduction. However, they do not reproduce the deterioration pattern of AD and do not distinguish between deteriorated and non-deteriorated regions (the letter S is in the second position both in $w^{\left(3\right)}$ and $w^{\left(4\right)}$). For this reason, they do not seem to be good models for the spreading of misfolded $\tau$ in presence of AD.

6. Discussion

In this study we had two main objectives. We were interested in: 1) comparing different modellistic possibilities concerning the spreading of $\tau$ protein in a brain with AD, when also $A\beta$ protein is present and the synergy between the two is considered; 2) producing clinical data which can be compared with the outputs of the models in order to verify their reliability. We believe that the lack of comparison between the results of a theoretical model and clinical data could be considered a deficiency that should be rectified in order to support the usefulness of mathematical models in AD [46]. To this end, we evolved the two proteins on appropriate networks, created from medical data of human connectomes. The need for different connectomes was dictated by both the physiology of the brain and the biological characteristics of the proteins themselves. We therefore considered for the evolution of $A\beta$, which travels only on short distances, an intrinsic proximity connectome, and a standard diffusion Laplacian on it. With regard to $\tau$ spreading, which is supposed to possibly travel with a prion-like mechanism and over long distances, we considered four distinct mathematical models on as many networks: a diffusion via Laplacian on a structural connectome and on a cumulative connectome, a spreading via convolution on two further different connectomes. We emphasize that the so-called intrinsic proximity connectome and the cumulative connectome have been introduced in this work for the first time. We identified from public multimodal data, with a careful statistical analysis, six regions relevant when it comes to AD and we evaluated $\tau$ concentrations in them, obtaining a degradation pattern that was crucial for us to verify the goodness of the models.

By comparison of the results of the simulations with clinical data, we saw that only the model using the cumulative connectome was able to reproduce correctly the clinical degradation pattern, that also includes a network that is not damaged by the disease. We hypothesize that this result is due to the fact that the cumulative connectome encodes information on brain regions connected along long distances. On the basis of these findings we feel we can say that the models are not all equivalent, and that the comparison with clinical data is a diriment element to be able to assess their reliability. Of course, the results obtained here relate to the specific case of AD, and it is plausible that for other diseases the same models give different performances. As a next step we are going to apply the same paradigm of analysis, i.e. close comparison with appropriate clinical data, to the models enriched by taking into account atrophy, to be more realistic and see if something relevant happens. We believe that a highly interdisciplinary study, such as the one reported in the paper, could be a very promising direction in which cutting-edge biomedical mathematical research could be heading.