Magazine Cover
Download article (PDF)

This article is published open access.

Waterscales: Mathematical and computational foundations for modelling cerebral fluid flow

  • Marie E. Rognes

    Simula Research Laboratory, Oslo, Norway

1 Introduction

On 18 October 2013, I was sitting on my living room sofa scrolling through the news, when one article caught my attention. The headline read something along the lines of “Sleep cleans the brain of toxins”, and described a recent original research study from the lab of neuroscientist Prof. Maiken Nedergaard [39 L. Xie, H. Kang, Q. Xu, M. J. Chen, Y. Liao, M. Thiyagarajan, J. O’Donnell, D. J. Christensen, C. Nicholson, J. J. Iliff et al., Sleep drives metabolite clearance from the adult brain. Science 342, 373–377 (2013) ]. In a series of experiments, Nedergaard and her team had injected a fluorescent dye (a so-called tracer or contrast agent) into the fluids surrounding the brain of mice, and observed that the tracer moved into (and out of) the brain many times faster in sleeping mice than in awake mice. Their findings revealed a fundamental interplay between sleep and brain clearance, but also highlighted how far from complete our understanding of the brain’s waterscape (Definition 1) is.

Definition 1. The brain’s waterscape is the circulation, flow and exchange of tissue fluid and associated solute transport through and around the brain.

Our brains are composed of very soft tissue consisting of neurons, glial cells, and interstitial space filled with interstitial fluid (ISF), penetrated by blood vessels, and surrounded by a bath of cerebrospinal fluid (CSF). Its well-being crucially relies on the transport of solutes: the influx of oxygen and nutrients, and clearance of metabolic waste. Due to the barrier layers between blood, tissue and CSF in the brain, the fluid exchange between these spaces is limited and regulated. Tracer experiments act as a proxy for quantifying this fluid flow and exchange, in lack of more direct medical techniques for measuring water movement in-vivo (Figure 1)1All clinical data presented and visualized here courtesy of P. K. Eide and G. Ringstad, Oslo University Hospital – Rikshospitalet.. In the nearly ten years since 2013, the brain’s waterscape has been attracting substantial and increasing interest from the neurocommunity; more about that to follow.

A year before the striking sleep study [39 L. Xie, H. Kang, Q. Xu, M. J. Chen, Y. Liao, M. Thiyagarajan, J. O’Donnell, D. J. Christensen, C. Nicholson, J. J. Iliff et al., Sleep drives metabolite clearance from the adult brain. Science 342, 373–377 (2013) ], Nedergaard had launched a new theory for describing the brain’s waterscape: the glymphatic system [22 J. J. Iliff, M. Wang, Y. Liao, B. A. Plogg, W. Peng et al., A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid β. Sci. Transl. Med. 4, DOI 10.1126/scitranslmed.3003748 (2012) ]. In the original glymphatic concept, CSF enters the brain via perivascular spaces surrounding brain arteries, mixes with ISF and flows in bulk and convectively through the tissue, and ISF is cleared from the brain along perivascular spaces surrounding brain veins. With my background in analysis of numerical methods for partial differential equations describing biological tissue, I was immediately intrigued by this setting. I also quickly realized that this was an underdeveloped area for applied mathematics, and that even fundamental concepts for multiscale modelling of the brain’s waterscape were missing. Fortunately, the European Research Council Mathematics Panel agreed, and funded my project proposal in the summer of 2016.

In the next few years, we very much realized that the brain’s waterscape was even less understood and more discussed that what I had originally thought, and that its mechanisms were surrounded by controversy. It seemed as though almost every piece of the glymphatic theory was disputed with different groups arguing about, e.g., the role of diffusion versus convection, the existence of any convective flow and, if existent, its directionality, the importance vs. nonimportance of aquaporin channels and glial cells, the pathways for clearance, the anatomy of the compartments involved, etc. The list could go on and on. For us, this made the need for mathematical and computational foundations even more obvious to, e.g., quantify observations, bridge between species, and test, reject or support hypotheses.

2 Quantifying brain solute transport

Transport of solutes is often described via the usual convection-diffusion-reaction equation: in a domain ΩRd\Omega\subset\mathbb{R}^{d} (d=1,2,3d=1,2,3) and over a time interval (0,T](0,T], find the concentration c=c(x,t)c=c(x,t) such that

ctdivDgradc+div(cφ)+rc=0,c_{t}-\operatorname{div}D\operatorname{grad}c+\operatorname{div}(c\varphi)+rc=0,

where the subscript tt denotes the time derivative, div\operatorname{div} and grad\operatorname{grad} are the divergence and gradient respectively, DD is the diffusion tensor, φ\varphi is a convective velocity field and rr is a kinetic reaction rate. In addition, cc may be prescribed on the boundary Ω\partial\Omega for a.e. 0<tT0<t\leqslant T, and initially c(x,0)=c0(x)c(x,0)=c_{0}(x) for xΩx\in\Omega. Numerical approximations of (1) are readily computed using finite difference and finite element methods, via, e.g., the FEniCS software [2 M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes and G. N. Wells, The FEniCS project version 1.5. Archive of Numerical Software 3, 9–23 (2015) ].

Figure 1. Contrast agent concentrations c1,c2c_{1},c_{2} in a human brain at 6 (left) and 24 (right) hours after injection in the lower spine.

2.1 Modelling the glymphatic theory via stochastic fields

Now, in our first step towards modelling human brain solute transport [7 M. Croci, V. Vinje and M. E. Rognes, Uncertainty quantification of parenchymal tracer distribution using random diffusion and convective velocity fields. Fluids Barriers CNS 16, 1–21 (2019) ], we took a just-published clinical study as a starting point [32 G. Ringstad, S. A. S. Vatnehol and P. K. Eide, Glymphatic MRI in idiopathic normal pressure hydrocephalus. Brain 140, 2691–2705 (2017) ]. Ringstad, Vatnehol and Eide [32 G. Ringstad, S. A. S. Vatnehol and P. K. Eide, Glymphatic MRI in idiopathic normal pressure hydrocephalus. Brain 140, 2691–2705 (2017) ] had injected a contrast agent into the fluid-filled spaces surrounding the spine, and then imaged its whereabouts within the brain 1–48 hours later (Figure 1). They concluded that it seems unlikely that diffusion alone explains the brain-wide distribution. We asked: “is it?” and were curious as to whether mathematical modelling could support, quantify and/or add to this statement.

Stochastic diffusion

First, in order to quantify the unlikeliness of diffusion as a sole mechanism, a key question was how to capture the uncertainty associated with the brain’s diffusivity DD as a random field. Importantly, we aimed to ensure that diffusivity samples were positive, with a literature-based expected value DgD^{*}_{g}, and appropriate variability. To this end, we represented DD as a continuous random field

D(x,ω)=0.25Dg+Df(x,ω),D^{*}(x,\omega)=0.25D^{*}_{g}+D_{f}^{*}(x,\omega),

where for each fixed xΩx\in\Omega, Df(x,ω)D^{*}_{f}(x,\omega) is a gamma-distributed random variable with a prescribed shape (k=3k=3) and scaling (θ=0.75Dg/k\theta=0.75D^{*}_{g}/k). To also enforce continuity and to effectively sample the diffusion field, we drew samples of DvD_{v}^{*} by first sampling a Matérn field (with a given smoothness and correlation length), and then mapping it onto a gamma random field via a copula (see, e.g., [7 M. Croci, V. Vinje and M. E. Rognes, Uncertainty quantification of parenchymal tracer distribution using random diffusion and convective velocity fields. Fluids Barriers CNS 16, 1–21 (2019) ] and references therein for further details). Next, we pieced together inputs from several sources: a human brain finite element mesh with nearly two million vertices and ten million cells, and a upwards-travelling contrast agent distribution on the mesh boundary based on an eye-norm estimation of the published data [32 G. Ringstad, S. A. S. Vatnehol and P. K. Eide, Glymphatic MRI in idiopathic normal pressure hydrocephalus. Brain 140, 2691–2705 (2017) ].

