Polytopal discontinuous Galerkin discretization of brain multiphysics flow dynamics

A comprehensive mathematical model of the multiphysics flow of blood and Cerebrospinal Fluid (CSF) in the brain can be expressed as the coupling of a poromechanics system and Stokes' equations: the first describes fluids filtration through the cerebral tissue and the tissue's elastic response, while the latter models the flow of the CSF in the brain ventricles. This model describes the functioning of the brain's waste clearance mechanism, which has been recently discovered to play an essential role in the progress of neurodegenerative diseases. To model the interactions between different scales in the porous medium, we propose a physically consistent coupling between Multi-compartment Poroelasticity (MPE) equations and Stokes' equations. In this work, we introduce a numerical scheme for the discretization of such coupled MPE-Stokes system, employing a high-order discontinuous Galerkin method on polytopal grids to efficiently account for the geometric complexity of the domain. We analyze the stability and convergence of the space semidiscretized formulation, we prove a-priori error estimates, and we present a temporal discretization based on a combination of Newmark's $\beta$-method for the elastic wave equation and the $\theta$-method for the other equations of the model. Numerical simulations carried out on test cases with manufactured solutions validate the theoretical error estimates. We also present numerical results on a two-dimensional slice of a patient-specific brain geometry reconstructed from diagnostic images, to test in practice the advantages of the proposed approach.


Introduction
In the brain, multiple fluid components play different roles: the blood supplies nutrients and oxygen and removes carbon dioxide, while the Cerebrospinal Fluid (CSF) and the glymphatic system [61,51,22] have the primary function of clearing the waste produced by brain activity.Also, it has been recently shown that waste clearance mechanisms play a major role in the evolution of neurodegenerative diseases [15,55,50,24].These physical systems are strongly interconnected with one another and with the cerebral matter: for example, a large amount of the CSF is generated in the choroid plexus thanks to the high concentration of blood capillary vessels in it [25]; also, an auxiliary function of the CSF is to protect the brain from impact against the skull and to compensate blood pulsatility in terms of flow rate and pressure inside the braincase.For this reason, the modeling of the fluid dynamics in the brain requires a multiphysics perspective, able to capture these interactions and their mutual interplay.
From the modeling viewpoint, the brain can be described as a porous material (white and gray matter) with fluids both filtrating through it and flowing in hollow regions (brain ventricles).In terms of mathematical modeling, this can be represented as the coupling between a poromechanics system -described, e.g., by Darcy's or Biot's equations [74] -and Stokes/Navier-Stokes' equations for the flow in the ventricles [58,48].Numerical methods for such Fluid-Poroelastic Structure Interaction (FPSI) systems can be found in the literature: the most studied is the Biot-Stokes' system [16,3,2,23], but also more complex models have been investigated, e.g.considering a multilayer structure of the porous medium [21,28], a non-linear constitutive model for either the structural or fluid part of the system [40], or Multi-compartment Poroelasticity (MPE) models coupled with Stokes' equations [41,18].This latter framework is the best suited to model brain poromechanics and CSF flow, since it can simultaneously represent cerebral tissue deformation, blood vessel networks at different scales, CSF flow, and waste clearance.Moreover, the dynamic formulation of MPE equations can describe the inertial forces associated with blood pressure pulsatility during a heartbeat, which affects vascular and tissue deformations and waste clearance [71,33,39,35].The analysis of conforming Finite Element methods for FPSI problems has been addressed in several of the abovementioned works, with a discussion on the inf-sup requirements on the different components of the system (poroelastic, fluid, and coupling conditions).However, aiming at an accurate representation of the poromechanics and flow with low dispersion and dissipation errors, high-order discretization methods provide a better choice [66,49].Polytopal Discontinuous Galerkin (PolyDG) methods fit in this framework, and they can also naturally account for complex geometries such as vessel networks and brain folds, thanks to their flexibility in terms of local refinement and agglomeration, hanging nodes treatment, and generality of mesh element shapes [19,31,11,30,65,29].In addition, discontinuous Galerkin schemes provide a general framework to embed physically-consistent interface conditions directly in the weak form.So far, these methods have been mostly employed to solve porous media and poroelasticity problems for geophysical applications [59,12,10] and for fluid dynamics and fluid-structure interaction problems [34,42,72,7,73,75,14].
In this work, we introduce a high-order PolyDG method to spatially discretize a coupled model encompassing dynamic Multiple-Network Poroelastic (MPE) equations for poromechanics [35] and Stokes' equations [8,14], with physically-consistent coupling conditions inspired from the mass and stress balance at the interface between the brain tissue and the CSF.Moreover, we take full advantage of the DG framework to embed, directly in the physically consistent formulation, the coupling conditions at the interface between the poroelastic and fluid regions.We analyze the well-posedness and convergence of the semidiscrete numerical method, and we employ a combination of Newmark's β-method and the θ-method for time discretization.
The paper is organized as follows.In Section 2, we introduce the multiphysics problem, both in strong and weak form, and discuss the coupling conditions.Section 3 describes the PolyDG space discretization; stability and convergence properties are analyzed in Section 4.Then, time discretization is introduced Section 5. Verification tests are discussed in Section 6, corroborating the theoretical results of Section 4. Section 7 demonstrates the capabilities of the method on a two-dimensional brain section considering physiological settings.

