Structure Preserving Polytopal Discontinuous Galerkin Methods for the Numerical Modeling of Neurodegenerative Diseases

Many neurodegenerative diseases are connected to the spreading of misfolded prionic proteins. In this paper, we analyse the process of misfolding and spreading of both $\alpha$-synuclein and Amyloid-$\beta$, related to Parkinson's and Alzheimer's diseases, respectively. We introduce and analyze a positivity-preserving numerical method for the discretization of the Fisher-Kolmogorov equation, modelling accumulation and spreading of prionic proteins. The proposed approximation method is based on the discontinuous Galerkin method on polygonal and polyhedral grids for space discretization and on $\vartheta-$method time integration scheme. We prove the existence of the discrete solution and a convergence result where the Implicit Euler scheme is employed for time integration. We show that the proposed approach is structure-preserving, in the sense that it guaranteed that the discrete solution is non-negative, a feature that is of paramount importance in practical application. The numerical verification of our numerical model is performed both using a manufactured solution and considering wavefront propagation in two-dimensional polygonal grids. Next, we present a simulation of $\alpha$-synuclein spreading in a two-dimensional brain slice in the sagittal plane. The polygonal mesh for this simulation is agglomerated maintaining the distinction of white and grey matter, taking advantage of the flexibility of PolyDG methods in the mesh construction. Finally, we simulate the spreading of Amyloid-$\beta$ in a patient-specific setting by using a three-dimensional geometry reconstructed from magnetic resonance images and an initial condition reconstructed from positron emission tomography. Our numerical simulations confirm that the proposed method is able to capture the evolution of Parkinson's and Alzheimer's diseases.


Introduction
A neurodegenerative disease is a process that causes the progressive death or function loss of neurons.Many different pathologies belong to this group and some of them are called proteinopathies because their aetiology involves misfolding and aggregation of prions into toxic and insoluble proteins [1].Typical examples of proteins that undergo this process are α-Synuclein, related to Parkinson's disease [2], and Amyloid-β, whose aggregation is a triggering mechanism of Alzheimer's disease [3].
Recently, the mathematical modelling of prion dynamics has been studied to elucidate how the physical processes at the basis of the agglomeration and diffusion processes can be related to complex brain structure and functioning.A mathematical description of the spreading at the macroscopic level can be a useful tool in clinical practice, where the use of positron emission tomography imaging (PET) is often considered too invasive and expensive [4].Moreover, for some pathologies, like α-sinucleopathies, there not exist satisfactory chemical ligands [5] that prevent diagnostic investigations, and this computed-assisted modelling is mandatory.
Concerning the numerical modelling of neurodegeneration, the most diffused mathematical description of this phenomenon is based on the Fisher-Kolmogorov (FK) equation (also known as Fisher-KPP) [6,7].This model is a nonlinear diffusion-reaction equation that is particularly suited to describe biological species' evolution [8,9].Many different numerical methods have been proposed to compute the approximate solution of the FK equation, also in the context of brain neurodegeneration.For example, we recall Finite Element Methods (FEM) [10,11], Boundary Elements Methods (BEM) [12], and Discontinuous Galerkin (DG) methods [13].
In the context of modelling neurodegenerative disorders, the solution c of the FK problem represents the (relative) concentration of misfolded proteins, which needs to be non-negative.It can be shown that in the continuous formulation, the solution of the FK equation has two equilibrium states: c = 1 and c = 0 [14].However, due to the unstable nature of the second equilibrium, at the discrete level, it is fundamental to construct a positivity-preserving numerical method to avoid numerical instabilities that lead to unphysical (negative) numerical solutions [13].For this reason, some works analyze the construction of suitable positivity-preserving methods both within the context of finite differences [15] and DG [16] methods.The latter work uses a change of variable based on the exponential transformation to ensure positivity and entropy preservation at the discrete level.
Starting from the high-order idea of [16] -limited to simplicial meshes -in this work we present and analyse a DG formulation on polygonal/polyhedral grids (PolyDG).The proposed approach presents several advantages and novelties: (i) The flexibility in the construction of the mesh, based on mesh agglomeration [17].This plays an important role, especially because of the complexity of the geometrical domain of the application at hand, i.e., the human brain; (ii) The freedom in the choice of discretization parameters, like the polynomial degree, which might locally change, from element to element [18].In the context of brain neurodegeneration, where the geometrical complexity of the domain is an issue, the use of elementwise approximation orders allows us to reduce the computational cost, without affecting the correctness of wavefront propagation; (iii) The use of higher-order time integration.Indeed, to the timescales of the brain neurodegeneration process (that typically means over decades), the use of low-order time integration methods is not convenient to catch the wave propagation correctly.For this reason, we adopt a second-order time integration strategy; (iv) Finally, we consider spatially varying and discontinuous physical parameters, which are fundamental to correctly describe the axonal diffusion of prionic proteins [10,19].
From the point of view of the analysis, we extend the proof of the existence of the numerical solution provided in [16] for the implicit Euler method, to the generic ϑ-method.The proof of existence is based on the use of the Leray-Shauder fixed point theorem and relies on the coercivity and continuity of the diffusion term.Even though the convergence of the fully discrete numerical solution to the analytical one is not theoretically proved, it is numerically demonstrated in the case ϑ = 0.5 (Crank-Nicolson (CN) scheme), with application to brain neurodegenerative diseases, and it is shown that it outperforms first-order advancing schemes.
Concerning the application to the modelling of neurodegenerative disorders, the typical solution of the FK model is a wavefront propagating inside the brain geometry.For this reason, we analyze also the capabilities of our method in approximating wavefronts, providing also a comparison with the DG method proposed in [13], which is proven to suffer possible instabilities due to the fact it does not preserve the positivity when low order polynomial degrees are employed.
The paper is organized as follows.In Section 2, we introduce the FK mathematical model and discuss its application to neurodegeneration.In Section 3, we introduce the PolyDG space discretization and the time discretization using the ϑ-method.Moreover, we show the coercivity and continuity of the variational forms in order to prove the existence of the discrete solution, and we discuss the extension of the convergence results of the fully discrete formulation.In Section 4, we present some convergence tests with a known exact solution and we discuss the accuracy of the proposed scheme in approximating travelling waves in a two-dimensional setting, making a comparison with the DG method of [13].Section 5 is dedicated to the application of the proposed method to α-Synuclein spreading in Parkinson's disease in a two-dimensional framework, employing agglomerated polygonal meshes, and Amyloid-β in Alzheimer's disease in a three-dimensional patient-specific geometry, with initial conditions reconstructed from PET images.Finally, in Section 6, we draw some conclusions and discuss future developments.