In short, we found that the uncertainty in the diffusion field magnitude had a substantial impact on the amount of contrast agent both in the grey and white matter. While the clinically observed distribution of contrast in the grey matter was well within the expected variability of diffusion, its distribution into the white matter was not. We could thus support the claim that diffusion was likely not sufficient as a sole mechanism, but at the same time highlight that diffusion was likely nearly sufficient. Further, we demonstrated that local variations (i.e., heterogeneity) in the diffusion field had little impact on expected values, and thus that contrast agent distribution in larger brain regions could be well approximated by average diffusion coefficients.

Stochastic convection

We next asked how we could represent the convective velocity φ=(φ1,φ2,φ3)\varphi=(\varphi_{1},\varphi_{2},\varphi_{3}) described by the glymphatic theory in a stochastic setting at the brain scale? Well,

  • φ\varphi varies more after a distance proportional to the mean distance between arterioles and venules (λ1\lambda\approx 1 mm);

  • the vasculature is random and independent of space in the sense that the presence of paraarterial or paravenous spaces are equally likely at any point in space (hence the expected value of each φi\varphi_{i} is zero);

  • φ\varphi is continuous and divφ=0\operatorname{div}\varphi=0; and

  • older experimental studies indicate an expected velocity magnitude E(φ)=vavg=0.17E(\lVert\varphi\rVert)=v_{\mathrm{avg}}=0.17 µm/s with some variability.

In turn, we defined the stochastic glymphatic circulation velocity field φ\varphi as the curl of three standard independent identically distributed Matérn fields with correlation length λ\lambda and scaled to align with the expected value and variability. Then, again using Monte Carlo simulations, we computed cc samples via (4), but now with non-zero φ\varphis.

Surprisingly (or perhaps not), we found that the expected tracer distribution with or without the glymphatic velocity were nearly identical and with little variability. Thus, on average, small fluctuations in the CSF/ISF velocity did not increase (or decrease) the influx or clearance of contrast agent in the brain on a larger scale. Indeed, further model variations, such as including a directional flow at the macroscale and allowing for local fluid influx and/or efflux, were required for the convective terms to effectively contribute to the solute transport [7 M. Croci, V. Vinje and M. E. Rognes, Uncertainty quantification of parenchymal tracer distribution using random diffusion and convective velocity fields. Fluids Barriers CNS 16, 1–21 (2019) ].

Remark 2. Despite its limitations, this uncertainty quantification study felt like a true breakthrough: we had taken advantage of an original mathematical approach and produced new physiologically relevant insights. Importantly, it was the first study by my team that had been accepted and published in Fluids and Barriers of the Central Nervous System, a domain-specific journal read by everyone and anyone interested in brain mechanics, including neurosurgeons and neuroradiologists, neuroscientists, biophysicists, bioengineers (and mathematicians).

Remark 3. Another, more mathematical, question is how to sample Matérn field efficiently, especially in the context of more advanced Monte Carlo methods over non-trivial geometries. As such, this physiological problem setting also guided us to develop new sampling and Monte Carlo algorithms [8 M. Croci, V. Vinje and M. E. Rognes, Fast uncertainty quantification of tracer distribution in the brain interstitial fluid with multilevel and quasi Monte Carlo. Int. J. Numer. Methods Biomed. Eng. 37, Paper No. e3412 (2021) ].

(a) Subject-specific mean diffusivity Dˉ\bar{D}^{*} extracted from diffusion tensor imaging.
(b) Optimal convective flow field φ\varphi (streamlines) between 6 and 24 hours post contrast injection (cf. Figure 1), estimated by high-dimensional inverse modelling [40].
Figure 2.

2.2 Identifying convective flow: An inverse problem

With these insights, we understood more about the roles and properties of diffusion and convection via CSF/ISF flow within the brain. We also understood that in order to get much further, we needed a closer collaboration with the clinicians and better access to the raw data. Luckily, Prof. P.-K. Eide and Dr. G. Ringstad at Oslo University Hospital were very willing to share their data and expertise with us and very enthusiastic about exploring how computational modelling could contribute to quantify and interpret their clinical findings.

We then found ourselves in the unprecedented position of having multi-modal magnetic resonance (MR) data, including T1- and T2-weighted structural images, T1-maps, diffusion tensor images and contrast MR images, available for several individual in different cohorts. MR data are characteristically at high spatial resolution, but only available at a sparse set of discrete time points, e.g., {1,6,24,48}\{1,6,24,48\} hours after contrast injection. These images combined allowed us to segment and construct subject-specific finite element meshes, interpolate subject-specific diffusion tensor fields (Figure 2 (a)), and map MR contrast signal intensities onto finite element concentration fields for each subject and each time point (Figure 1).

Remark 4. In fact, to handle this data set, we developed and published an open-source pipeline going from magnetic resonance brain images to finite element simulations, accompanied by an open-access introductory booklet [25 K.-A. Mardal, M. E. Rognes, T. B. Thompson and L. M. Valnes, Mathematical modeling of the human brain: From MRI to finite element simulation. Springer (2022) ].

An optical flow problem?

Now, given this data set, could we identify and quantify the convective flow of ISF/CSF within the brain? Or the task in more mathematical terms: given non-negative scalar fields c1,c2c_{1},c_{2}, identify a transport field φ\varphi that maps c1c_{1} onto c2c_{2} in a (to be defined) suitable sense.

This is a classical problem setting, known in computer vision as the optical flow problem [21 B. K. Horn and B. G. Schunck, Determining optical flow. Artif. Intell. 17, 185–203 (1981) ]. The idea is that given c1,c2 ⁣:ΩRc_{1},c_{2}\colon\Omega\to\mathbb{R} at t1,t2t_{1},t_{2} for 0t1<t20\leqslant t_{1}<t_{2}, find an optimal φ ⁣:Ω×[t1,t2]Rd\varphi\colon\Omega\times[t_{1},t_{2}]\to\mathbb{R}^{d} that minimizes an objective functional JJ,

minφJminφt1t2Ωf2dxdt+α2R2.\min_{\varphi}J\equiv\min_{\varphi}\int_{t_{1}}^{t_{2}}\int_{\Omega}f^{2}\,\mathrm{d}x\,\mathrm{d}t+\alpha^{2}R^{2}.

In the original optical flow setting, it is assumed that φ\varphi maps c1c_{1} to c2c_{2} by convective transport alone – corresponding to

f=ct+div(φc),f=c_{t}+\operatorname{div}(\varphi c),

and the initial and final conditions c(t1)=c1c(t_{1})=c_{1}, c(t2)=c2c(t_{2})=c_{2}. The choice of regularization term RR (and also ff for that matter) gives variations of the method; one option is to minimize the kinetic energy:

R2=cφφL1([t1,t2])×L1(Ω)t1t2 ⁣Ωcφφdxdt.R^{2}=\lVert c\varphi\cdot\varphi\rVert_{L^{1}([t_{1},t_{2}])\times L^{1}(\Omega)}\equiv\int_{t_{1}}^{t_{2}}\!\int_{\Omega}c\varphi\cdot\varphi\,\mathrm{d}x\,\mathrm{d}t.

This formulation yields the L2L^{2}-Monge–Kantorovich mass transfer problem, analysed across three centuries (see [3 J.-D. Benamou and Y. Brenier, A computational fluid mechanics solution to the Monge–Kantorovich mass transfer problem. Numer. Math. 84, 375–393 (2000) ] and references therein), and used by Tannenbaum, Ratner et al. [30 V. Ratner, L. Zhu, I. Kolesov, M. Nedergaard, H. Benveniste and A. Tannenbaum, Optimal-mass-transfer-based estimation of glymphatic transport in living brain. In Medical Imaging 2015: Image Processing, SPIE (2015) ] early on to study mouse brain transport via in-vivo imaging.

However, when we applied this method to our human image data set, we quickly faced several challenges. First, by considering transport by convection alone (3), we ignore the substantial contribution from diffusion. Second, we found the method sensitive to the choice of the regularization parameter α>0\alpha>0. And third, c0c\approx 0 locally, for which (2) is not obviously well posed. We therefore needed to develop a more physiologically accurate and more numerically robust approach.

Bilinear flow control

Still targeting the task of “given c1,c2 ⁣:ΩRc_{1},c_{2}\colon\Omega\to\nobreak\mathbb{R} at t1,t2t_{1},t_{2}, identify a convective flow field φ\varphi defined over the time interval [t1,t2][t_{1},t_{2}]”, we instead considered the bilinear optimal control problem of finding minimizers