Mathematical model
We introduce the mathematical model consisting of the coupling between a Multiple-Network Poroelasticity system and Stokes' system: the former, introduced in [35], accounts for the poromechanics of the brain tissue and its interaction with different fluids compartments flowing in its pores, while the latter describes the flow of the CSF in the brain ventricles.To ease the presentation, we consider a simplified configuration whose two-dimensional representation is depicted in Fig. 1.The poroelastic medium occupies a portion Ω el ⊂ R d (d = 2, 3) of the domain, while the CSF flows according to Stokes' equations in the remaining portion Ω f ⊂ R d .The interaction between the two systems occurs at the interface Σ = Ω el ∩ Ω f , which we suppose to be a (piecewise) smooth (d − 1)−manifold.The source of CSF comes from an exchange of mass with the poroelastic medium and it exits from the domain at Γout, while the rest of the domain boundary is a solid wall Γw.We denote by Ω the interior of Ω el ∪ Ω f .Given a final observation time T > 0, we introduce the CSF velocity u : and the network pressures pj : Ω el × [0, T ] → R, j ∈ J, where J is a given set of labels denoting the different fluid network compartments [35].In particular, for the application at hand, we consider J = {A, C, V, E}, where A, C, V correspond to the arterial, capillary, and venous blood compartments, respectively, while E denotes the extracellular CSF permeating the brain tissue.(1) with corresponding physiological values.
The coupled problem reads as follows: where the linear elastic and fluid (viscous) stress tensors are defined as are supposed sufficiently regular.The parameters of the models are explained in Table 1, with their typical physiological values extracted from [57,35,32].The initial and boundary conditions of problem (1) are defined as with suitable definition of the data function p out : Γout × (0, T ] → R, that represents the external normal stress at the outlet, and of the initial conditions d0 : On the interface Σ, we introduce the following coupling conditions, based on physiological considerations: Condition (3a) expresses the balance of total normal stress.Due to the blood-brain barrier [1,38], we assume that mass exchange between the poroelastic domain and the CSF only occurs through compartment E, as expressed by (3b)-(3c).Consistently, the normal stress of the CSF fluid is balanced by the pressure of compartment E, as in (3d), while we assume the tangential stress on the fluid to be negligible (cf.(3e)).Similar assumptions were made in [32], although we do not make use of the Beavers-Joseph-Saffman condition [64].
Aiming at solving problem (1) with the Finite Element method, we introduce its weak formulation.For the sake of generality, let Γ D,d , ΓD,u, ΓD,P j , with j ∈ J, denote the portions of ∂Ω where Dirichlet boundary conditions on d, u, pj are imposed, respectively.Then, we introduce the following functional spaces: where H 1 (Ω) denotes the classical Sobolev space of order 1 over L 2 (Ω).For the problem at hand, all Dirichlet boundary conditions on ∂Ω are homogenous.We denote by (•, •)Ω the L 2 -product over Ω and we define the following forms and functionals over the spaces introduced above: Remark 2.1 (Derivation of the interface form J) The interface form J : → R introduced above naturally arises during the derivation of the weak form of problem (1).We test (1a)-(1c) against functions w ∈ W and qj ∈ Qj, with j ∈ J, over Ω el , and (1d) against v ∈ V over Ω f .Then, integrating by parts and summing all the contributions yield the following boundary terms on the interface: Using the interface conditions (3a)-(3b) and then (3c)-(3d)-(3e), we can rewrite (4) as follows: where we also used that a ⊗ b : Denoting by L 2 (0, T ; H), H 1 (0, T ; H) the time-dependent Bochner spaces associated to a Sobolev space H, and setting the weak formulation of problem (1) reads as follows:  for all (w, {qj}j∈J , v, q) ∈ D × P × V × Q, and d(0) = d0, ∂td(0) = ḋ0, u(0) = u0, pj(0) = pj0 ∀j ∈ J.
3 Semidiscrete formulation based on a polytopal discontinuous Galerkin method In this section, we introduce a space discretization of problem (6) based on discontinuous Finite Element methods on polytopal grids.