The mathematical model
In this section, we present the FK equation to describe the reaction and diffusion of misfolded proteins.Given the final time T > 0, the problem depends on the time t ∈ (0, T ] and space x ∈ Ω ⊂ R d (d = 2, 3) variables.The unknown is the relative concentration of the misfolded protein c = c(x, t), taking values in the interval [0, 1].A detailed derivation of the model can be found in [10].The problem in its strong formulation reads as follows: Find c = c(x, t) such that: where α = α(x) is the reaction parameter, representing the local conversion rate of the proteins from the healthy to the misfolded state, modelling also on the clearance mechanisms [20,21], and Due to the physical meaning of the solution c, we aim to construct a positivity-preserving numerical scheme.Following [16], we apply the exponential transformation c = e λ , where λ = λ(x, t) becomes the new unknown of the problem.As a result, we obtain the following strong formulation of the problem: Find λ = λ(x, t) such that: The homogeneous Neumann boundary condition in problem (1) reflects a homogeneous Neumann boundary condition also in problem (2).Concerning the initial condition and the Dirichlet boundary term we impose that c 0 = e λ0 and c D = e λD , respectively.We make the following assumption on the data regularity.
Assumption 1 (Data's regularity).We assume the following regularity on the data appearing in (1):

Numerical discretization
This section presents the discretization of the continuous problem (2), which is based on the polygonal discontinuous Galerkin method for the space discretization and the ϑ−method for the time advancement.

Discrete setting and preliminary estimates
Let T h be a polytopic mesh partition of the domain Ω, being the collection of disjoint polygonal/polyhedral elements K.For each element K ∈ T h , |K| denotes the Hausdorff measure of the element, and h K denotes its diameter.We set h = max K∈T h h K .Given two neighboring elements K 1 , K 2 ∈ T h , their interface is defined as the intersection of their (d − 1)−dimensional facets.In the case of d = 2, the interface is a collection of line segments and the set of all of them is denoted with F h .In the case d = 3, the interface can be a generic polygon; for this reason, we introduce a decomposition of the polygon in planar triangles collected in the set F h .Finally, we decompose F h into the union of interior faces (F I h ) and boundary faces (F B h ), i.e.F h = F I h ∪ F B h .Moreover, we assume that F B h can be further split according to the corresponding boundary condition: Assumption 2 (Mesh Regularity [22]).The mesh sequence {T h } h satisfies the following properties: 1. Shape Regularity: ∀K ∈ T h it holds : 3. Submesh Condition: There exists a shape-regular, conforming, matching simplicial submesh T h such that: • The family { T h } h is shape and contact regular.
Remark 1.We remark that most of the analysis is valid also under milder assumptions on the mesh [23]; however in this work, we need to refer to the ones in Assumption 2. The technical point is the validity of (6) that holds under mesh assumptions of Assumption 2.3.However, we notice that from the numerical results of Sections 4 and 5, the assumption seems not to be needed in practice.
Concerning the space discretization, we introduce the following discontinuous finite element spaces with an elementwise variable polynomial degree: is the space of polynomials of total degree p K ≥ 1 over a mesh element K. Concerning the physical data, we assume D ∈ W DG h,p and α ∈ W DG h,p .We introduce the following trace operators [24].Let F ∈ F I h be a face shared by the elements K ± and let n ± be the unit normal vector on face F pointing exterior to K ± , respectively.Then, for sufficiently regular scalar-valued functions v and vector-valued functions q, we define: The superscripts ± denote the traces of the functions on F taken within the interior to K ± .In an analogous way, on the face F ∈ F D h associated with the cell K ∈ T h with n outward unit normal on ∂Ω, we define: Let us introduce the following broken Sobolev spaces for an integer r ≥ 1: . We define the following penalization function η : where ζ(p, h, D) is defined as We point out that in Equation ( 3), we are considering both the harmonic average operator {•} H and the arithmetic average operator {•} A on F ∈ F I h and η 0 is a parameter at our disposal (to be chosen large enough to ensure stability).Moreover, we are defining Finally, we can define the following DG-norm: Remark 2. The choice of using this combination of harmonic and arithmetic averages is fundamental to obtaining the coercivity and continuity bounds of Propositions 1 and 2 below.
Finally, we recall the result of inverse trace inequality [25]:

PolyDG semi-discrete formulation
To construct the semi-discrete formulation, we first introduce the interior penalty DG discretization of the nonlinear diffusion term A : where ∇ h • is the elementwise gradient [26] and η is defined as in Equation ( 3).The semi-discrete PolyDG formulation reads as follows: For any t ∈ (0, T ], find λ h (t) ∈ W DG h,p such that: where λ 0h ∈ W DG h,p is a suitable approximation of λ 0 ∈ W .We next show some preliminary estimates that will be needed in the forthcoming well-posedness and convergence analysis.
Proposition 1 (Coercivity of A ).The form A , defined in Equation (7), satisfies for all v ∈ W DG h,p : under the assumption on the penalty parameter value η 0 ≥ 16C 2 I D 0 , where d 0 and D 0 are defined as in Assumption 1, and C I is the inverse trace inequality constant of relation (6).
Proof.Taking u = v = w in Equation ( 7), we have: By treating each term separately, we obtain for (I) the following estimate: Then we control the term (II) by means of the Young's inequality: where β F > 0 is a parameter we define as follows: In ( 13) d 0 and D 0 are defined as in Assumption 1 and C I is the inverse trace inequality constant of relation (6).
Let us recall the following relation: Then, by applying the inverse trace inequality and relation ( 14) we obtain: Inserting the above estimates in Equation ( 9), we obtain: For F ∈ F I h , the second integral on the rhs of Equation ( 15) is positive provided that: The same bound can be obtained on F ∈ F D h .By taking η 0 ≥ 16C 2 I D 0 the positivity is guaranteed, and by exploiting the following relation max{(e v ) + , (e v ) − } max e ∥v∥ L ∞ (K + ) , e ∥v∥ L ∞ (K − ) ≥ 1 we obtain: where ζ has been defined in (4).
Proposition 2 (Continuity of A ).The form A , defined in Equation (7), satisfies for all u, v ∈ W DG h,p : with µ := max 1, , where d 0 and D 0 are defined as in Assumption 1, C I is the inverse trace inequality constant in (6), and η 0 is the penalty constant introduced in 4.
Proof.From Equation ( 7) we obtain: By treating each term separately, we obtain for (I) the following estimate using the regularity assumption on D in Assumption 1: Then, we control the term (II) by means of the Young's inequality: The bound on term (III) follows thanks to the Young's inequality: where γ F > 0 is defined as follows: Let us recall the following relation: Then, by applying the inverse trace inequality in relation (6) and relation (23) we obtain: From the above estimates it follows: Finally, we estimate the term (IV) by applying the inverse trace inequality in relation (6) and relation (23): Finally, putting together all the previous bounds, we obtain: and the proof is complete.

Fully discrete formulation
To discretize Equation ( 8) in time, we consider the ϑ−method scheme.We remark that due to the nonlinear nature of the strong formulation with the change of variable, we need a nonlinear solver, and therefore using an implicit scheme for time integration does not affect the computational cost.In this section, we consider homogeneous Dirichlet conditions λ D = 0 for simplicity in the calculations.However, the results can be extended to the nonhomogeneous case with proper regularity assumptions on λ D .Let {t ℓ } Nt ℓ=0 be the uniform partition of the time interval [0, T ] into N t intervals with length dt = T Nt , namely, 0 = t 0 < t 1 < ... < t Nt = T and t ℓ = ℓT Nt for ℓ = 0, ..., N t .Let us introduce a parameter ε > 0.Then, the fully discrete formulation of problem (8) reads: given the initial condition λ 0 = λ 0 , find λ k+1 h for k = 0, ..., N t − 1, such that: The introduction of the additional regularizing terms proportional to the parameter ε > 0 is fundamental to prove the existence of the solution via the Leray-Schauder fixed-point theorem [16].However, from a numerical point of view, the presence of a ε > 0 is not really needed and it can be chosen equal to 0 in the simulations (see Section 4).We next prove that formulation (25) admits a solution.
Proposition 3 (Existence of a solution).Let ε > 0. Given λ k h ∈ W DG h,p , then the fully discrete formulation in Equation (25) admits a solution λ k h ∈ W DG h,p , provided that Assumptions 1 and 2 hold and the penalty constant η 0 defined as in (4) is chosen sufficiently large.
Proof.The proof is based on the application of the Leray-Schauder theorem.For clarity, we subdivide the proof into 3 steps.
Step 1: Definition of the operator Φ.First of all, let us introduce the fixed point operator Φ : h,p such that Φ(w, σ) = v with v ∈ W DG h,p being the unique solution of the linear problem: Step 2: Compactness of Φ. Φ is well defined by the Lax-Milgram lemma, thanks to the coercivity and continuity on W DG h,p of the left-hand side of (26) and to the continuity of the right-hand side of (26).Finally, we observe that Φ(w, 0) = 0. Due to the finite dimension of the space W DG h,p , these properties are enough to prove also the compactness of the operator.
Step 3: Uniform bound for all the fixed points.To prove the property of uniform bound we take v ∈ W DG h,p and σ ∈ [0, 1] such that v = Φ(v, σ).First of all, let us notice that we can bound the right-hand side of (26) by using the coercivity of A and the existence of a constant ).Indeed, there holds Then, by introducing the function s(x) = x(log(x) − 1) + 1 ≥ 0 and exploiting its convexity we obtain: Thus, using also the fact that −s(e v ) ≤ 0 and relation (27) we obtain: Using Equation ( 17) and the Young's inequality with suitable coefficients ϵ 1 and ϵ 2 , we get: By applying the Leray-Schauder theorem [14] we derive the existence of a solution for problem (25), and the proof is complete.