minc,φ(c(t2)c22+α2R2),\min_{c,\varphi}\bigl(\lVert c(t_{2})-c_{2}\rVert^{2}+\alpha^{2}R^{2}\bigr),

subject to c(t1)=c1c(t_{1})=c_{1}, and such that (1) holds (with r=0r=0) over (t1,t2](t_{1},t_{2}]. We let R=φL2+φH1R=\lVert\varphi\rVert_{L^{2}}+\lvert\varphi\rvert_{H^{1}}, though that choice can easily be reevaluated. Note that we deliberately allow for non-divergence-free velocity fields in light of our previous results (Section 2.1) indicating that local fluid influx and efflux may be key for effective convective transport. This problem setting is very similar to that studied in terms of well-posedness and existence by Glowinski et al. [20 R. Glowinski, Y. Song, X. Yuan and H. Yue, Bilinear optimal control of an advection-reaction-diffusion system. SIAM Rev. 64, 392–421 (2022) ] modulo choice of regularization and incompressibility (divergence-free) constraint (or lack thereof).

There are two main approaches to solving the constrained minimization problem (4): either (i) introducing a Lagrange multiplier for the PDE-constraint and solving the resulting nonlinear Euler–Lagrange problems directly (the all-at-once approach), or (ii) iteratively solving the optimization problem by solving the PDE-constraint (1) with c1c_{1} as initial condition for a series of φ\varphi by means, e.g., a quasi-Newton method (the reduced approach). After numerical experiments and comparisons, we found that the latter approach (using an L-BFGS algorithm) was more robust for our uses in the sense that this method successfully converged for most patients and time intervals between brain scans and yielded consistent results for different mesh resolutions and values of the regularization parameters.

Remark 5. Our previous work on automated solution of high-dimensional inverse problems and accompanying software dolfin-adjoint [18 P. E. Farrell, D. A. Ham, S. W. Funke and M. E. Rognes, Automated derivation of the adjoint of high-level transient finite element programs. SIAM J. Sci. Comput. 35, C369–C393 (2013) , 2 M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes and G. N. Wells, The FEniCS project version 1.5. Archive of Numerical Software 3, 9–23 (2015) ] proved to be very useful here.

Ultimately, we thus decided to estimate net (time-averaged) velocity fields φ ⁣:ΩRd\varphi\colon\Omega\to\mathbb{R}^{d} in the contrast agent influx phase (t1,t2)=(6,24)(t_{1},t_{2})=(6,24) h and clearance phase (t1,t2)=(24,48)(t_{1},t_{2})=(24,48) h for each patient by solving (4) (using subject-specific finite element meshes, diffusion tensors and concentration observations) via a reduced approach using a fixed maximal number of iterations (Figure 2 (b)). Interestingly, for all subjects and time intervals, we identified persistent velocity fields with brain-averaged velocity magnitudes (i.e., flow speeds) of \sim1–8 µm/min. These flow speeds correspond to bulk flow rates of 0.1–0.8 µL/(g min) – a range which compares fascinatingly well with estimated bulk flow rates of 0.1–0.3 µL/(g min) from early (1970s–80s) tracer experiments in cats [1 N. J. Abbott, Evidence for bulk flow of brain interstitial fluid: Significance for physiology and pathology. Neurochem. Int. 45, 545–552 (2004) ]. Interestingly, we can also estimate the net fluid influx/efflux by computing, e.g., the brain-wide average of divφ\operatorname{div}\varphi, and find that water dominantly flows into the brain (e.g., by filtration from the blood stream) at a rate of \sim1–4 min.

Open problem 6. Our findings suggest that high-dimensional inverse modelling offers a powerful avenue of investigation for the brain’s waterscape. However, a question that remains is how well-posed the optimization problem is with respect to the uncertainty in the input data and choice of regularization functional. Moreover, clearly there are “terra incognita” (or perhaps rather “mare incognita”) for which there is little information in the medical images, i.e., little or no contrast agent present, and where velocity field estimates intuitively would be associated with substantial uncertainty. Quantifying the uncertainty and information level required for reliable flow estimates in this setting is a open problem.

So, using stochastic and inverse modelling as described, we have quantified the relative contributions of diffusion and convection in brain solute transport, suggesting that these modes coexist and co-contribute. The next question is: given that there is a persistent convective fluid flow within the brain parenchyma, what could be the mechanisms and the brain mechanics allowing for such a flow?

3 The brain as a pulsating poroelastic medium

From the viewpoint of mechanics, the brain is an almost surprisingly soft elastic medium, enclosed and protected by the stiffer meninges and much stiffer skull, and permeated by a number of fluid networks including the blood-filled vasculature (arteries, capillaries and veins), the ISF-filled extracellular space between brain cells, and potential CSF- or ISF-filled perivascular spaces surrounding the vasculature. With every heart beat (the cardiac rhythm) and breath (the respiratory cycle), your brain expands and contracts with displacements of the order of 1 mm and volume changes of the order of 1 cm. The precise interplay and transfer of forces and sources between these compartments, and how these change with age and disorders, still lacks understanding and quantification however.

Figure 3. Top: Brain displacements and arterial blood network pressure predicted by a pulsatile fluid influx and the MPET theory (5) (snapshot in time). Below: A posteriori error indicators (sagittal view) on initial (left) and adaptively refined final (right) mesh [33].

In pioneering work, Tully and Ventikos [34 B. Tully and Y. Ventikos, Cerebral water transport using multiple-network poroelastic theory: Application to normal pressure hydrocephalus. J. Fluid Mech. 667, 188–215 (2011) ] introduced the multiple-network poroelastic theory (MPET), a theory originating in fractured geological reservoirs, in the context of modelling brain solid and fluid mechanics. Specifically, the quasi-static MPET equations enforce balance of momentum and mass and read as: for a given number of networks JNJ\in\mathbb{N}, ΩRd\Omega\subset\mathbb{R}^{d} (d=1,2,3d=1,2,3), and a time interval II find the displacement field u ⁣:Ω×IRdu\colon\Omega\times I\to\mathbb{R}^{d} and the jj-th network pressure field pj:Ω×Ip_{j}:\Omega\times I for j=1,,Jj=1,\ldots,J such that

div(σ(u))+1emjαjgradpj=f,\displaystyle-\operatorname{div}(\sigma(u))+\sum\nolimits_{1emj}\alpha_{j}\operatorname{grad}p_{j}=f,
t(sjpj+αjdivu)divκjgradpj+Tj(p)=gj.\displaystyle\partial_{t}(s_{j}p_{j}+\alpha_{j}\operatorname{div}u)-\operatorname{div}\kappa_{j}\operatorname{grad}p_{j}+T_{j}(p)=g_{j}.

Here σ(u)\sigma(u) is the elastic stress-strain relationship and TT describes the transfer of mass between fluid networks. In the linear and isotropic case,

σ(u)=2με(u)+λdiv(u)I,ε(u)12(gradu+graduT),\displaystyle\sigma(u)=2\mu\varepsilon(u)+\lambda\operatorname{div}(u)I,\quad\varepsilon(u)\equiv\smash[b]{\frac{1}{2}}(\operatorname{grad}u+\operatorname{grad}u^{\mathsf{T}}),
Tj(p)=1emiξji(pjpi).\displaystyle T_{j}(p)=\sum\nolimits_{1emi}\xi_{ji}(p_{j}-p_{i}).

with Lamé parameters μ>0\mu>0 and λ>0\lambda>0. Each network jj is equipped with its Biot–Willis coefficient αj(0,1]\alpha_{j}\in(0,1] with jαj1\sum_{j}\alpha_{j}\leqslant 1, its storage coefficient sj>0s_{j}>0, and its hydraulic conductivity tensor κj\kappa_{j}, while the network interactions are described via the exchange coefficients ξji\xi_{ji}. We also note that the fluid (Darcy) velocity vjv_{j} in network jj is defined by

vj=κjgradpj.v_{j}=-\kappa_{j}\operatorname{grad}p_{j}.

For simplicity, consider the basic boundary conditions u=0u=0, pj=0p_{j}=0 on Ω\partial\Omega and initial conditions pj(0)=pj,0p_{j}(0)=p_{j,0} here.

In the case J=1J=1, the MPET equations reduce to Biot’s equations, whose properties have been studied for decades (and still actively are). The general MPET equations however seemed to have received much less attention from the mathematical community.