Notation
Let T h,el , T h,f be polytopal meshes discretizing the domains Ω el , Ω f , respectively.We define as faces of an element K ∈ T h,el ∪ T h,f the (d − 1)-dimensional entities corresponding to the intersection of ∂K with either the boundary of a neighboring element or the domain boundary ∂Ω: • for d = 2, the faces are always straight line segments; • for d = 3, the faces are generic polygons.We assume that each face can be decomposed into triangles.
With this definition, we denote by F el , F f the sets of element faces corresponding to each physical domain.We partition them into internal faces for the elastic displacement, the pressure of the j-th compartment, and the fluid velocity, respectively), and interface faces F Σ ⊂ Σ.In the latter, we assume that the meshes T h,el , T h,f are aligned with Σ, namely that there is no gap or overlap between them, although hanging nodes are permitted.
We introduce the symmetric outer product v ⊙ n = 1 2 (v ⊗ n + n ⊗ v) and, for regular enough scalar-, vector-and tensor-valued functions q, v, τ , we define the following average and jump operators: where n + , n − are defined as in Fig. 2 -left.
where n is the unit normal vector pointing outward to the element K to which the face F belongs.
• On a face F ∈ F Σ shared by two elements K el ∈ T h,el and K f ∈ T h,f : where n el , n f are defined as in Fig. 2 -right.

PolyDG semidiscrete problem
For a given integer m ≥ 1, we introduce the following piecewise polynomial spaces: Moreover, we denote by H s (T el ), H s (T f ) the broken Sobolev spaces of order s over the mesh of the poroelastic and fluid domains, namely H s (T⋆) = {q ∈ L 2 (Ω⋆) : q|K ∈ H s (K) ∀K ∈ T⋆} for ⋆ = el, f.We then introduce the following forms, for d, w where ∇ h , ε h , div h denote the element-wise gradient, symmetric gradient, and divergence operators, respectively, and the stress tensors σ el , σ f are implicitly defined in terms of these piecewise operators.The parameters η, ζj, γv, γp appearing in these forms are defined as follows [14,35]: where {h}H denotes the harmonic average on K ± (with {h}H = hK on Dirichlet faces) 2 are the L 2 -norms of the symmetric second-order tensors appearing in the elasticity and Darcy equations, for each K ∈ T el , and η, ζj ∀j ∈ J, γ v , γ p are penalty constants to be chosen large enough.

The form J(•, •, •) is the piecewise discontinuous version of the interface form
where we have used the identity a•b = a ⊙ b : I, ∀a, b ∈ R d and then the definition of the average and jump operators introduced in Section 3.1.
Finally, the semidiscrete formulation reads as follows: Problem ( 10) is supplemented with suitable initial conditions d h (0), ḋh (0), {p j,h (0)}j∈J , u h (0) that are projections of the initial data introduced in (1) onto the corresponding DG spaces.The bilinear forms appearing in (10) are defined as