Convergence of the discrete solution
In this section, we prove the convergence of the solution to the PolyDG fully discrete formulation in Equation ( 25) with ϑ = 1 (Implicit Euler method) to the solution of the continuous problem.An additional assumption we make in this proof is the forcing-free model f = 0.The result follows by extending the convergence theorem proved in [16] to the case pf polytopal/polyhedral meshes and high-order approximations.
Let us start introducing the notion of entropy S : [0, T ] → R of the system [27], namely To prove the convergence of the numerical solution, we need to show that the discrete entropy S k h = Ω (e λ k h (λ k h − 1)+1) decays as k → ∞ [16], and that the DG norm (see equation ( 5)) of the discrete solution is uniformly bounded.
Remark 3. The analysis in this section is performed only for the case ϑ = 1.The treatment of the general case ϑ ∈ [0, 1] is not straightforward, due to the presence of the components from the previous timestep that cannot be easily treated and prevent to recover the decay of the discrete entropy.Nevertheless, as it will be demonstrated in the numerical result sections, the scheme exhibits optimal convergence rates for any ϑ ∈ [0, 1].The extension of the analysis to the case ϑ ̸ = 1 is under investigation and will be the subject of future research.Proposition 4. Let Assumptions 1 and 2 hold and let η 0 defined as in Equation (4) be chosen sufficiently large.Let λ k+1 h be the solution to problem (25) with ϑ = 1, ε > 0 and homogeneous forcing term f = 0. Then: where S 0 h is the initial discrete entropy.Proof.Let us consider the problem (25) with ϑ = 1 and φ h = λ k+1 h : .
The proof follows the same steps as in [16] and it makes use of Propositions 1 and 4, as well as of the extensions of variational inequalities valid for polygonal/polyhedral meshes.