3.1 MPET as an elliptic-parabolic problem

To analyse approximations of the MPET equations (Figure 3), we realized that the framework of coupled elliptic-parabolic problems as introduced by Ern and Meunier [17 A. Ern and S. Meunier, A posteriori error analysis of Euler–Galerkin approximations to coupled elliptic-parabolic problems. M2AN Math. Model. Numer. Anal. 43, 353–375 (2009) ] provided an ideal starting point. Specifically, this framework considers variational problems in space and time of the form “given bilinear forms aa, bb, cc, and dd and input data ff and gg, find uH1(I;Va)u\in H^{1}(I;V_{a}) and pH1(I;Vd)p\in H^{1}(I;V_{d}) such that, for almost every tIt\in I,

a(u,v)b(v,p)=f,vvVa,\displaystyle a(u,v)-b(v,p)=\langle f,v\rangle_{*}\quad\forall v\in V_{a},
c(pt,q)+b(ut,q)+d(p,q)=g,qqVd,\displaystyle c(p_{t},q)+b(u_{t},q)+d(p,q)=\langle g,q\rangle_{*}\quad\forall q\in V_{d},

with the initial condition p(0)=p0p(0)=p_{0}”. Crucially, under natural assumptions on the bilinear forms and the underlying (Hilbert) spaces, and subsequently spatial and temporal discretizations, Ern and Meunier [17 A. Ern and S. Meunier, A posteriori error analysis of Euler–Galerkin approximations to coupled elliptic-parabolic problems. M2AN Math. Model. Numer. Anal. 43, 353–375 (2009) ] derived optimal stability estimates, a priori error estimates and a posteriori error estimates for such problems.

Introducing Va=H01(Ω;Rd)V_{a}=H^{1}_{0}(\Omega;\mathbb{R}^{d}), Vd=H01(Ω;RJ)V_{d}=H^{1}_{0}(\Omega;\mathbb{R}^{J}) and denoting p=(p1,,pJ)p=(p_{1},\dots,p_{J}) and analogously for α\alpha, the MPET equations (5) take the form of a coupled elliptic-parabolic problem (7) via the identifications

a(u,v)=σ(u),ε(v),\displaystyle a(u,v)=\langle\sigma(u),\varepsilon(v)\rangle,
b(u,p)=αp,divu,\displaystyle b(u,p)=\langle\alpha\cdot p,\operatorname{div}u\rangle,
c(p,q)=j=1Jsj1empj,qj,\displaystyle c(p,q)=\sum_{j=1}^{J}\langle s_{j}1emp_{j},q_{j}\rangle,
d(p,q)=j=1Jκjgradpj,gradqj+Tj,qj,\displaystyle d(p,q)=\sum_{j=1}^{J}\langle\kappa_{j}\operatorname{grad}p_{j},\operatorname{grad}q_{j}\rangle+\langle T_{j},q_{j}\rangle,

Clearly, VaV_{a} and VdV_{d} as defined are Hilbert spaces, dense in La=L2(Ω;Rd)L_{a}=L^{2}(\Omega;\mathbb{R}^{d}) and Ld=L2(Ω;RJ)L_{d}=L^{2}(\Omega;\mathbb{R}^{J}), respectively, aa is symmetric, continuous and coercive if μ>0\mu>0, 2λ+dμ>02\lambda+d\mu>0), bb is continuous for αj(0,1]\alpha_{j}\in(0,1] and cc is symmetric, continuous and coercive for sj>0s_{j}>0. By assuming that the exchange coefficients are positive, symmetric ξji=ξij>0\xi_{ji}=\xi_{ij}>0, and bounded, and that the conductivities are bounded (including from below: κj>κmin>0\kappa_{j}>\kappa_{\mathrm{min}}>0), shuffling terms and the Poincaré inequality gives that dd is symmetric, coercive and continuous, then stability and uniqueness follows for (weak) solutions of the MPET equations (5) [17 A. Ern and S. Meunier, A posteriori error analysis of Euler–Galerkin approximations to coupled elliptic-parabolic problems. M2AN Math. Model. Numer. Anal. 43, 353–375 (2009) ].

Introducing a conforming finite element discretization of VaV_{a} and VdV_{d} relative to a mesh Th\mathcal{T}_{h} with mesh size hh and a first-order implicit time discretization with time steps τn\tau_{n} to define discrete approximations uhτu_{h\tau} and phτp_{h\tau}, a priori and a posteriori error estimates (Figure 3) also follow.

Theorem 7 (Eliseussen, Rognes and Thompson [12 E. Eliseussen, M. E. Rognes and T. B. Thompson, A-posteriori error estimation and adaptivity for multiple-network poroelasticity, preprint, arXiv:2111.13456 (2021) ]).

For solutions (u,p)(u,p) and approximations (uhτ,phτ)(u_{h\tau},p_{h\tau}) of the MPET equations, and for each time step tnt_{n}, n{0,1,,N}n\in\{0,1,\dots,N\},

uuhτL(0,tn;H1)+pphτL(0,tn;L2)+pphτL2(0,tn;H1)η1+η2+η3+η4+Eh0(u0,p0)+Eh(f,g),\begin{split}&\lVert u-u_{h\tau}\rVert_{L^{\infty}(0,t_{n};H^{1})}+\lVert p-p_{h\tau}\rVert_{L^{\infty}(0,t_{n};L^{2})}+\lVert p-p_{h\tau}\rVert_{L^{2}(0,t_{n};H^{1})}\\ &\qquad\lesssim\eta_{1}+\eta_{2}+\eta_{3}+\eta_{4}+\mathcal{E}_{h0}(u_{0},p_{0})+\mathcal{E}_{h}(f,g),\end{split}

where η1,η2,η3,η4\eta_{1},\eta_{2},\eta_{3},\eta_{4} are a posteriori computable quantities involving appropriately weighted cell and facet residual contributions, Eh0\mathcal{E}_{h0} and Eh\mathcal{E}_{h} are determined by the approximation of the initial data and source terms, respectively.

3.2 The incompressible MPET equations

As the human brain is composed of \approx 80 % water, it is widely considered (nearly) incompressible (λ\lambda\to\infty in (5)). It therefore seems highly relevant to ask if the variational formulation (7) with forms given by (8) and the estimates of Theorem 7 are robust with respect to λ\lambda? The answer is no: the continuity of the linear elasticity form aa (and thus the error estimates) degenerate as λ\lambda\to\infty. Indeed, it is illuminating to take a look at the structure of the MPET equations in this scenario. Set sj=0s_{j}=0 to reveal the extreme cases. As λ\lambda\to\infty, divu0\operatorname{div}u\to 0, the MPET equations decouple into a coupled Darcy flow system for p=(p1,,pJ)p=(p_{1},\dots,p_{J}) and an elliptic equation for uu,

divκjgradpj+Tj(p)=gj,\displaystyle-\operatorname{div}\kappa_{j}\operatorname{grad}p_{j}+T_{j}(p)=g_{j},
div2με(u)=f1emjαjgradpj.\displaystyle\notag-\operatorname{div}2\mu\varepsilon(u)=f-\sum\nolimits_{1emj}\alpha_{j}\operatorname{grad}p_{j}.

We can thus expect finite element discretization of the MPET equations in the incompressible limit to be wrought with analogous challenges as for the linear elasticity equations.

How should we remedy this situation? An appealing solution, first introduced in the context of Biot’s equations, is to introduce the total pressurep0p_{0} as a new variable, p0=λdivuαpp_{0}=\lambda\operatorname{div}u-\alpha\cdot p. Then the action of the MPET operator transforms into

( ⁣div2με ⁣grad0divλ1λ1α0λ1αTtC~+λ1ααTt)(up0p),\begin{pmatrix}-\!\operatorname{div}2\mu\varepsilon&-\!\operatorname{grad}&0\\ \operatorname{div}&\lambda^{-1}&\lambda^{-1}\alpha\\ 0&\lambda^{-1}\alpha^{\mathsf{T}}\partial_{t}&\tilde{C}+\lambda^{-1}\alpha\alpha^{\mathsf{T}}\partial_{t}\end{pmatrix}\begin{pmatrix}u\\ p_{0}\\ p\end{pmatrix},