Algebraic formulation
We introduce suitable sets of basis functions such that span{φ i el } Denoting by uppercase letters the d.o.f.vectors corresponding to the problem unknowns, the (formal) algebraic form of (10) is the following: Given D0, Ḋ0, U0, Pj0, j ∈ J, find D, U , P , Pj, j ∈ J such that The matrices and vectors of ( 12) are defined as follows (where ⋆ = el, f):

A priori analysis of the semidiscrete problem
For the analysis contained in this section, we consider a generic set J of NJ ∈ N0 compartments and we assume that all the physical parameters of the model (defined in Table 1) are piecewise constant.For r ≥ 1, we introduce the following broken norms [35,14]: We introduce the energy norms at time t ∈ (0, T ] , and set For the sake of simplicity, in the inequalities appearing in the following analysis, we neglect the dependencies on the model parameters and the polynomial degree m, and we use the notation x ≲ y to indicate that ∃C > 0 : x ≤ Cy, where C is independent of the mesh discretization parameters.Moreover, we make the following assumption on the mesh [11,30]: Assumption 4.1 For each h > 0, the two meshes T h,el , T h,f are aligned with Σ, namely there is no gap nor overlap between them (hanging nodes are allowed).Moreover, denoting by T h the union of T h,el and T h,f , we consider a sequence of meshes {T h } h>0 satisfying the regularity requirements of [35]: • {T h } h>0 is h-uniformly polytopic-regular, namely for each K ∈ T h there exists a set where hK is the diameter of K, and the cardinality of {S F K } F ⊂∂K does not depend on h > 0. • There exists a shape-regular simplicial covering T h of T h such that, for each • The mesh size satisfies a local bounded variation property: -dimensional measure and all the hidden constants independent of h and the number of faces of K1 and K2.
Under these assumptions, we can prove the following stability result: Theorem 4.1 (Stability estimate) Let Assumption 4.1 hold true and let us also assume that the penalty constants defined in (8) are chosen sufficiently large.Then, the semidiscrete solution (d h , {p j,h }j∈J , u h , p h ) of (10) satisfies the following inequality for each t ∈ (0, T ]: where, according to the initial conditions of (1), Proof.Let us fix a time t ∈ (0, T ] and consider the following test functions in ( 6): w = ∂td(t), v = u(t), q = p(t), qj = pj(t) ∀j ∈ J.With this choice of test functions, the terms J cancel out, as well as the B forms in the L terms Proceeding as in [35], we integrate ( 16) w.r.t.time and we can employ integration by parts in time in the following terms to obtain Using these identities and the definition of the bulk and broken norms introduced above, we can rewrite the left-hand side of ( 16) (integrated over time) as follows, and use Cauchy-Schwarz's and Young's inequalities on the right-hand side: We now consider continuity and coercivity results for the bilinear forms appearing in (17) with respect to the broken norms (13).These results, proven in [35,14], are reported in Appendix A (Lemma A.1) and they include an inf-sup condition for the form B f , with a constant β f,h independent of the mesh elements size.According to [14], if γv is large enough, there exists α > 0 such that α = O(β −2 f,h ) and Analogously, the coercivity results of the forms A el , Aj, Cj and the continuity of We can use these results on the left-hand side of ( 17) to obtain The definitions of || • ||EN,t and || • ||EN,0 and the fact that ρ el , cj, ρ f are constant conclude the proof.We now proceed to derive an a-priori error estimate for the error introduced by the space discretization of the problem.For any t ∈ (0, T ], let (d, {pj}j∈J , u, p) denote the weak solution of problem (6) and let (d h , {p j,h }j∈J , u h , p h ) denote the semidiscrete solution of problem (10), obtained with sufficiently large stability parameters as defined in (8).
Using the inverse trace inequality [13,11] we can state the following continuity result for the interface forms, whose proof is provided in Appendix B: Lemma 4.1 Under Assumption 4.1, the following inequalities hold: We introduce the following additional norms for non-discrete functions: Denoting by EK : H s (Ω) → H s (R d ) the Stein extension operator for a Lipschitz domain Ω defined in [69], the following interpolation result can be stated (cf.[35,14,8,6]): Lemma 4.2 Under Assumption 4.1, the following estimates hold: , where K ⊇ K, for each K ∈ T h , are shape-regular simplexes as in Assumption 4.1.
Combining the results above, we can prove the following optimal convergence result, whose proof is provided in Appendix C. Theorem 4.2 (A priori error estimate) Under the same assumptions of Theorem 4.1, if the solution of problem (6) is sufficiently regular, and the initial conditions d 0 , ḋ0 , u 0 , p 0 j , j ∈ J are sufficiently regular, the following estimate holds for each t ∈ (0, T ] and for each m ≥ 1: ||(e d , {e P j }j∈J , e u , e where e d = d − d h , e P j = pj − p j,h ∀j ∈ J, e u = u − u h , e P f = p − p h , and K ⊇ K, for each K ∈ T h , are shape-regular simplexes as in Assumption 4.1.

Fully discrete problem
We introduce a uniform partition {tn} N n=0 of the interval (0, T ], with constant timestep ∆t = tn+1 − tn, for all n = 0, . . .N − 1. Starting from the algebraic form (12) of problem (10), we discretize the elastic momentum equation (first row) with Newmark's β-method, whereas we employ the θ-method to discretize all the compartment pressure equations and the fluid problem.For Newmark's discretization, we introduce two auxiliary vector variables Z n , A n representing the expansion coefficients of the approximate elastic velocity and acceleration at time tn [35].The resulting algebraic problem has the form where and the expression of the matrices A1, A2 are reported in Appendix D.

Verification tests -convergence
In this section, we verify the theoretical error bounds of Theorem 4.2.In all tests, we consider only one pressure compartment, that is j ∈ J = {E}.The constant coefficients of the penalty parameters defined in (8) are set as Regarding time discretization, all the results were obtained with β = 0.25 and γ = 0.5 for the Newmark scheme and θ = 0.5 for the θ-method, both for Stokes' problem and the pressure compartment of the poro-elastic system.
The tests have been implemented in the 2D version of lymph [9], an open-source MATLAB library for the solution of multiphysics problems with the PolyDG method, developed at MOX, Department of Mathematics, Politecnico di Milano.

Test case 1: steady solution
In this test, we consider a manufactured steady solution of problem ( 1)-( 2)-( 3), so that we can verify the convergence of the semidiscrete formulation without spoiling the results with time discretization error.In particular, we introduce the following exact solution on the 2D domain Ω = Ω el ∪Ω f = (−1, 0)×(0, 1)∪(0, 1)×(0, 1), with interface Σ = {0}×(0, 1) depicted in Fig. 1: The source terms of (1) corresponding to the exact solution above are With these data, we performed spatial convergence tests setting αE = 0.5 and all the physical parameters ρ el , µ el , λ el , cE, µE, KE, β e E , ρ f , µ f equal to 1.In Fig. 3 we report the computed error.From the results, we can clearly observe that the error decays w.r.t.h with a rate of order m, as predicted by our theoretical estimate (19).Moreover, we can observe spectral convergence of the error w.r.t. the polynomial approximation degree m.

Test case 2: unsteady solution
We consider the rectangular domain Ω = Ω el ∪ Ω f = (−1, 0) × (0, 1) ∪ (0, 1) × (0, 1) of Section 6.1 and we introduce the following functions of time: A time-dependent exact solution of problem ( 1)-( 2)-( 3) can be manufactured combining these functions with the steady solution introduced in Section 6.1: d(x, y, t) = g el (t)d steady (x, y), pE(x, y, t) = g el (t)p steady E (x, y), u(x, y, t) = gu(t)u steady (x, y), p(x, y, t) = gp(t)p steady (x, y), (23) with suitable definitions of the source terms and boundary data.Again, we performed convergence tests setting αE = 0.5 andall other physical parameters to 1.We choose a time step ∆t = 10 −3 s and a final time T = 5 ∆t.We report the computed errors in Fig. 4 (log-log scale).As predicted by Theorem 4.2, we observe that the error in the energy norm decays at a rate proportional to h m , for any m ≥ 1.Moreover, even though not covered by our theoretical analysis, we observe exponential convergence of the error for fixed h and increasing m.Finally, we observe that the value chosen for ∆t is small enough to prevent error saturation for the sequence of meshes considered in this test.

Numerical results on a 2D slice of the brain
In this section, we demonstrate the capability of the proposed method for the solution of a realistic problem on a 2D slice of the brain and its ventricles, shown in Fig. 5.The geometry of the problem is based on structural Magnetic Resonance Images (MRI) available in the OASIS-3 database (https://oasis-brains.org)[56].By means of Freesurfer (https://surfer.nmr.mgh.harvard.edu/)[44], a three-dimensional brain geometry was segmented, and then sliced along the sagittal plane by VMTK (vmtk.org)[5].The resulting triangular 2D meshes of the cerebral tissue and of the brain ventricles surrounded by it are composed of 25847 and 3286 elements, respectively (see Fig. 5).Yet, the flexibility of the PolyDG method allows to employ meshes with elements of generic shape: by agglomeration, we can thus considerably reduce the number of mesh elements while retaining the same geometrical detail of the original triangular grid.We agglomerate the grid by means of ParMETIS (https://github.com/KarypisLab/ParMETIS)[54], obtaining the polygonal mesh shown in Fig. 5, consisting of 910 elements in T el and 101 in T f .This agglomeration is performed separately for the two physical domains, to preserve geometrical accuracy at the interface Σ.
Aiming at reproducing conditions in the physiological regime, we set the data and parameters of problem (1) as follows.We split the Dirichlet boundary Γw into Γ w,el representing the dura mater membrane surrounding the brain tissue and the boundary Γ w,f of the corpus callosum (the whole in the center of Fig. 5).In the poroelastic problem, we consider only the compartment related to the extracellular CSF, namely J = {E}, and we assume no flow  (∇pE • n el = through the dura mater Γ w,el .The source of extracellular CSF is given by gE = 2 • 10 −3 π sin(2πt)s −1 , representing the variations in the interstitial CSF due to blood pulsation [57,35,32], with a period of 1 s representing a heartbeat.No additional external forces act on the poroelastic tissue or the ventricle flow, that is f el = f f = 0.No slip conditions u = 0 are enforced on Γ w,f , while a no-stress condition p out = 0 is imposed on the obex Γout of the fourth ventricle, where the CSF would flow into the central canal of the spinal cord.Zero initial conditions are imposed on all variables and the values of the remaining physical parameters of the problem are reported in Table 1, with the slight modifications αE = 0.49, β e E = 0 m 2 • N −1 • s −1 : these values are in the physiological range of the brain function [57,35,32].For the PolyDG method we employ a polynomial degree m = 2 for all variables, while for the time advancement based on Newmark's (β = 0.25, γ = 0.5) and Crank-Nicolson's (θ = 0.5) methods we use a time step ∆t = 10 −2 s from t0 = 0 s to T = 1 s.
Selected snapshots of the computed solutions are reported in Figs. 6 and 7. Due to the choice of a zero outlet pressure p out , pressures pE, p should be interpreted as pressure differences (w.r.t. the fourth ventricle outflow) rather than absolute pressure values.We notice that the amplitude of the brain displacement d is always below 0.2 mm, thus justifying the choice of a linear elasticity modeling of the tissue, at least in the current settings.Due to such small displacements, the interstitial pressure pE at the interface and the ventricle CSF pressure p are substantially equal through the whole simulation timespan, in accordance with the interface conditions (3).The ventricle pressure p undergoes very small variations, which are better observable in Fig. 7.We can see that the CSF circulates in the third ventricle around the corpus callosum, whereas its flow through the cerebral aqueduct is always directed outwards, thus fulfilling its clearance function.This clearance occurs notwithstanding the fact that the source term gE is negative in the second half of the timespan, and consequently the interstitial pressure pE is lower than the outlet pressure p out = 0: this effect is due to the inertial terms in the problem formulation, thus highlighting their importance in the model.
All the observations above are in general agreement with measurements and computational results in the literature, particularly in the magnitude of the displacement d and in the pressure difference between the third ventricle and the outflow [25,32].Some discrepancies can be observed in the overall range of pE and the magnitude of u, but they do not exceed one order of magnitude.These discrepancies may be due to the fact that the computational domain of the simulation presented here is a 2D slice of the brain, and does not include the lateral ventricles: simulations in three dimensions, capturing the complete geometry of the brain, are envisaged to address this issue.

Conclusions
We introduced a discontinuous Galerkin method for the discretization of a multiphysics fluid-mechanics model of the brain, encompassing a Multi-compartment Poro-Elastic (MPE) modeling of the tissue, perfused by blood and CSF, and Stokes' flow of the CSF in the brain cavities.Our numerical method is particularly suitable for capturing the high geometrical complexity of the brain, thanks to the possibility of supporting high-order approximations and agglomerated grids.We proved stability and optimal error bounds for the semidiscrete formulation under standard assumptions in the PolyDG framework.Specific attention was paid to the modeling and numerical treatment of the interface conditions between the MPE and fluid domain.We implemented the method in our PolyDG library lymph [9].
We performed verification tests, validating the optimal order of convergence of our method with respect to the mesh element size and spectral convergence with respect to the polynomial degree.Then, we showed the capability of our computational model to represent the physiological clearance function of the CSF in physiological settings, on a two-dimensional sagittal slice of the brain and of its ventricles.
Several directions for further development of the present work should be addressed.First, in terms of modeling, nonlinear hyperelastic rheology could be considered to account for the extremely soft nature of the brain tissue [63,27,4].This would require more detailed fluid-structure interaction conditions on the tissue-fluid interface, which have not been fully studied in the brain ventricles, but have been widely investigated in other biological systems such as the heart and blood vessels [67,43,52,26,46].Second, to thoroughly analyze the effect of the complex brain geometry on waste clearance, the model could be applied to full three-dimensional geometries of different subjects.To this aim, suitable geometry reconstruction techniques based on MRI data should be employed, such as those developed in [68,62,70,45,47,20,17].Third, to better represent the effect of blood pulsation on CSF flow (here modeled by a homogeneous CSF generation term), the blood compartments of the MPE model should be coupled to models of the main brain arteries [53,60].Finally, the coupling of the model proposed here with a suitable description of the generation, aggregation, and transport of misfolded proteins such as amyloid-β and tau, would allow the investigation of their clearance in different physiological and pathological conditions [36,37].This would entail the development of numerical methods able to account for very different time scales, since the model presented here is defined on the scale of seconds (characteristic of blood pulsation) while prion aggregation and neurodegeneration occur over a timespan of decades.

A Continuity and coercivity in the discrete spaces
We collect the continuity and coercivity results of [35,14] in the following Lemma, which is instrumental to the proofs of stability (Theorem 4.1) and convergence (Theorem 4.2) of our numerical method (10).
Lemma A.1 Under Assumption 4.1, the forms of (7) are continuous over the discrete spaces: Moreover, provided that the penalty constants are chosen sufficiently large, the following coercivity and inf-sup inequalities hold (for any j ∈ J): where β f,h is the discrete inf-sup constant [14].
B Proof of Lemma 4.1 . For any F ∈ F Σ h , we denote by K F el and K F f the elements of T h,el and T h,f , respectively, sharing F in their boundary.Then, where, in line TR , we have used the inverse trace inequality (18) and the fact that η|K ∼ h −1 K ∀K ∈ T h,el and γp|K ∼ hK ∀K ∈ T h,f .An analogous argument (using γ ) allows to control |J(q h , w, v)|, thus concluding the proof.

C Proof of Theorem 4.2
In this section, we prove the optimal convergence estimate stated in Theorem 4.2.This result is based on the inverse trace inequality (18), on Lemma 4.1 and on the following continuity results (proved in [10,8,14]), which extend Lemma A.1 to consider also non-discrete functions: Lemma C.1 Under the same assumptions of Lemma A.1, the following inequalities hold: We are now ready to prove our optimal convergence result.h .An analogous notation is introduced for the errors e P j ∀j ∈ J, e u , e P f .We fix a time t ∈ (0, T ] and we proceed similarly to [35], namely we subtract equation (10) h , e u ) + S(e P f , e According to the splitting introduced above, we separate the terms depending only on the discrete approximation errors from those involving the interpolation errors.Then, integrating in time from 0 to t and proceeding as in the proof of Theorem 4.1, the coercivity inequalities of Lemma A.1 yield @ @ @ @ @ @ @ J(e P E h , ∂te d h , e u h ) − @ @ @ @ @ @ @ J(e where no contribution arises from the initial conditions, due to the hypotheses.Using Lemmas 4. Regarding the interface terms, we can observe that the choice of η, γv, γp made in (8) implies that η, γv, γ −1 p scale as hK in each element K ∈ T h .Thus, thanks to Lemma 4.2, we obtain ∥η and similar optimal estimates for the other interpolation errors at the interface appearing in the last two lines of (25).By Cauchy-Schwarz's and Young's inequalities, all the terms involving the discrete errors e d h , e h , e u h , e P f h on the right-hand side of (25) can be moved to the left-hand side, whereas the interpolation errors e d I , e

D Algebraic form of the fully discrete problem
The matrices A1, A2 of (20) have the following form:

3 . 5 • 10 − 3
kg • m −3 density of the solid tissue ρ f 1000 kg • m −3 density of the CSF µ el 216 Pa first Lamé parameter of the solid λ 505 Pa second Lamé parameter of the solid µ j Pa • s viscosity of the fluid in compartment j ∈

Figure 2 :
Figure 2: Polygonal elements sharing an internal face (left) or a face on the interface Σ (right).

Figure 4 :
Figure 4: Test case 2. Left: computed errors in the energy norm (14) versus h, with different polynomial degrees m (same for all variables).Right: computed errors in the energy norm (14) versus m, with h = 0.273 corresponding to N = 80 polygons.

Figure 5 :
Figure 5: Triangular mesh of the 2D slice reconstructed from MRI (top left; zoom on top right to see the triangular elements), polygonal computational mesh obtained by agglomeration (bottom left), and topographic anatomy of the fluid domain Ω f (bottom right).

Figure 6 :
Figure 6: Brain simulation of Section 7. Discrete tissue displacement d h and pressures p E,h , p h in the whole domain Ω el ∪ Ω f at different times during a heartbeat.Same scale for p E,h , p h .

Figure 7 :
Figure 7: Brain simulation of Section 7. Velocity u h and pressure p h in the brain ventricles Ω f at different times during a heartbeat.The cerebral aqueduct is indicated by the box.
Proof.(Theorem 4.2) Using the interpolators defined in Lemma 4.2, we split the error e d = e d I − e d h into an interpolation error e d I = d − dI ∈ [H m+1 (T h,el )] d and an approximation error e d h = d h − dI ∈ W DG Aj + Cjj ∀j ∈ J.
1/2 e P E I ∥ F Σ K ∥EK pE∥ H m+1 ( K) , (26)N,t that is completely analogous to(26)can be proven by resorting to Lemma 4.2, the triangle inequality ||(e d , {e P j }j∈J , e u , e P f )|| 2 EN,t ≤ ||(e d h , {e I }j∈J , e u I , e P f I h }j∈J , e u h , e P f h )|| 2 EN,t + ||(e d I , {e P j I }j∈J , e u I , e P f I )|| 2 EN,t concludes the proof.