Numerical results: verification
In this section, we aim at verifying the accuracy of the method presented in section 3.

Test case 1: convergence analysis in two dimensions
For the numerical tests in this section, we use Lymph library [28] to solve the FK equation (d = 2).We define the domain Ω = (0, 1) 2 , which we discretize by means of a polygonal mesh obtained by using PolyMesher [29].Concerning the time discretization, we use a timestep ∆t = 10 −6 and the final time T = 2 × 10 −5 .We consider the following manufactured exact solution: λ(x, y, t) = log (cos(πx) cos(πy) + 2)e −t . ( We fix the physical parameters as follows: D = I and α = 0.1.The forcing term and the Dirichlet boundary conditions are derived accordingly.To solve the resulting nonlinear system we adopt a Newton method with tolerance equal to ϵ = 10 −10 . In Figure 1, we report the computed errors in both the DG and L 2 norms at the final time.We have performed the convergence test keeping fixed the polynomial order of the space approximation p K = p = 1, ..., 6 ∀K ∈ T h and using different mesh refinements (N el = 30, 100, 300, 1000).It can be observed that the slope of error decrease is equal to the polynomial degree p for the DG-norm and equal to p + 1 for the L 2 -norm.
In Figure 2a, we report the computed errors with respect to the timestep ∆t, using both Crank-Nicolson (ϑ = 0.5) and the Implicit Euler (ϑ = 1) schemes.The space discretization is computed on a mesh of N el = 1000 elements and with polynomial degree p = 6.As expected, the use of the Crank-Nicolson method leads to a second-order convergence whereas the error decays with a first-order rate if the implicit Euler scheme is employed.We remark that the case ϑ = 1 is fully covered by our theoretical analysis, whereas the proof of convergence for ϑ = 1/2 is under investigation.
A convergence analysis with respect to the polynomial order p is also performed on a coarse mesh of 30 elements and with a time integration based on the Crank-Nicolson scheme with timestep ∆t = 10 −6 .The results are reported in Figure 2b, where can observe exponential convergence can be observed.

Test case 2: Travelling waves in two dimensions
In this section, we exployt the positivity-preserving PolyDG formulation to simulate a traveling-wave solution of the FK equation in 2D, with the aim of comparing the formulation we propose in this article, with the (non positivity-preserving) scheme introduces in [13].The manufactured solution is of the form: By substituting it in Equation ( 1) with f = 0, we obtain the following equivalent system of ordinary differential equations: where we have used the assumption of isotropic diffusion tensor D = dI.In particular, we fix d = 10 −3 , α = 1 and penalty parameter η 0 = 1.Concerning the wave's parameters we take the speed v = 0.1 and the initial data ψ(0) = 1 and χ(0) = −10 −2 .We consider a rectangle Ω = (0, 5) × (0, 1) as domain.We present the results of two simulations, with different final times T = 5 and T = 10 and timestep ∆t = 10 −2 .Concerning the nonlinear Newton solver, we fix a tolerance ϵ = 10 −6 .
In Tables 1 and 2, we report the computed errors in the L 2 and DG norms computed at the final times T = 5 and T = 10, respectively.In particular, we compare the results obtained by using our positivity-preserving method (25) and the DG method proposed in [13], with a semi-implicit treatment of the nonlinearity and a penalty parameter η = 10.We can observe that, also using low order polynomials (p = 1), our method is able to correctly represent the wave propagation front and leads to smaller errors (one order of magnitude).On the contrary, the method in [13] fails to correctly simulate the wavefront because it does not preserve the positivity of the solution and the equilibrium c = 0 is unstable.
Moreover, from the results of Table 1 and Table 2, we can observe for p = 4 and T = 10 that the proposed positivity-preserving scheme does not lead to a reduction of the error compared with the results obtained for p = 3.Indeed, we can observe in Figure 3 that for p = 3 we have the formation of some small oscillations around the equilibrium c = 1.This is probably due to Newton's iterations that might be badly conditioned for large values of polynomial degrees.The effect of this problem cannot be observed in the method of [13], but in this case, the positivity of the solution cannot be guaranteed.Table 2: Comparison of the computed errors in the DG-norm based on employing the proposed positivity-preserving scheme and the DG method of [13], for different polynomial degrees, different mesh sizes, and different final times.

Numerical results: brain applications
In this section, we present the numerical results obtained in two different test cases: a two-dimensional simulation of a sagittal section of a brain and a three-dimensional simulation of brain geometries reconstructed from Magnetic Resonance Images (MRI).
In the prions' spreading applications, the diffusion tensor is typically modelled as the superimposition of an extracellular diffusion effect with magnitude d ext and an axonal diffusion with magnitude d axn [10]; for this reason, in this section, we assume that D has the following structure: where n = n(x) is the axonal fibres direction in the point x ∈ Ω and d ext , d axn ≥ 0. The axonal direction is derived from Diffusion Weighted Imaging (DWI) and represents the principal orientation of the connections between the neurons (axons).Most of the spreading of the prions seems to happen through the axons [10], however, due to the brain structure, this is true only in white matter, while in grey matter, the diffusion can be considered to be isotropic.
In order to construct the axonal component of the diffusion tensor D, we derive the diffusion tensor from DWI medical images by using Freesurfer and Nibabel [30].The principal eigenvector n of the tensor is then computed elementwise to find the diffusion tensor in Equation (37).We refer to [31] for more details on the reconstruction of

Test case 3: spreading of α-Synuclein in a two-dimensional brain section
In this section, we address a numerical simulation of the spreading of α-Synuclein in Parkinson's disease on a polygonal agglomerated grid of a sagittal 2D brain section.The geometry is segmented from a structural MRI of a brain from the OASIS-3 database [32] by means of Freesurfer [33].The construction of the final mesh of a slice of the brain is performed by using VMTK [34].The resulting triangular mesh is composed of 43 402 triangles, and each element of the mesh is labelled to be in white or grey matter, according to the MRI segmentation, as in Figure 4a.However, the generality of the PolyDG method allows us to use mesh elements of any shape and the use of a smaller number of elements allows saving computational cost.For this reason, by using ParMETIS [35], we agglomerate the initial triangular mesh into a polygonal mesh of 534 elements, as shown in Figure 4b.In particular, the agglomeration procedure is performed in a segregated way for the white and the grey matter, in this way we are sure to correctly describe both the domain boundary and the interface between grey/white matters.Finally, in Figure 4c, we report the axonal directions computed in the white matter starting from DWI.
Concerning the physical parameters, we fix the reaction coefficient α = 0.45/year in grey matter, and α = 0.9/year in white matter [36].Moreover, we impose a constant isotropic diffusion d ext = 8 mm 2 /year, and axonal diffusion which is 10 times faster than the isotropic one in the white matter (d axn = 80 mm 2 /year) and is negligible in the grey matter (d axn = 0 mm 2 /year) [36].In this simulation, we fix ∆t = 0.01 years and p = 1, moreover the penalty parameter η 0 = 1.
The simulation of α-Synuclein diffusion in Parkinson's disease starts from an initial condition, with concentration located in the dorsal motor nucleus [37].In Figure 5, we report both the initial condition and the computed solution at different times t = 0, 5, 10, 15, 20, 25 years.First of all, it can be observed that the directions of protein propagations are coherent with the medical literature [38].Indeed, the activation of brain regions follows the Braak staging theory [37].Moreover, we can notice that the heterogeneity of the reaction parameters causes an earlier activation of the white matter in general, which is clearly visible in the frontal cortex at time t = 20.By making a comparison with the literature results of [13], we have that the reduced reactivity and diffusion inside grey matter   causes a slowing of the disease progression times, starting with the same initial condition and an agglomerated mesh with comparable refinement level.
In Figure 6a, we report the average concentration of misfolded protein e λ h (t) inside the brain with respect to the time t.Moreover, we compute the average concentrations in white and grey matter separately.As we can observe, in the first years, the increase in the concentration is almost equivalent for the two regions, after 14 years we have a clear distinction.In particular, the higher reactivity and diffusion of the white matter tissue causes a faster increase in the concentration.Moreover, we compute the activation time of the pathology as: where χ is the indicator function and c crit = 0.95 is the critical value of α-Synuclein concentration.We report the computed activation time in Figure 6b.From a pathological perspective, high concentrations of α-Synuclein alter the electric signal transport.The indicator (38) measures the time after which a region of the brain can be affected    by pathological electric stimuli.The result is qualitatively similar to the literature results [10,13].Comparing the result with respect to [13], we can notice a longer activation time, due to the reduced reactivity and diffusion in grey matter, introduced in this work.

Test case 4: spreading of Amyloid-β in a three-dimensional brain
In this section, we present a numerical simulation of the spreading of the Amyloid-β on a three-dimensional domain, reconstructed starting from an MRI taken from OASIS-3 database [32].The medical images are associated with an 83 years old patient, which is diagnosed to be affected by Alzheimer's disease at the moment of the acquisition.The geometry is segmented by means of Freesurfer [33] and then is used to construct a mesh grid of 323'014 tetrahedral elements, using SVMTK library [31].The resulting mesh is reported in Figure 7a.The problem is solved with the use of a FEniCS code [39] (version 2019).
Concerning the parameters of the model, in this test case, for simplicity, we do not make any distinction between white and grey matters, choosing α = 0.9/year, d ext = 8 mm 3 /year, and d axn = 80 mm 3 /year [36,13].
To set up the initial condition for the FK problem in a patient-specific setting, we estimate the function λ 0 (x) of Amyloid-β protein at the initial time t = 0. To do that, we project the clinical data derived from PET images with Pittsburgh compound B (PET-PiB) [4].The PET-PiB adopts a radioligand, which identifies the presence of Amyloid-β plaques inside the brain parenchyma (for the specifics about the acquisition techniques of the image used in this work we refer to [32]).We report the result of the initial concentration rescaled between 0 and 1 and projected on the mesh grid in Figure 7b.In particular, we can observe the presence of large damaged regions (c ≃ 1) in the brainstem and in the thalamus.
Starting from pathology in an advanced state, we set up a simulation with a final time T = 2 years and a timestep ∆t = 0.01 years.Concerning the space discretization we adopt the DG method for p = 1.The nonlinear solver for the resulting system is based on the relaxed Newton method with absolute tolerance equal to 10 −10 and relaxation parameter ω = 0.75.
The results are reported in Figure 8 for different times (t = 0, 0.5, 1, 1.5, 2 years).The solution is visualized on many slices inside the brain geometry on the three different planes: horizontal, coronal and sagittal.The results show a propagation of the Amyloid-β concentration inside the parenchyma, following the typical paths of the pathology [4].In particular, we can observe a late activation of the cerebellum in the slice along the coronal plane of Figure 8 (middle line).This is coherent with Braak's stages of Alzheimer's pathology, which show the presence of Amyloid-β accumulation only in the last stages of the pathological development [40].Moreover, coherently to the clinical stage of the pathology we are simulating (due to the presence of evident symptoms from the patient's documentation), we can find a generalised misfolding after a few years from the PET acquisition and this is also coherent with what we expected in the disease evolution [4].

Conclusions
In this work, we have proposed a positivity-preserving DG method on polygonal and polyhedral grids for the solution of the FK model.The main applicative motivation is the modelling of neurodegeneration caused by the spreading of prionic proteins, such as α-synuclein protein in Parkinson's disease and amyloid-β in Alzheimer's disease.We have analyzed the existence of the discrete solution by means of the Leray-Schauder theorem and we have discussed the convergence of the numerical scheme.
Numerical tests have been presented both in two and three dimensions.In particular, we have analyzed the convergence in space both with respect to the mesh size and the polynomial order of the method on polygonal grids.Then, we have discussed the convergence in time, by making a comparison between implicit Euler and Cranck-Nicolson schemes.Finally, we have performed a numerical simulation to test the capabilities of the proposed formulation to approximate propagating wavefronts in two dimensions.In this test, we have compared the proposed positivity-preserving method with the polyDG method introduced in [13], highlighting the advantages and disadvantages of both formulations.Finally, we have presented two applications of the proposed scheme in the framework of neurodegenerative diseases.In particuarl, we have performed a simulation of α-synuclein spreading on a slice of a real brain in the sagittal plane, constructing a polygonal agglomerated mesh that preserves the quality of both domain boundaries and the interface between white matter and grey matter.Moreover, starting from initial amyloid-β concentrations derived from PET images, we have simulated the spreading of amyloid-β in a three-dimensional brain in a patientspecific Alzheimer's disease setting.The results obtained in both patient-specific settings are coherent with the clinical literature, showing that the proposed approach is a valuable instrument that can be employed for patientspecific computed-assisted simulations of the evolution of Parkinson's and Alzheimer's neurodegenerative disorders.
A possible future development of this work consists in extending the convergence analysis to the general ϑmethod, by proving a discrete entropy decay.Another possibility can be the use of PET images at different times of the disease to calibrate the physical parameters of the Fisher-Kolmogorov model, for example by means of inverse uncertainty quantification methods.
is the diffusion tensor, denoting the spreading of misfolded protein.Finally, f = f (x, t) is the forcing term modelling the external addition of mass.Concerning the boundary conditions, we impose null flux at the boundary Γ N of the domain, while c D fixes the value of concentration on a part of the boundary Γ D , where {Γ D , Γ N } form a partition of ∂Ω, namely, Γ D ∪ Γ N = ∂Ω, Γ D ∩ Γ N = ∅, and |Γ D | > 0.
where |F | is the Hausdorff measure of the face F .

Figure 2 :
Figure 2: Test case 1: computed errors and convergence rates in convergence concerning timestep (left) and polynomial degree (right) with N el = 30.

Figure 3 :
Figure 3: Test case 2:Comparison of the numerical solutions computed both using the proposed positivitypreserving DG scheme and the DG method presented in[13], with the exact solution at time T = 10.

Figure 4 :
Figure 4: Brain section with the specification of the white-gray matter subdivision (a), the agglomerated mesh (b), and the axonal directions (c).

Figure 5 :
Figure 5: Patterns of α-synuclein concentration at different stages of the pathology.
Mean value of the concentration inside the brain, white, and grey matter (b) Activation time on the brain section

Figure 6 :
Figure 6: Test case 3: mean value of the concentration (a) and activation time (b).
(a) Three-dimensional tetrahedral brain mesh (b) Initial condition from PET images on brain slices in horizontal plane

Figure 7 :
Figure 7: Test case 4: three-dimensional brain mesh (a) and projected initial condition from PET image (b).
where F D h and F N h are the boundary faces contained in Γ D and Γ N , respectively.The last assumption implies that any F ∈ F B h is contained in either Γ D or Γ N .

Table 1 :
[13]arison of the computed errors in the DG-norm based on employing the proposed positivity-preserving scheme and the DG method of[13], for different polynomial degrees, different mesh sizes, and different final times.