This formulation is indeed robust in the incompressible limit with energy estimates uniform in λ\lambda and (subsequently a priori error estimates for stable discretizations):

Theorem 8 (Lee et al. [23 J. J. Lee, E. Piersanti, K.-A. Mardal and M. E. Rognes, A mixed finite element method for nearly incompressible multiple-network poroelasticity. SIAM J. Sci. Comput. 41, A722–A747 (2019) ]).

Given sufficiently regular ff, gjg_{j} and initial conditions I0I_{0}, solutions uu, pjp_{j} to system (5) satisfy an energy estimate of the form

ε(u)L(I,L2(Ω))+1emj(pjL(I,sjL2(Ω))+pjL2(I,κjH01(Ω)))I0+f,f˙,gj,g˙j,\begin{split}&\lVert\varepsilon(u)\rVert_{L^{\infty}(I,L^{2}(\Omega))}+\sum\nolimits_{1emj}(\lVert p_{j}\rVert_{L^{\infty}(I,s_{j}L^{2}(\Omega))}+\lVert p_{j}\rVert_{L^{2}(I,\kappa_{j}H^{1}_{0}(\Omega))})\\ &\qquad\lesssim I_{0}+\lVert f,\dot{f},g_{j},\dot{g}_{j}\rVert,\end{split}

with inequality constant independent of the parameter λ\lambda.

Preconditioning by congruence

When next turning to the efficient solution of the MPET equations via preconditioned iterative methods, an interesting puzzle appeared with a charming solution motif [27 E. Piersanti, J. J. Lee, T. Thompson, K.-A. Mardal and M. E. Rognes, Parameter robust preconditioning by congruence for multiple-network poroelasticity. SIAM J. Sci. Comput. 43, B984–B1007 (2021) , 28 E. Piersanti, M. E. Rognes and K.-A. Mardal, Parameter robust preconditioning for multi-compartmental Darcy equations. In Numerical mathematics and advanced applications—ENUMATH 2019, Lect. Notes Comput. Sci. Eng. 139, Springer, Cham, 703–711 (2021) ]. To illustrate, consider a coupled Darcy flow with exchange system like (9) written as

Ap=(KΔ+E)p=b,Ap=(-K\Delta+E)p=b,

where K=diag(κ1,,κJ)K=\mathrm{diag}(\kappa_{1},\dots,\kappa_{J}) and

Eij=ξijfor ij,Eii=ξj,ξj=1emiξji.E_{ij}=-\xi_{ij}\quad\textrm{for}\ i\neq j,\qquad E_{ii}=\xi_{j},\quad\xi_{j}=\sum\nolimits_{1emi}\xi_{ji}.

Block-diagonal preconditioners BB of this system easily yield condition numbers that grow with the ratio between the exchange ξji\xi_{ji} and conductance κj\kappa_{j} coefficients [28 E. Piersanti, M. E. Rognes and K.-A. Mardal, Parameter robust preconditioning for multi-compartmental Darcy equations. In Numerical mathematics and advanced applications—ENUMATH 2019, Lect. Notes Comput. Sci. Eng. 139, Springer, Cham, 703–711 (2021) ]. But, can the task of constructing block-diagonal preconditioners be simplified by a linear change of variables pp~p\mapsto\tilde{p} with a matrix PP (1emp=Pp~p=P\tilde{p})?

By definition, a matrix AA is diagonalizable by congruence if and only if there exists a matrix PP such that PTAPP^{\mathsf{T}}AP is diagonal. Further, two matrices AA and BB are simultaneously diagonalizable by congruence if there exists a matrix PP such that both PTAPP^{\mathsf{T}}AP and PTBPP^{T}BP are diagonal. Now, if AA is a real, symmetric and positive definite and BB is a real symmetric matrix, then AA and BB are simultaneously diagonalizable by congruence. In our case, KK is symmetric and positive definite as long as κj>0\kappa_{j}>0 and EE is symmetric; therefore, KK and EE are simultaneously diagonalizable by congruence.

Next, in general, if C=A1BC=A^{-1}B has JJ distinct eigenvalues with eigenvectors vjv_{j}, then the matrix P={vj}j=1JP=\{v_{j}\}_{j=1}^{J} constructed from these eigenvectors diagonalizes AA and BB by congruence. This observation gives a direct construction for PP in terms of the eigenvectors of K1EK^{-1}E. The transformed system is then block-diagonal and easily preconditioned with condition numbers that are uniform for a range of exchange and conductance parameters [28 E. Piersanti, M. E. Rognes and K.-A. Mardal, Parameter robust preconditioning for multi-compartmental Darcy equations. In Numerical mathematics and advanced applications—ENUMATH 2019, Lect. Notes Comput. Sci. Eng. 139, Springer, Cham, 703–711 (2021) ]. Interestingly, the construction extends to the MPET system [27 E. Piersanti, J. J. Lee, T. Thompson, K.-A. Mardal and M. E. Rognes, Parameter robust preconditioning by congruence for multiple-network poroelasticity. SIAM J. Sci. Comput. 43, B984–B1007 (2021) ].

Figure 4. Intracranial pressure and fluid flow in the brain and surrounding CSF subject to a pulsatile fluid (blood) influx in the brain – a fluid-structure interaction problem [6, Figure 4].

3.3 Brain-CSF interactions

Remark 9. The total pressure formulation of the MPET equations is also advantageous when considering the fluid-structure interaction between the poroelastic brain and the CSF that surrounds it and clinical applications (Figure 4) [6 M. Causemann, V. Vinje and M. E. Rognes, Human intracranial pulsatility during the cardiac cycle: A computational modelling framework. Fluids Barriers CNS 19, 863–874 (2022) ]. The intracranial pressure, a key clinical measure, can be interpreted as the total pressure in the parenchyma and the hydrostatic pressure in the CSF.

Open problem 10. In-vivo imaging options for brain pulsatility (deformations, stresses) are currently limited, but new techniques such as, e.g., magnetic resonance encephalography (MREG) are coming into play. The quantification of MREG images remains somewhat enigmatic however, and we envision that computational fluid-structure poroelastic brain simulations could provide new interpretations.

3.4 Impermeability and minimal Biot–Stokes stability

The permeability of the brain’s extracellular space is estimated to be of the order 10 nm, which is many orders of magnitude lower than, e.g., the permeability of the blood circulation. Thus, another limit of interest for studying brain fluid flow is the impermeable case κ0\kappa\to 0. This is a particularly interesting case, often attempted addressed by introducing the Darcy velocity (as defined by (6)) as a separate field.

Taking the (time-discrete) Biot equations (i.e., (5) with J=1J=1) and s1=0s_{1}=0 as an example; these form a generalized saddle-point system that in the limit as κ=0\kappa=0 reduces to the incompressible Stokes equations in terms of uu and pp. On the other hand as noted previously (cf. (9)), in the limit as λ\lambda\to\infty, the Biot equations decouple into an elliptic equation for uu and a (mixed) Darcy equation for (zz and) pp. These observations hint at close relations between Stokes, Darcy, Biot and MPET equations, and suggest that, for combinations of finite element spaces U×W×PU\times W\times P, that the pairing U×PU\times P is Stokes stable and W×PW\times P Darcy stable in the discrete Babuška–Brezzi sense. However, in the limit at κ0\kappa\to 0, it is highly non-trivial to ensure uniform Darcy stability [24 K.-A. Mardal, M. E. Rognes and T. B. Thompson, Accurate discretization of poroelasticity without Darcy stability: Stokes–Biot stability revisited. BIT 61, 941–976 (2021) ].

Figure 5. Brain blood vessels are surrounded by fluid-filled perivascular spaces. Are these dual networks amenable for geometrical model reduction? Pulsatile Stokes flow in perivascular spaces can accurately and efficiently be represented by geometrically reduced models over the centreline geometry [9, Figure 6].

Fortunately, we could show that a revised notion, coined minimal Stokes–Biot stability is sufficient.

Definition 11 (Minimal Stokes–Biot stability [24 K.-A. Mardal, M. E. Rognes and T. B. Thompson, Accurate discretization of poroelasticity without Darcy stability: Stokes–Biot stability revisited. BIT 61, 941–976 (2021) ]). A family of discrete spaces {Uh×Wh×Qh}h\{U_{h}\times W_{h}\times Q_{h}\}_{h} with UhUU_{h}\subset U, WhWW_{h}\subset W and QhQQ_{h}\subset Q is called minimally Stokes–Biot stable if and only if

  1. the bilinear form aa is continuous and coercive on UhU_{h};

  2. {Uh×Qh}h\{U_{h}\times Q_{h}\}_{h} are Stokes stable in the discrete Babuška–Brezzi sense;

  3. divWhQh\operatorname{div}W_{h}\subseteq Q_{h} for each hh.

Theorem 12 (Mardal, Rognes and Thompson [24 K.-A. Mardal, M. E. Rognes and T. B. Thompson, Accurate discretization of poroelasticity without Darcy stability: Stokes–Biot stability revisited. BIT 61, 941–976 (2021) ]).

Minimally stable Stokes–Biot discrete solutions at each time step mm converge,

umuhm1+zmzhmK1,τ;div+pmphmhcM1+τM2\lVert u^{m}-u_{h}^{m}\rVert_{1}+\lVert z^{m}-z_{h}^{m}\rVert_{K^{-1},\tau;\operatorname{div}}+\lVert p^{m}-p_{h}^{m}\rVert\lesssim h^{c}M_{1}+\tau M_{2}

for cNc\in\mathbb{N} (depending on regularity, order, etc.), where

M1,M2(u,z,p,tu,ttu,tp,tz).M_{1},M_{2}\lesssim(\lVert u\rVert,\lVert z\rVert,\lVert p\rVert,\lVert\partial_{t}u\rVert,\lVert\partial_{tt}u\rVert,\lVert\partial_{t}p\rVert,\lVert\partial_{t}z\rVert).

Now, returning to the physiology at hand, the moderate pressure differences observed clinically and the nearly vanishing permeability of the extracellular space suggest that extracellular fluid flow is negligible. Could there be alternative, more permeable pathways within the brain?

4 Perivascular fluid flow and model reduction

As the brain lacks lymphatic vessels, fluid-filled spaces surrounding blood vessels – the perivascular spaces (PVSs) – are conjectured to act as higher-permeability proxy pathways; “highways” for fluid efflux and metabolic solute clearance [31 M. L. Rennels, T. F. Gregory, O. R. Blaumanis, K. Fujimoto and P. A. Grady, Evidence for a ‘paravascular’ fluid circulation in the mammalian central nervous system, provided by the rapid distribution of tracer protein throughout the brain from the subarachnoid space. Brain Res. 326, 47–63 (1985) , 22 J. J. Iliff, M. Wang, Y. Liao, B. A. Plogg, W. Peng et al., A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid β. Sci. Transl. Med. 4, DOI 10.1126/scitranslmed.3003748 (2012) ]. At this mesoscale (\approx mm, s), the vasculature and perivasculature form slender dual networks, running along the brain’s surface and diving into the brain. While experimental studies indicate that solutes move rapidly in PVSs, a lingering question is “what are the mechanisms, characteristics and forces underlying PVS fluid flow”?

Flow in perivascular spaces appears synchronized with the cardiac cycle. This observation has led to the widespread notion that perivascular flow is driven by arterial wall motion. However, can pulsatile wall motions be a driver for net flow at this scale? To investigate, we created computational models of an annular perivascular segment surrounding an image-based bifurcating blood vessel (Figure 5) and induced laminar (Stokes) flow in the moving PVS by (i) travelling pulse waves on the (inner) blood vessel wall, (ii) pulsatile pressure differences at inlet and outlets, and (iii) static pressure differences. Wall pulsations induce substantial pulsatile flow, but only negligible net flow (due to their long wavelength of \sim100 mm) [11 C. Daversin-Catty, V. Vinje, K.-A. Mardal and M. E. Rognes, The mechanisms behind perivascular fluid flow. PLOS ONE 15, Paper No. e0244442 (2020) ]. On the other hand, even a small static hydrostatic pressure difference can induce net flow of relevant magnitudes. Such hydrostatic pressure differences could be induced by experimental procedures [37 V. Vinje, A. Eklund, K.-A. Mardal, M. E. Rognes and K.-H. Støverud, Intracranial pressure elevation alters CSF clearance pathways. Fluids Barriers CNS 17, 1–19 (2020) ], but are also of comparable magnitude as transient pressure differences measured in the human brain [38 V. Vinje, G. Ringstad, E. K. Lindstrøm, L. M. Valnes, M. E. Rognes, P. K. Eide and K.-A. Mardal, Respiratory influence on cerebrospinal fluid flow: A computational study based on long-term intracranial pressure measurements. Sci. Rep. 9, 1–13 (2019) ].

Figure 6. Fluid flow velocity profiles in idealized (A) and image-based (B) periarterial and perivenous geometries [36, Figure 4]. The slow perivenous flow results in delayed tracer transport around veins.

Another disputed point is whether there is a net efflux of ISF along perivenous spaces. The original glymphatic theory emphasizes this concept, as injected tracers were observed along arteries, but not around veins at early time points, and along veins at later time points [22 J. J. Iliff, M. Wang, Y. Liao, B. A. Plogg, W. Peng et al., A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid β. Sci. Transl. Med. 4, DOI 10.1126/scitranslmed.3003748 (2012) ]. However, could an alternative explanation for this be connected to the geometry of brain surface arteries and veins? We created pressure-driven Stokes flow models as well as convection-diffusion models of tracer transport from optical coherence tomography images of perivascular spaces surrounding arteries and veins to investigate [36 V. Vinje, E. N. Bakker and M. E. Rognes, Brain solute transport is more rapid in periarterial than perivenous spaces. Sci. Rep. 11, 1–11 (2021) ]. PVS flow speeds were 2–6 times higher in the periarterial geometries than in perivenous geometries (Figure 6). Interestingly, these differences in flow speeds due to geometrical differences (area, shape) lead to delayed tracer transport by about 25–30 min, in agreement with the originally observed delayed tracer distribution surrounding veins.

4.1 Perivascular spaces as topologically 1d-networks

Modelling the interplay between blood vessels and tissue via geometrically reduced models has been an active and important research topic for decades, cf., e.g., [5 L. Cattaneo and P. Zunino, Computational models for fluid exchange between microcirculation and tissue interstitium. Netw. Heterog. Media 9, 135–159 (2014) ] and references therein. In the brain, the characteristic dual networks formed by the cerebral vasculature and perivasculature define a new setting for geometrically reduced fluid flow modelling [9 C. Daversin-Catty, I. G. Gjerde and M. E. Rognes, Geometrically reduced modelling of pulsatile flow in perivascular networks. Front. Phys. (Lausanne), DOI 10.3389/fphy.2022.882260 (2022) ].

Consider a perivascular domain Ω\Omega defined as the union of generalized annular cylinder branches Ωi\Omega^{i}. Assume that each such cylinder has a well-defined, oriented, and topologically one-dimensional centreline Λi\Lambda^{i} with coordinate ss, and let Λ=iΛi\Lambda=\bigcup_{i}\Lambda^{i} (Figure 5). The set of bifurcation points i.e., the points at which the centrelines of branches meet is denoted B\mathcal{B}. Finally, assume that the boundaries pulsate with a wall speed ww in the normal direction.

Under a set of assumptions on the moving geometries (axial symmetry, radial boundary movement and fixed centreline) and flow and pressure fields (axial symmetry, constant cross-section pressure, axial velocity profile), we can derive a reduced system of equations for the cross-section flux qq and average cross-section pressure pp over time t>0t>0 on each centreline Λi\Lambda^{i} with inner and outer radius R1R_{1}, R2R_{2} (denoting qΛiq|_{\Lambda^{i}} by qiq^{i} and pΛip|_{\Lambda^{i}} by pip^{i}),

ρAiqtiμAiqssi+μαiAiqi+psi=0\displaystyle\smash[b]{\frac{\rho}{A^{i}}q_{t}^{i}-\frac{\mu}{A^{i}}q_{ss}^{i}+\mu\frac{\alpha^{i}}{A^{i}}q^{i}+p_{s}^{i}}=0
on Λi,\displaystyle\textrm{on}\ \Lambda^{i},
qsi=fi\displaystyle q_{s}^{i}=f^{i}
on Λi,\displaystyle\textrm{on}\ \Lambda^{i},

where

fi(s)2π(R1i(s,t)w(R1,s,t)+R2i(s,t)w(R2,s,t)).f^{i}(s)\equiv 2\pi\bigl(R_{1}^{i}(s,t)w(R_{1},s,t)+R_{2}^{i}(s,t)w(R_{2},s,t)\bigr).

In (11), ρ\rho is the density of the fluid and μ\mu its viscosity, Ai=Ai(s,t)A^{i}=A^{i}(s,t) denotes the cross-section area, while α^i=α^i(s,t)\hat{\alpha}^{i}=\hat{\alpha}^{i}(s,t) is a lumped flow parameter that depends on the domain geometry and the choice of axial velocity profile. We also define the (one-dimensional) normal stress induced by q^\hat{q} and p^\hat{p} as

σμAqsp,\sigma\equiv\smash[t]{\frac{\mu}{A}}q_{s}-p,

which corresponds to an average of the axial (ss-)component of the normal stress in each cross-section.

At the bifurcation points bBΩb\in\mathcal{B}\subset\Omega, we impose conservation of flux and continuity of normal stress,

qp(sp)=qd1(sd1)+qd2(sd2),\displaystyle q^{p}(s^{p})=q^{d_{1}}(s^{d_{1}})+q^{d_{2}}(s^{d_{2}}),
σp(sp)=σd1(sd1)=σd2(sd2),\displaystyle\sigma^{p}(s^{p})=\sigma^{d_{1}}(s^{d_{1}})=\sigma^{d_{2}}(s^{d_{2}}),

where Λp\Lambda^{p} and Λd1\Lambda^{d_{1}}, Λd2\Lambda^{d_{2}} represent centrelines of the parent and daughter branches meeting at the bifurcation point bb, and ss^{\cdot} the respective coordinates of bb.

This system is particularly amenable for a variational finite element formulation imposing the flux conservation condition weakly using a Lagrange multiplier. Relative to a finite element mesh T\mathcal{T} of the centreline Λ\Lambda composed of mesh segments Ti\mathcal{T}^{i} (one for each centreline branch Λi\Lambda^{i}), we define the

  • flux space VhV_{h} as the space of continuous piecewise quadratics over Ti\mathcal{T}^{i} for each ii,

  • pressure space QhQ_{h} as the space of continuous piecewise linears on T\mathcal{T}, and

  • Lagrange multiplier space Rh=RBR_{h}=\mathbb{R}^{\lvert\mathcal{B}\rvert}.

The discrete problem then reads as follows: for each discrete time tkt^{k} and time step τ\tau, find qhVhq_{h}\in V_{h}, phQhp_{h}\in Q_{h} and λhRh\lambda_{h}\in R_{h} solving

a((qh,ph,λh),(ψ,φ,ξ05em))=Lk((ψ,φ,ξ05em))a\bigl((q_{h},p_{h},\lambda_{h}),(\psi,\varphi,\xi 05em)\bigr)=L^{k}\bigl((\psi,\varphi,\xi 05em)\bigr)

for all ψVh\psi\in V_{h}, φQh\varphi\in Q_{h}, and ξRh\xi\in R_{h}. Here, aa is defined by

a((q,p,λ),(ψ,φ,ξ05em))=iIΛi(Ciqiψi+τμAiqsiψsi+qsiφiτψsipi)ds+bBλb[ψ]b+ξb[q]b,\begin{split}&a\bigl((q,p,\lambda),(\psi,\varphi,\xi 05em)\bigr)\\ &\qquad=\mathop{\smash[b]{\sum_{i\in I}}}\int_{\Lambda^{i}}\biggl(C^{i}q^{i}\psi^{i}+\frac{\tau\mu}{A^{i}}q_{s}^{i}\psi_{s}^{i}+q_{s}^{i}\varphi^{i}-\tau\psi_{s}^{i}p^{i}\biggr)\,\mathrm{d}s\\ &\qquad\qquad+\sum_{b\in\mathcal{B}}\lambda^{b}[\psi]^{b}+\xi^{b}[q]^{b},\end{split}

where Ci=Ai1(ρ+τμαi)C^{i}=A_{i}^{-1}(\rho+\tau\mu\alpha^{i}), and λb\lambda^{b} (and ξb\xi^{b}) is simply the λ\lambda (or ξ\xi1em) corresponding to the point bb, and the natural jump is

[ψ]b=ψp(b)ψd1(b)ψd2(b).[\psi]^{b}=\psi^{p}(b)-\psi^{d_{1}}(b)-\psi^{d_{2}}(b).

This formulation provides an inexpensive and reasonably accurate framework for estimating pulsatile perivascular fluid flow in large networks, and thus establishes a mathematical foundation for future computational studies of perivascular flow and transport (Figure 5) [9 C. Daversin-Catty, I. G. Gjerde and M. E. Rognes, Geometrically reduced modelling of pulsatile flow in perivascular networks. Front. Phys. (Lausanne), DOI 10.3389/fphy.2022.882260 (2022) ].

Open problem 13. Originally, I was planning on coupling the vasculature and perivascular spaces with the tissue via topologically one-dimensional dual networks embedded in the three-dimensional volume, see e.g., [5 L. Cattaneo and P. Zunino, Computational models for fluid exchange between microcirculation and tissue interstitium. Netw. Heterog. Media 9, 135–159 (2014) ]. Due to the controversy regarding the existence of flow in brain PVSs and extracellular space, we did not go further with this approach, and as such it remain as (at least) an interesting mathematical problem.

5 Bridging electrochemistry and fluid mechanics

When you are thinking about the brain, the first thing that comes to mind is probably not brain fluid flow. Indeed, neuroscience is possibly foremost occupied with the electrical signalling of the brain. These signals are however generated by electric potential differences, which in turn are induced by ion concentration differences [29 R. Rasmussen, J. O’Donnell, F. Ding and M. Nedergaard, Interstitial ions: A key regulator of state-dependent neural activity? Prog. Neurobiol. 193, DOI 10.1016/j.pneurobio.2020.101802 (2020) ]. On the other hand, differences in ion concentrations also induce water movement by osmosis. Thus, up till now, we have only considered one (mechanical) piece of the puzzle, while ignoring electrical and chemical effects. The challenge is that

the coupling between electrochemical effects, fluid transport and elastic deformation is particularly difficult and is only little understood, specifically in the brain [19 A. Goriely, M. G. Geers, G. A. Holzapfel, J. Jayamohan, A. Jérusalem, S. Sivaloganathan, W. Squier et al., Mechanics of the brain: Perspectives, challenges, and opportunities. Biomech. Model. Mechanobiol. 14, 931–965 (2015) ].

By a pure serendipity, just as we were getting ready to address this challenge, Mori [26 Y. Mori, A multidomain model for ionic electrodiffusion and osmosis with an application to cortical spreading depression. Phys. D 308, 94–108 (2015) ] published an elegant and comprehensive multidomain tissue-level model ionic electrodiffusion and osmotic water flow in biological tissue. In this framework, the tissue domain Ω\Omega is composed on RR coexisting compartments (for instance neuronal, glial, and extracellular spaces) with K\lvert K\rvert interacting ionic species (such as sodium Na, potassium K, chloride Cl). The coupled, time-dependent, nonlinear system of PDEs describes the volume fractions αr\alpha_{r}, ion concentrations [k]r[k]_{r} for kKk\in K, electric potentials φr\varphi_{r}, mechanical pressures prp_{r}, and fluid velocities uru_{r} for each compartment r=1,,Rr=1,\dots,R (see [26 Y. Mori, A multidomain model for ionic electrodiffusion and osmosis with an application to cortical spreading depression. Phys. D 308, 94–108 (2015) , 13 A. J. Ellingsrud, N. Boullé, P. E. Farrell and M. E. Rognes, Accurate numerical simulation of electrodiffusion and water movement in brain tissue. Math. Med. Biol. 38, 516–551 (2021) ]). For instance, in compartments r=1,,R1r=1,\dots,R-1, the movement of each ion concentration [k]r[k]_{r} is governed by t(αr[k]r)+divJrkJrRk\partial_{t}(\alpha_{r}[k]_{r})+\operatorname{div}J_{r}^{k}\propto J_{rR}^{k}, where JrkJ_{r}^{k} is the ion flux density and JrRkJ_{rR}^{k} is the transmembrane water flux between the compartments rr and RR. Coupling (chemical) diffusion, drift due to the electric field, as well as fluid flow within each compartment, yields ion fluxes of the form

Jrk=Drk[k]rDrkzkψ[k]rgradφr+αrur[k]r;J_{r}^{k}=-D_{r}^{k}[k]_{r}-\frac{D_{r}^{k}z^{k}}{\psi}[k]_{r}\operatorname{grad}\varphi_{r}+\alpha_{r}u_{r}[k]_{r};

see the aforementioned references for a complete description.

Figure 7. Accurate methods for simulating the interplay between chemical, electrical and mechanical quantities of interest (here during wave propagation through a 10 mm line of brain tissue) is an essential foundation for studying computationally the role of the brain’s waterscape in brain signalling and vice versa [13].

The combination of diffusion, nonlinear reaction and convection dynamics makes this a highly interesting system to approximate numerically using, e.g., finite element or spectral methods [13 A. J. Ellingsrud, N. Boullé, P. E. Farrell and M. E. Rognes, Accurate numerical simulation of electrodiffusion and water movement in brain tissue. Math. Med. Biol. 38, 516–551 (2021) ]. Moreover, the underlying physiology induces sharp wave fronts that require a very fine spatial and fine temporal resolution for all current methods (Figure 7). As such, the numerical simulation of ionic electrodiffusion and water movement remains a real challenge, and solution approaches that retain accuracy at lower computational expense could enable further uptake in the community.

Remark 14. As the Waterscales project nears its end, it has been rewarding to observe its trajectory from roots in numerical analysis, rapid impact of poroelasticity and fluid dynamics in brain mechanics, and stretching toward impact in neuroscience in the long run. Our most recent work [15 A. J. Ellingsrud, D. B. Dukefoss, R. Enger, G. Halnes, K. Pettersen and M. E. Rognes, Validating a computational framework for ionic electrodiffusion with cortical spreading depression as a case study. eNeuro 9, DOI 10.1523/ENEURO.0408-21.2022 (2022) ] took us back to the project’s origins, by realizing a collaboration with medical experimentalists at the GliaLab at the University of Oslo involved in the development of the original glymphatic theory, while its publication in eNeuro points towards a closer relation with the neuroscience community.

Remark 15. Another question is whether the tissue level is the appropriate scale for modelling the interplay between electrical, chemical and mechanical effects. Indeed, we can envision a paradigm for modelling excitable tissue at the level of cells instead, with explicit representation of cell membranes and spaces. This is however a story for another day [35 A. Tveito, K.-A. Mardal and M. E. Rognes, Modeling excitable tissue: The EMI framework. Simula SpringerBriefs Comput., Rep. Comput. Physiol. 7, Springer, Cham (2021) , 14 A. J. Ellingsrud, C. Daversin-Catty and M. E. Rognes, A cell-based model for ionic electrodiffusion in excitable tissue. In Modeling excitable tissue: The EMI framework, edited by A. Tveito, K.-A. Mardal and M. E. Rognes, Simula SpringerBriefs Comput., Rep. Comput. Physiol. 7, Springer, Cham, 14–27 (2021) ].

6 Computational abstractions and algorithms

Over the last decades, we have witnessed a tremendous increase in the physical, mathematical, numerical and computational complexity of models for physiological processes. To keep up, there has been a corresponding change in how algorithms and software frameworks for the numerical solution of (partial) differential equations are designed. In the context of finite element methods, the FEniCS Project provides a computational platform combining high-level abstractions for the specification of discrete variational formulations with automated code generation of optimized kernels [2 M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes and G. N. Wells, The FEniCS project version 1.5. Archive of Numerical Software 3, 9–23 (2015) , 18 P. E. Farrell, D. A. Ham, S. W. Funke and M. E. Rognes, Automated derivation of the adjoint of high-level transient finite element programs. SIAM J. Sci. Comput. 35, C369–C393 (2013) ]. However, FEniCS support for multidomain, multiscale and/or multiphysics problems has been sorely lacking.

A core Waterscales aim was to remedy this situation by designing numerical algorithms and software abstractions that allow for high-level specification and high-performance forward and reverse solution of models with multiscale features. We addressed this goal by designing and introducing mathematical software concepts together with lower-level algorithms for expressing, representing, and solving systems of PDEs coupled across interfaces or subdomains (Figure 8) [10 C. Daversin-Catty, C. N. Richardson, A. J. Ellingsrud and M. E. Rognes, Abstractions and automated algorithms for mixed domain finite element methods. ACM Trans. Math. Software 47, Article No. 31, DOI 10.1145/3471138 (2021) ]. These tools enable automated assembly and solution of a wide range of mixed finite element variational formulations, such as, e.g., the finite element spaces and formulations involved in the reduced perivascular flow models (11), interactions across the cell membrane in geometrically resolved models of excitable tissue [35 A. Tveito, K.-A. Mardal and M. E. Rognes, Modeling excitable tissue: The EMI framework. Simula SpringerBriefs Comput., Rep. Comput. Physiol. 7, Springer, Cham (2021) , 16 A. J. Ellingsrud, A. Solbrå, G. T. Einevoll, G. Halnes and M. E. Rognes, Finite element simulation of ionic electrodiffusion in cellular geometries. Front. Neuroinf., DOI 10.3389/fninf.2020.00011 (2020) ] or fluid-structure interfaces [6 M. Causemann, V. Vinje and M. E. Rognes, Human intracranial pulsatility during the cardiac cycle: A computational modelling framework. Fluids Barriers CNS 19, 863–874 (2022) ]. All algorithms are publicly and openly available via the FEniCS Project software [2 M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes and G. N. Wells, The FEniCS project version 1.5. Archive of Numerical Software 3, 9–23 (2015) , 10 C. Daversin-Catty, C. N. Richardson, A. J. Ellingsrud and M. E. Rognes, Abstractions and automated algorithms for mixed domain finite element methods. ACM Trans. Math. Software 47, Article No. 31, DOI 10.1145/3471138 (2021) ].

Figure 8. Computational representations of mixed-domain and mixed-dimensional problems take advantage of mappings within (simplicial) complexes such as meshes and submeshes, faces and stars [10, Figure 4].

Open problem 16. High-level computational specification and automated solution of, e.g., coupled n×(n2)n\times(n-2)-dimensional problems (n2n\geqslant 2), such as to represent interactions between the (peri)vasculature and CSF or tissue, comes with an additional set of challenges including, e.g., non-conforming averaging and extension operators and computational geometry classification problems, and remains an open problem.

7 Brain clearance and neurodegeneration

All-in-all, why are we really so interested in solute transport in the human brain? Well, there are two intertwined clinical reasons – either to transport treatment drugs into the brain, or to understand how metabolic waste clears from the brain – in order to better comprehend and hopefully treat neurodegenerative diseases such as, e.g., Alzheimer’s disease.

Figure 9. Stages of toxic protein concentrations (coronal view upper row, sagittal view lower row) governed by (12). Courtesy of G. Brennan, Mathematical Institute, University of Oxford.

With new and quantitative physiological insights at hand, we could now mathematically describe the role of dynamical clearance in the propagation of toxic proteins across the brain. Representing the brain’s connectome (i.e., fiber bundle pathways) by a graph (V,E)(\mathcal{V},\mathcal{E}), we study the distribution of the protein concentrations pip_{i} and regional clearance values λi\lambda_{i} (iVi\in\mathcal{V}), evolving via the graph Laplacian LL that describes anisotropic diffusion and nonlinear reactions as

tpi=ρj=1NLijpj+(λcritλi)piαpi2,\displaystyle\partial_{t}p_{i}=-\rho\sum_{j=1}^{N}L_{ij}p_{j}+(\lambda_{\mathrm{crit}}-\lambda_{i})p_{i}-\alpha p_{i}^{2},
tλi=βipi(λiλi),\displaystyle\partial_{t}\lambda_{i}=-\beta_{i}p_{i}(\lambda_{i}-\lambda_{i}^{\infty}),

with initial conditions on pi,λip_{i},\lambda_{i}; here αi,βi\alpha_{i},\beta_{i} are kinetic parameters, and λicrit\lambda_{i}^{\mathrm{crit}} and λi\lambda_{i}^{\infty} critical and minimal clearance values, respectively [4