Geometrically Reduced Modelling of Pulsatile Flow in Perivascular Networks

Flow of cerebrospinal fluid in perivascular spaces is a key mechanism underlying brain transport and clearance. In this paper, we present a mathematical and numerical formalism for reduced models of pulsatile viscous fluid flow in networks of generalized annular cylinders. We apply this framework to study cerebrospinal fluid flow in perivascular spaces induced by pressure differences, cardiac pulse wave-induced vascular wall motion and vasomotion. The reduced models provide approximations of the cross-section average pressure and cross-section flux, both defined over the topologically one-dimensional centerlines of the network geometry. Comparing the full and reduced model predictions, we find that the reduced models capture pulsatile flow characteristics and provide accurate pressure and flux predictions across the range of idealized and image-based scenarios investigated—at a fraction of the computational cost of the corresponding full models. The framework presented thus provides a robust and effective computational approach for large scale in-silico studies of pulsatile perivascular fluid flow and transport.

1 Introduction

Flow of cerebrospinal fluid (CSF) in perivascular spaces (PVSs) is a key transport mechanism in and around the brain [1–3]. A PVS is a space or potential space along or around a blood vessel through which fluid and particles can pass [4]. Such spaces appear along blood vessels on the brain surface (surface or pial PVSs) or along blood vessels within the brain parenchyma (parenchymal PVSs). While their shape and structure, and to some extent existence, remain disputed [4–8], PVSs are typically represented as (elliptic) annular structures surrounding the blood vessels. As such, surface and parenchymal PVSs form structural networks, dual to and in close interaction with the vascular network, and the surrounding brain tissue and/or subarachnoid space.

Mathematical and computational models are playing an increasingly important role in understanding and predicting PVS flow characteristics [9]. Theoretical models have quantified the resistance in PVS networks [10], while detailed numerical simulations can predict perivascular fluid velocities and pressures in idealized [11–17] and image-based geometries [18]. However, computational fluid dynamics simulations rapidly become prohibitively expensive for large, three-dimensional PVS networks. A natural question is therefore whether reduced models can accurately capture PVS flow and transport characteristics and magnitudes. Of particular interest and relevance are geometrically-reduced models for which the computational domain is reduced from an initial three-dimensional representation to a network of topologically one-dimensional branches. Such models have been subject to active research over the last decades in the context of the vasculature, arterial blood flow, and tissue perfusion [19–29]. For the one-dimensional arterial blood flow models, see e.g., the seminal work of Olufsen [19], the vasculature is typically represented by a branching network of centerlines, and the model variables are the time-varying cross-section flux and vascular area. The corresponding PVS flow setting has received less attention from the mathematical and numerical community on the other hand.

In this work, we introduce a geometrically-reduced mathematical model and numerical solution techniques for the time-dependent flow of an incompressible viscous fluid such as CSF in surface PVS networks. The cross-section flux and average pressure are the primary model variables. We consider different computational scenarios including PVS flow induced by a systemic pressure gradient, by cardiac pulse wave-induced movement of the inner vascular wall and by vasomotion in idealized or image-based model geometries. We evaluate the accuracy and efficiency of the reduced models by qualitative and quantitative comparison with the full three-dimensional model analogues.

The reduced models provide accurate approximations of the cross-section average pressure, cross-section flux and net flow in all geometries considered with relative model discrepancies in the peak flux between 0 and 35% and in the peak pressure between 0 and 52%. For realistic three-dimensional geometries, the reduced model reduces the computational costs (memory and runtime) by factors of 50 −200× with higher factors expected for larger scale networks.

2 Materials and Methods

2.1 PVS Geometries (3D and 1D)

In general, we consider a perivascular tree-like domain Ω consisting of a network of branching generalized annular cylinders Ωⁱ, with Ω ⊆∪_i∈IΩⁱ, spatial coordinates x ∈ Ω and time t ≥ 0. The boundary is denoted ∂Ω, with boundary normal n. We assume that each generalized annular cylinder Ωⁱ has a well-defined and oriented, topologically one-dimensional centerline Λⁱ with coordinate s. We set Λ = ∪_i∈IΛⁱ. Along s, we define the cross-sections Cⁱ = Cⁱ(s, t) of Λⁱ with area Aⁱ = Aⁱ (s, t). We denote the inner radius of Ωⁱ by $R_1^i$ and the outer radius of Ωⁱ by $R_2^i$; these radii will in practice vary with s, t and the angular coordinate θ. We denote the set of bifurcation points i.e. the points at which the centerlines of branches meet by $\mathcal{B}$.

We introduce three specific geometries of increasing complexity: from an axisymmetric cylinder (A) to an image-based perivascular geometry without any bifurcations (B) and one with a bifurcation (C) (Figure 1 and Table 1). The image-based geometries (B) and (C) are constructed from non-pathological artery segments from the Aneurisk dataset repository [30], and thus define high-fidelity 3D representations of human brain surface arteries. In each of these geometries, the PVS domain is defined by creating a generalized annular cylinder surrounding the vascular segment with the vascular wall as the inner surface of the PVS. The width of the PVS is set proportional to the blood vessel diameter (by factor of 0.95) and scaled (to a mouse scale) [18, 31]. Three-dimensional PVS flow in geometries A and C have been studied previously [18] and will be used for comparison. We define as PVS inlets and outlets (∂Ω_in and ∂Ω_out) the PVS ends surrounding the vascular inlets and outlets, respectively, noting however that fluid may flow both in and out of both the inlet and outlets. We denote the inner PVS wall (boundary) by ∂Ω_inner and outer wall by ∂Ω_outer.

FIGURE 1

FIGURE 1. Overview of the full three-dimensional and topologically one-dimensional reduced model domains. The idealized geometry (A) (the axisymmetric PVS) is a single 1 mm long axisymmetric annular cylinder represented by its two-dimensional angular cross-section. Geometry (B) (the image-based PVS) is generated from a cerebral artery segment (Aneurisk dataset repository, case id C0092) and represents an image-based perivascular space without bifurcation. Geometry (C) (the bifurcating image-based PVS) is generated from a middle cerebral artery (MCA M1–M2) segment (Aneurisk dataset repository, case id C0075) and represents an image-based perivascular space including a bifurcation.

TABLE 1

TABLE 1. Geometrical or numerical PVS domain characteristics for domains A, B, C.

The 3D PVS construction and the 1D centerline extraction are performed using PVS-meshing-tools [32], largely based on VMTK [33]. The extracted centerline comes with underlying data including the branch lengths and vessel radii. The centerline radius refers to the radius of the maximal inscribed circle of the vessel cross-sections. The meshing of both 3D and 1D PVS domains is performed within PVS-meshing-tools [32] using meshio [34] and GMSH [35]. The centerline meshes consist of topologically one-dimensional intervals embedded in three dimensions. The bifurcation points $b\in \mathcal{B}\subset \mathrm{\Omega }$ are explicitly labeled within each centerline mesh. Each branch is also separately tagged and given a consistent orientation. This procedure allows for the identification of bifurcation points as the outlet of one (parent) centerline and the inlet of other (daughter) centerlines, and a split of the full perivascular network into oriented mesh branches.

2.2 Stokes Flow in a Deforming Perivascular Domain

Flow of CSF in surface PVSs is reported to be laminar, with low Reynolds numbers (10^–4−10^–2) and moderate Péclet numbers (10²–10⁴) for 1 μm spherical particles transported at low Reynolds number [31], a mean flow speed of up to 60 μm/s, and parabolic flow profiles [31]. We therefore model the flow of an incompressible, viscous fluid flowing at low Reynolds and Womersley numbers via the time-dependent Stokes equations over a time-dependent domain Ω = Ω(t) representing the PVS. The fluid velocity v = v (x, t) for x ∈ Ω(t) at time t and the CSF pressure p = p (x, t) then solve the following system of time-dependent partial differential equations (PDEs) [18, 36]:

$\rho {\partial }_{t}v-\mu {\nabla }^{2}v+\nabla p=0\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\left(t\right),(\mathrm{1a})$

$\nabla \cdot v=0\text{\hspace{0.17em}in\hspace{0.17em}}\mathrm{\Omega }\left(t\right),(\mathrm{1b})$

where ρ is the fluid density and μ is the dynamic fluid viscosity. To model CSF at body temperature, we set the fluid density to ρ = 10³ kg/m³ and the dynamic viscosity to μ = 0.697 × 10^–3 Pa s. As in our previous full models of perivascular flow [18], the initial PVS mesh defines the reference domain Ω(0), and we assume that Ω(t) at time t > 0 is given by a deformation d of the reference domain: Ω(0)↦Ω(t) with x = d (X, t), X ∈ Ω(0), x ∈ Ω(t). We denote the domain velocity associated with d by w (thus $\dot{d}=w$).

2.3 Boundary Conditions, Initial Conditions and Periodicity

At the PVS ends, we prescribe a traction condition corresponding to a known, applied pressure $\tilde{p}=\tilde{p}(x,t)$:

${\sigma }_n\equiv \left(\mu \nabla u-pI\right)\cdot n=-\tilde{p}n\text{\hspace{0.17em}on\hspace{0.17em}}\partial {\mathrm{\Omega }}_{\text{in}}\text{\hspace{0.17em}and\hspace{0.17em}}\partial {\mathrm{\Omega }}_{\text{out}}.(2)$

We either prescribe (i) zero pressure at both ends $\tilde{p}=0$, or (ii) a constant-in-time pressure gradient $\mathrm{\Delta }\tilde{p}>0$ by setting ${\tilde{p}}_{\text{in}}=L_{\text{out}}\mathrm{\Delta }\tilde{p}$ at the inlet, letting ${\tilde{p}}_{\text{out}}=0$ at the outlet furthest from the inlet with distance L_out, and setting ${\tilde{p}}_{\text{out}}$ at any other outlets such that the average pressure gradient over each branch path $({\tilde{p}}_i-{\tilde{p}}_{\text{out}})/L_{\text{out}}$ is constant and equal to the prescribed pressure gradient $\mathrm{\Delta }\tilde{p}$ mmHg/m. This static pressure difference can represent e.g., a hydrostatic pressure difference, a venous pressure differential, or some other systemic pressure difference. Other types of boundary conditions could also be considered, see e.g., a discussion of compliance conditions in [18], or [37].

On the inner and outer PVS walls (along the length of the PVS), we set the fluid velocity v to match a known, prescribed domain velocity w = w (x, t). For the inner PVS wall, we either (i) consider a rigid wall and set v = w = 0, or (ii) impose a pulsating wall displacement:

$d{|}_{\partial {\mathrm{\Omega }}_{\text{inner}}}\left(X,t\right)=A\left(X,t\right)n,(3)$

with reference to the initial (fixed) mesh with coordinates X and prescribe $v=w=\dot{d}$. To represent wall motion induced by the cardiac pulse wave, we let the amplitude A be defined by the juxtaposition of an experimentally-observed wall motion time series [31] either applied uniformly along the length of the PVS or as a travelling wave along the PVS length with wave speed c = 1 m/s and frequency 10 Hz. We refer to [18] for the detailed description. To represent wall motion due to vasomotion, we consider a similar set-up but with a travelling sinusoidal wave in time with a frequency of 0.1 Hz and wave length λ = 8 mm [38], and an amplitude A of 7.5% of the initial inner radius R₁. We note that for all models, the wall moves in the normal (radial) direction only. For the outer PVS wall ∂Ω_outer, we set v = w = 0.

The system starts at rest with v = w = 0 at t = 0. The system reaches the periodic steady state nearly immediately, and we report results starting from the first cycle.

2.4 Model Reduction Assumptions

We define a reduced, topologically one-dimensional, model approximation of the full PVS flow model [(1) with the given boundary and initial conditions] under the following stipulations. For each branch Ωⁱ(t) with centerline Λⁱ and local coordinate system (s, r, θ), where s represents the path length (or axial coordinate), r is the radial coordinate and θ is the angular coordinate, we suppose that:

(I) Axial symmetry. Fields and input parameters are independent of the angular coordinate θ;

(II) Radial displacements. Boundaries displace in the radial direction only;

(III) Fixed centerline. The centerline Λ is fixed in time and defines the axial direction;

(IV) Constant cross-section pressure. The pressure field is independent of the angular and radial coordinates i.e., p = p (s, t);

(V) Axial velocity profile The axial velocity v_s, i.e., the velocity component in the axial direction can be decomposed in the form

$v\cdot s=v_s=v_s\left(s,r,t\right)=\widehat{v}\left(s,t\right)v_{\text{vp}}\left(r\right),(4)$

where v_vp is a given velocity profile varying radially only, $\widehat{v}$ is to be determined.

For the velocity profile v_vp, we here choose a normalized annular Poiseuille flow:

$v_{\text{vp}}\left(r\right)=\frac{v_{\text{poise}}\left(r\right)}{v_{\text{poise}}\left(\frac{R_1+R_2}{2}\right)},v_{\text{poise}}\left(r\right)=\left(1-\frac{r^2}{R_1^2}+\frac{R_2^2-R_1^2}{R_1^2\ln \left(R_2/R_1\right)}\ln \left(r/R_1\right)\right).(5)$

This velocity profile is parabolic in r (as for Poiseuille flow in a cylinder) with a logarithmic correction that accounts for the annulus.

In particular, the domain velocity w is assumed independent of the angular coordinate θ. Note that we do not assume other velocity components (than the axial) to necessarily be zero. We emphasize that these assumptions will in general not be satisfied by realistic geometries and flows. Thus, the reduced model defines a model approximation associated with a certain modelling error.

2.5 Reduced Model Equations

Under the assumptions (I–V), the full PVS flow model can be reduced to the following system of time-dependent differential equations: find the cross-section flux $\widehat{q}=\widehat{q}(s,t)$ and the cross-section average pressure $\widehat{p}=\widehat{p}(s,t)$ such that for each centerline Λⁱ (denoting $\widehat{q}{|}_{{\mathrm{\Lambda }}^i}={\widehat{q}}^i$ and $\widehat{p}{|}_{{\mathrm{\Lambda }}^i}={\widehat{p}}^i$):

$\frac{\rho }{A^i}{\partial }_t{\widehat{q}}^i-\frac{\mu }{A^i}{\partial }_{ss}{\widehat{q}}^i+\mu \frac{{\alpha }^i}{A^i}{\widehat{q}}^i+{\partial }_s{\widehat{p}}^i=0\text{\hspace{0.17em}on\hspace{0.17em}}{\mathrm{\Lambda }}^i,(\mathrm{6a})$

${\partial }_s{\widehat{q}}^i={\widehat{f}}^i\text{\hspace{0.17em}on\hspace{0.17em}}{\mathrm{\Lambda }}^i,(\mathrm{6b})$

hold.

${\widehat{f}}^i\left(s\right)\equiv 2\pi R_1^i\left(s,t\right)w\left(R_1,s,t\right)\cdot n{|}_{\partial {\mathrm{\Omega }}_{\text{inner}}}+2\pi R_2^i\left(s,t\right)w\left(R_2,s,t\right)\cdot n{|}_{\partial {\mathrm{\Omega }}_{\text{outer}}}.(7)$

Moreover, Aⁱ = Aⁱ (s, t) denotes the cross-section area, while ${\widehat{\alpha }}^i={\widehat{\alpha }}^i(s,t)$ is a lumped flow parameter that depends on the domain geometry and the choice of velocity profile v_vp:

${\alpha }^i\left(s,t\right)\equiv \frac{1}{A{\overline{\stackrel{̄}{v}}}_{\text{vp}}\left(s\right)}\left(2\pi R_1^i\left(s,t\right){\partial }_r{v}_{\text{vp}}\left(R_1^i\left(s,t\right)\right)-2\pi R_2^i\left(s,t\right){\partial }_r{v}_{\text{vp}}\left(R_2^i\left(s,t\right)\right)\right),(8)$

and where ${\overline{\stackrel{̄}{v}}}_{\text{vp}}$ is the velocity profile integrated over each cross-section:

${\overline{\stackrel{̄}{v}}}_{\text{vp}}\equiv {\int }_{C\left(s\right)}{v}_{\text{vp}}r\mathrm{d}r\mathrm{d}\theta .(9)$

We also define the (one-dimensional) normal stress induced by $\widehat{q}$ and $\widehat{p}$:

$\widehat{\sigma }\equiv \frac{\mu }{A}{\partial }_s\widehat{q}-\widehat{p},(10)$

which corresponds to an average of the axial (s-)component of the normal stress in (2) over each cross-section

At the bifurcation points $b\in \mathcal{B}\subset \mathrm{\Omega }$, we impose the following two conditions representing conservation of flux and continuity of normal stress, respectively:

${\widehat{q}}^p\left(s^p\right)={\widehat{q}}^{d_1}\left(s^{d_1}\right)+{\widehat{q}}^{d_2}\left(s^{d_2}\right),(11)$

${\widehat{\sigma }}^p\left(s^p\right)={\widehat{\sigma }}^{d_1}\left(s^{d_1}\right)={\widehat{\sigma }}^{d_2}\left(s^{d_2}\right),(12)$

where Λ^p and ${\mathrm{\Lambda }}^{d_1}$, ${\mathrm{\Lambda }}^{d_2}$ represent the centerlines of the parent and two daughter branches, respectively, associated with the bifurcation point b and s^⋅ = ι^⋅(b) where ι^⋅ denotes the map from three-dimensional bifurcation point to the one-dimensional centerline coordinate for each branch Ω^⋅.

System (6) defines a set of equations for each branch centerline Λ_i and is closed by the bifurcation conditions (11, 12), together with boundary conditions at the PVS inlet and outlets, as well as initial conditions for the cross-section flux. Specifically, in place of the traction condition (2), we prescribe the corresponding pressure difference for the (average) normal stress $\widehat{\sigma }$ cf. (10). In this manner, the (one-dimensional) solutions ${\widehat{q}}^i$ and ${\widehat{p}}^i$ of the reduced model (6) define approximations of the (three-dimensional) axial flux and pressure solving (1) integrated or averaged over each cross-section:

$\begin{array}{ll}\hfill {\widehat{q}}^i\left(s\right)& \approx {\int }_{C^i\left(s\right)}{v}_s\left(s,r,t\right)r\mathrm{d}r\mathrm{d}\theta \equiv A^i\left(s\right){\bar{q}}^i\left(s\right),\hfill \\ \hfill {\widehat{p}}^i\left(s\right)& \approx \frac{1}{A^i\left(s\right)}{\int }_{C^i\left(s\right)}p\left(s,r,t\right)r\mathrm{d}r\mathrm{d}\theta .\hfill \end{array}$

The factor r originates from integrating in cylindrical coordinates. We note that the wall velocity w, which defines a boundary condition for the full PVS model (1), enters as a body force in the reduced model (6).

2.6 Numerical Solution and Software

We solve the full PVS (1) via a previously developed and verified arbitrary Lagrangian-Eulerian (ALE) formulation and finite element discretization [18]. This solver builds on the standard FEniCS finite element software suite [39], and is openly available [40].

To compute numerical solutions to the reduced model (6), we consider a first-order implicit Euler scheme in time and a higher-order finite element method in space. The finite element mesh $\mathcal{T}$ of the centerline Λ is composed of mesh segments ${\mathcal{T}}^i$, one for each centerline branch Λⁱ. Each mesh segment is a mesh consisting of intervals embedded in ${\mathbb{R}}^3$. We label the set of bifurcation points $\mathcal{B}$, inlet points $\mathcal{I}$ and outlet points $\mathcal{O}$, and define the following finite element spaces:

• The flux space V_h is the space of continuous piecewise quadratics over ${\mathcal{T}}^i$ for each i.

• The (average) pressure space Q_h is the space of continuous piecewise linears on $\mathcal{T}$.

• The Lagrange multiplier space $R_h={\mathbb{R}}^B$ where B is the number of bifurcation points.

The flux is thus solved on each mesh segment representing the PVS network branches and may be discontinuous across bifurcations. We impose the flux conservation condition (11) weakly using a Lagrange multiplier formulation. The pressure is solved on the whole mesh and is continuous at bifurcations by construction.

For each discrete time t^k, given ${\widehat{q}}_h^{k-1}$ at the previous time t^k−1 and time step Δt = t^k − t^k−1, we solve for the approximate cross-section flux ${\widehat{q}}_h^k\in V_h$, average pressure ${\widehat{p}}_h^k\in Q_h$ and a Lagrange multiplier [corresponding to the normal stress (10) at the bifurcation points] ${\lambda }_h^k\in R_h$ solving

$a\left(\left({\widehat{q}}_h^k,{\widehat{p}}_h^k,{\lambda }_h^k\right),\left(\psi ,\varphi ,\xi \right)\right)=L^k\left(\left(\psi ,\varphi ,\xi \right)\right),(13)$

for all finite element test functions ψ ∈ V_h, ϕ ∈ Q_h, and ξ ∈ R_h. The left-hand side bilinear form a is defined by:

$\begin{array}{l}a\left(\left(q,p,\lambda \right),\left(\psi ,\varphi ,\xi \right)\right)=\hfill \\ \sum _{i\in I}{\int }_{{\mathrm{\Lambda }}^i}\frac{1}{A^i}\left(\rho +\mathrm{\Delta }t\mu {\alpha }^i\right)q^i{\psi }^i+\frac{\mathrm{\Delta }t\mu }{A^i}{\partial }_s{q}^i{\partial }_s{\psi }^i+{\partial }_s{q}^i{\varphi }^i-\mathrm{\Delta }t{\partial }_s{\psi }^i{p}^i\mathrm{d}s+\sum _{b\in \mathcal{B}}{\lambda }^b{\left[\psi \right]}^b+{\xi }^b{\left[q\right]}^b,\hfill \end{array}(14)$

where λ^b (or ξ^b) is simply the entry of the vector λ (or ξ) corresponding to bifurcation point b, and we define the natural jump:

${\left[\psi \right]}^b={\psi }^p\left(b\right)-{\psi }^{d_1}\left(b\right)-{\psi }^{d_2}\left(b\right).(15)$

The right-hand side linear form L is:

$L^k\left(\psi ,\varphi ,\xi \right)=\sum _{i\in I}{\int }_{{\mathrm{\Lambda }}^i}\frac{\rho }{A^i}{\widehat{q}}_h^{k-1,i}{\psi }^i+f^i{\varphi }^i\mathrm{d}s-\sum _{x\in \mathcal{I}}\mathrm{\Delta }t{\tilde{p}}_{\text{in}}\left(x\right){\psi }^{i_I}\left(x\right)+\sum _{x\in \mathcal{O}}\mathrm{\Delta }t{\tilde{p}}_{\text{out}}\left(x\right){\psi }^{i_O}\left(x\right),(16)$

where the superscript i_I (i_O) in the inlet (outlet) terms above refers to the unique centerline branch associated with the inlet (outlet) points.

The numerical solver for the reduced model was implemented in the well-established FEniCS Project finite element software [39]. The solver, and in particular the definition of the partially continuous flux space, builds on mixed-domain features [41] and relies on the latest development version of FEniCS. All data and source code are available via Zenodo [42].

2.7 Overview of Computational Models, Output Functionals and Model Error Measures

An overview of the six computational models considered is given in Table 2. Each model is labeled with reference to its domain (A, B, or C) followed by a number indicating the driving forces included: (1) a given pressure drop, (2) wall movement due to cardiac pulsations and (3) wall movement due to vasomotion. For each model, we consider the full three-dimensional version as well as the reduced model.

TABLE 2

TABLE 2. Overview of computational models parameterized by domain, prescribed pressure gradient $\mathrm{\Delta }\tilde{p}$ and wall motion pattern (see Methods).

To compare the solutions from the full and reduced models, we consider the following quantities of interest. For each domain, we define a set of cross-sections as follows. For domain A, we define the left-most end as the inlet (s = 0) and define an upper cross-section. For domain B, we consider the inlet and outlet ends of the PVS, as well as upper and lower cross-sections. For domain C, we consider the inlet at s = 0, and the two outlets, as well as three additional cross sections near the inlet, on the largest daughter branch relatively close to the bifurcation, and near the outlet of the other daughter branch.

We then compute for each cross section C(s) the averaged pressure ${\overline{\stackrel{̄}{p}}}_h(s,t)$ and the cross-section flux ${\overline{\stackrel{̄}{q}}}_h(s,t)$:

${\overline{\stackrel{̄}{p}}}_h\left(s,t\right)=\frac{1}{|C|}\sum _k{w}_k{p}_h\left(x_k,t\right),{\overline{\stackrel{̄}{q}}}_h\left(s,t\right)=\sum _k{w}_k{v}_h\left(x_k,t\right)\cdot n_C\left(x_k\right)(17)$

for a quadrature scheme with points x_k and weights w_k defined over C and an approximation |C| of the cross-section area. Here n_C is the normal vector of the cross-section. The averaging is implemented by using the Frenet frame associated with Λ to map from an annular cylinder in a reference domain onto the cross-section, similar to the implementation of the averaging operator in fenics_ii [43].

With this in hand we define the percentagewise relative model discrepancy E_q(t) in the flux by

$E_q\left(t\right)=\frac{\Vert {\overline{\stackrel{̄}{q}}}_h\left(t\right)-{\widehat{q}}_h\left(t\right){\Vert }_{L^2\left(\mathrm{\Lambda }\right)}}{\Vert {\overline{\stackrel{̄}{q}}}_h\left(t\right){\Vert }_{L^2\left(\mathrm{\Lambda }\right)}}\cdot 100\%(18)$

and similarly for the pressure E_p. We typically compute this quantity if the flow is driven by a constant pressure gradient. In this case the fluid starts at rests and then quickly develops to stationary, annular Poiseuille flow. We then compute E_q(T) and E_p(T), where T denotes the final time.

For pulsatile flow, we typically compare the percentagewise relative error in peak pressure $e_p^{\text{peak}}(s^{\prime })$ and peak cross-section flow $e_q^{\text{peak}}(s^{\prime })$ at some cross section C (s′, t), where

$e_q^{\text{peak}}\left(s^{\prime }\right)=\frac{|{\max }_{t\in \left(0,T\right)}{\widehat{q}}_h\left(s^{\prime },t\right)-{\max }_{t\in \left(0,T\right)}{\overline{\stackrel{̄}{q}}}_h\left(s^{\prime },t\right)|}{|{\max }_{t\in \left(0,T\right)}{\overline{\stackrel{̄}{q}}}_h\left(s^{\prime },t\right)|}\cdot 100\%(19)$

and $e_p^{\text{peak}}(s^{\prime })$ is similarly defined. Finally, we compare the net fluxes Q of the full and reduced model, where Q associated with the velocity v = v (x, t) can be computed as:

$Q={\int }_0^T{\int }_{\partial {\mathrm{\Omega }}_{\text{in}}}v\cdot n\mathrm{d}x\mathrm{d}t,(20)$

and the corresponding quantity associated with the flux $\widehat{q}=\widehat{q}(s)$ by:

$Q={\int }_0^T\sum _{x\in \mathcal{I}}\widehat{q}\left(x\right)\mathrm{d}t(21)$

where the integration in time is over one period [0, T].

3 Results

The prescribed pressure gradient and the pulsating PVS walls each induce pressure gradients and fluid flow in the different PVS geometries. For each of the models (Table 2), we compare the simulation results from the full PVS (1) defined over the three-dimensional model domains and the reduced system (6) defined over the topologically one-dimensional domains, quantify the discrepancies between the models and the computational costs.

3.1 Reduced Model Exactly Predicts Pressure-Driven Axisymmetric Flow Characteristics

Flow in an axisymmetric annular cylinder of length ℓ driven by a constant pressure difference Δp (Model A1) is described by the analytic expression:

$\begin{array}{ll}\hfill \widehat{q}\left(s,t\right)& =A\frac{\mathrm{\Delta }p}{\mu \alpha \ell }\left(1-\exp \left(-\frac{\mu \alpha t}{\rho }\right)\right),\hfill \\ \hfill \widehat{p}\left(s,t\right)& =\frac{\mathrm{\Delta }p}{\ell }s+\widehat{p}\left(0\right),\hfill \end{array}(22)$

where α is the lumped flow parameter given by (8) and which is constant in time and space in this case. For the velocity profile (5) defined over geometry A (cf. Table 1), α = 7325.3/m², and μα/ρ = 5105.7/s. Thus, the time-dependency is negligible after only a few milliseconds, and the flow develops near-instantaneously to steady-state Poiseuille flow.

Both the full and reduced models reproduce the exact annular Poiseuille flow characteristics of this case (Figure 2A). The numerical difference between the analytic and computed reduced solutions for the cross-section flux $\widehat{q}$ and average pressure $\widehat{p}$ is at machine precision ($\Vert \widehat{q}(T)-{\widehat{q}}_h(T)\Vert =1.7\times 10^{-14}$ and $\Vert \widehat{p}(T)-{\widehat{p}}_h(T)\Vert =2.6\times 10^{-17}$) (T = 1 s). In general, the total error is the sum of the model error and the numerical error associated with the space-time discrete approximation (13). For Model A1, the model error is zero as the model reduction assumptions (I–V) are exactly fulfilled by the geometry and flow pattern. As the total error also vanishes, we note that the numerical error is also negligible for this case.

FIGURE 2

FIGURE 2. PVS flux and pressure in an axisymmetric annular cylinder induced by a constant pressure difference or cardiac wall motion (Models A1, A2). (A) Model A1: A constant pressure gradient induces annular Poiseuille flow in both the full axisymmetric model (upper panel) and the reduced model (lower panel): snapshot of steady solution at T = 0.1. (B–E) Model A2: Inner wall pulsations induce bidirectional and oscillatory flow. (B) Snapshot of the full model solutions at peak outflux (t = 0.05). Different cross-sections are marked in green (at the inlet) and blue (in the interior). (C) Pressure (upper panel) and cross-section flux ${\overline{\stackrel{̄}{q}}}_h$ (lower panel). (D) Cross-section flux predicted by the full model (dotted line) and the reduced model (solid line) at inlet versus time. (D) As for (C) but at the interior cross-section marked in (B).

3.2 Reduced Model Accurately Captures Axisymmetric PVS Wall Pulsations

Next, we examine the PVS flow and pressure generated by uniform axisymmetric pulsations of the inner PVS wall (Model A2, Figures 2B–E). The inner wall movement changes the inner domain radius R₁ in time. The fluid is pushed out at the both ends as the PVS width decreases, and flows back in at both ends as the PVS width returns to baseline. This behaviour is reproduced by both the full (Figure 2B [18]) and reduced models (Figure 2C). We note that the reduced model assumptions (IV, V) do not hold in this scenario as the PVS axial velocity profile is no longer identical to the Poiseuille velocity profile, and the pressure is not perfectly constant on each cross-section. Comparing the full and reduced cross-section fluxes ${\overline{\stackrel{̄}{q}}}_h$ and ${\widehat{q}}_h$, we observe however that the two models still agree closely (Figures 2D,E), both at the inlet and at an interior cross-section. Moreover, the time-profile of the reduced and full cross-section flux approximations are very similar (both at the inlet and at the interior cross-section, Table), though with small (Δt s) shifts in time. The peak outfluxes occur at the inlet and outlet; the peak outflux for the full model is 1.53 × 10^–3 μm³/s, and 1.47 × 10^–3 μm³/s for the reduced model (Figure 2D). The peak pressure occurs in the middle of the domain; the peak pressure for the full model is 0.194 and 0.193 Pa for the reduced model. Using (19) the relative error in the peak cross-section flux at the inlet is $e_q^{\text{peak}}(0)=4.0\%$ and in the peak cross-section (average) pressure $e_p^{\text{peak}}(0)=0.8\%$. There is thus a small discrepancy between the two models, as expected by the violation of the reduced model assumptions.

3.3 Radial Geometry Variations Induce Small Model Errors

In contrast to the axisymmetric geometry A, the image-based geometries B and C express angular and axial variations in radius. The inner and outer radii of these geometries vary along the length of the domain (with s) and depend on the angular coordinate θ, with the latter violating model assumption I. To study the resulting model error in isolation, we again examine the pressure-driven flow predicted in full and reduced models but now of geometry B (Model B1, Figure 3). The full numerical approximation of the pressure is nearly constant over each cross-section. On the other hand, the velocity profile varies between cross-sections and with the angular coordinate within each cross-section (Figure 3A). Therefore, we expect a larger model error in the reduced model compared to the previous case(s). At steady state (t = 0.5), the reduced pressure approximation $\widehat{p}$ varies nearly linearly along the length of the domain as expected, and the reduced flux approximation $\widehat{q}$ is essentially constant along the centerline with value $\widehat{q}=4.28\times 10^{-4}\mu$L/s. Computing the corresponding cross-section flux from the full model, we find values ranging from 3.5 × 10^–4 to 5.31 × 10^–4 μL/s. The relative model discrepancy (18) in the pressure is E_p = 2.6% and for the flux E_q = 12.6%.

FIGURE 3

FIGURE 3. In an image-based perivascular segment with varying radii, a pressure difference between inlet and outlet induces a pressure field that is nearly constant on each cross-section, but a velocity field that varies with the radial, angular and axial coordinates. (A) Full pressure and velocity approximations in the domain (left) along with close-up views of the pressure (middle) and velocity magnitude (right) at two cross-sections; (B) Reduced average cross-section pressure (left) and cross-section flux approximations (right).

3.4 Reduced Model is Robust with Respect to Wall Motion Amplitude and Frequency

Cardiac wall motion and vasomotion may drive pulsatile perivascular flow with different flow characteristics. To evaluate the model discrepancy induced by different physiological drivers, we compare the full and reduced models over an image-based PVS segment driven by wall motion induced by the cardiac pulse wave (Model B2) and by vasomotion (Model B3). The cardiac pulse wave induces wall motion at a higher frequency (10 Hz) travelling at a higher wave speed (1000 mm/s), while vasomotion creates pulsations at lower frequencies (0.1 Hz) and at a lower wave speed (0.8 mm/s). Both models include angularly, axially and temporally varying radii, and we expect model assumptions I, IV-V to not hold.

Both pairs of models induce pulsatile bidirectional flow in and out of the PVS segment in synchrony with the pulsating wall (Figure 4, Supplementary Video S1) with peak pressure magnitude in the middle of the segment, and conversely, low velocities in the middle of the domain and higher velocities near the PVS ends. Both model scenarios lead to pressure fields that are nearly constant on each cross-section (Figures 4C, 5), but with angularly varying velocity profiles (Figures 4B, 5).

FIGURE 4

FIGURE 4. Cardiac wall motion induces substantial pulsatile pressures and velocities in an image-based perivascular space segment, with the reduced model accurately capturing flow, pressure and transport characteristics. (A) Snapshot of the pressure and velocity at time of peak pressure (t = 0.05 s); (B) Velocity at upper and lower cross-section [zoom of (A)]; (C) Pressure at upper and lower cross-sections [zoom of (A)]; (D) Cross-section flux from reduced model (left) and full model (right); (E) cross-section average pressure from reduced model (left) and full model (right); (F) full and reduced model cross-section fluxes at the lower cross-section over time; (G) full and reduced model pressures at the lower cross-section over time.

FIGURE 5

FIGURE 5. Vasomotion induces higher domain deformations but lower wall velocities, pressure differences and cross-section fluxes. (A) Snapshots of the full model pressure and velocity at different time points with cross-section velocities (top). (B) Average pressure (upper panel) and flux (lower panel) for the full and reduced models at upper and lower cross-sections over time. The values at the different cross-sections are slightly shifted in time due to the travelling vasomotion. The pressure model discrepancy dominates the flux differences.

For the cardiac wall motion, the overall cross-section average of the full pressure ${\overline{\stackrel{̄}{p}}}_h$ ranges from −0.05 to 0.26 Pa, while the full cross-section flux ${\overline{\stackrel{̄}{v}}}_h$ ranges from −1.54 × 10^–3 to 1.95 × 10^–3 μL/s. The reduced model accurately captures the temporal and spatial characteristics of the full model (Figures 4D–G). For the reduced model, the overall cross-section pressure ${\widehat{p}}_h$ ranges from −0.06 to 0.29 Pa, while the cross-section flux ${\widehat{q}}_h$ is between −1.61 × 10^–3 and 2.23 × 10^–3 μL/s. Comparing the full and reduced pressure and flux over time at an interior, lower cross-section with axial coordinate s′ (Figures 4B,C), we observe that the reduced model slightly overestimates the peak pressure and flux when compared to the full model (Figures 4F,G). Using (19) at this cross-section we find that the relative error of the peak pressure is $e_p^{\text{peak}}(s^{\prime })=19.0\%$ and the relative error of the peak flux is $e_q^{\text{peak}}(s^{\prime })=1.2\%$. One shall note that the space discretization of the initial 3D model has a non negligible impact on these relative errors. Indeed, using a finer 3D mesh composed of 333000 tetrahedrons instead of the initial 63000 lowers the relative error of the peak pressure $e_p^{\text{peak}}(s^{\prime })$ to 12.1%. The relative error on peak flux is not significantly impacted but $e_q^{\text{peak}}(s^{\prime })$ was already very small.

For the vasomotion scenario, the domain movement is larger compared to the cardiac wall motion, but the wall velocity is lower (peak wall speed of 0.001 vs. 0.005 mm/s). The resulting peak (in terms of magnitude) cross-section pressure is −0.012 Pa and peak cross-section flux is 9.14 × 10^–5 μL/s (Figure 5, Supplementary Video S2). These are one-to-two orders of magnitude lower than for the cardiac wall motion scenario. Comparing the full and reduced models in two interior (upper and lower) cross-sections, we observe that the cross-section pressure ${\widehat{q}}_h$ matches pulsatile behaviour of the average cross-section pressure in the full model ${\overline{\stackrel{̄}{q}}}_h$ (Figure 5B) but that the peak amplitude is higher. The largest model differences in pressure at the lower cross-section is at the peak pressure; there the relative difference in peak pressure is $e_p^{\text{peak}}(s^{\prime })=52.7\%$. The similar observations hold for the flux, but the model discrepancies are lower: the relative difference in peak flux is $e_q^{\text{peak}}(s^{\prime })=15.6\%$. Moreover, the full and reduced models agree on a pressure phase shift of 0.5s. In agreement with our previous findings, the reduced pressure approximation displays a greater model discrepancy with higher predicted pressure variations in the reduced model (Figure 5B).

3.5 Reduced Model Captures Flow and Transport Characteristics Through Bifurcations

Now, we turn to compare the full and reduced model predictions of physiologically realistic perivascular flow in an image-based PVS surrounding a vascular bifurcation (Model C12). The prescribed pressure difference between inlet and outlets as well as the cardiac wall motion induces pulsatile flow with a net flow component [18] (Figure 6A, Supplementary Video S3). We note that the domain radii vary both angularly and axially, also for the initial domain, and also that the presence of a bifurcation region induces non-Poiseuille/non-Womersley-type velocity profiles. Comparing the full and reduced average pressure and flux at the time of peak velocity (Figures 6F,J), we note that the reduced model captures the qualitative and quantitative flow and pressure characteristics. The bifurcation conditions are satisfied at the bifurcation point b (Figures 6C,G) with a parent branch flux $\widehat{q}(b){|}_{{\mathrm{\Lambda }}^P}=-5.52\times 10^{-4}\mu$L/s and daughter branch fluxes $\widehat{q}(b){|}_{{\mathrm{\Lambda }}^{d_1}}=-7.8\times 10^{-5}\mu$L/s and $\widehat{q}(b){|}_{{\mathrm{\Lambda }}^{d_2}}=-4.74\times 10^{-4}\mu$L/s. The uneven flux distribution is induced by the smaller average width of one of the daughter vessels. The predicted stress $\widehat{\sigma }$ is continuous (data not shown).

FIGURE 6

FIGURE 6. Flow through a bifurcating PVS (Model C12) (A) Snapshot of pressure and velocity from full model at peak velocity (t = 0.05). (B) Full versus reduced cross-section flux at inlet (in) and outlets (out1 and out2) over time. (C) Snapshot of reduced cross-section pressure at peak velocity. (D) Snapshot of average cross-section pressure at the same time. (E) Pressure at upper, middle and lower cross-sections [zoom of (A)]. (F) Full versus reduced cross-section average pressure at cross-sections. (G) Snapshot of reduced cross-section flux at peak velocity (t = 0.05). (H) Snapshot of cross-section flux from the full model at the same time. (I) Flux at upper, middle and lower cross-sections [zoom of (A)]. (J) Full versus reduced cross-section flux at cross-sections.

The reduced peak cross-section flux (over time) at the inlet is −1.5 × 10⁻³ μL/s, and 1.2 × 10⁻³ μL/s and 7.9, ×, 10⁻⁴ μL/s at the outlets (Figure 6B).

Comparing the relative difference in peak flux at the inlet and outlets, we note that the discrepancy is largest at larger daughter outlet (s_out) with a relative difference $e_p^{peak}(s_{out})=12.7\%$. Comparing the full and reduced peak pressures at the upper (s_u), middle (s_m) and lower (s_l) cross-sections, we find relative differences $e_p^{peak}(s_u)=1.1\%$, $e_p^{peak}(s_m)=4.1\%$, and $e_p^{peak}(s_l)=16.4\%$. The analogous numbers for the fluxes are $e_q^{peak}(s_u)=1.0\%$, $e_q^{peak}(s_m)=33.5\%$, and $e_q^{peak}(s_l)=5.1\%$. Thus, the relative differences in peak flux are larger near the bifurcation region (Figures 6F,J).

The net flow is a key quantity of interest for the physiological relevance of perivascular flow and transport. The net flow per cycle in the full model is 3.5 × 10⁻⁵, and 2.9 × 10⁻⁵ μL for the reduced model, corresponding to a relative difference of 17%.

3.6 Reduced Models Offer Orders of Magnitude Saving in Computational Resources

Accurate direct three-dimensional simulations of pulsatile perivascular fluid flow in large, deforming vascular networks involve a significant computational cost. The expense is dominated by solving large linear systems of equations at each time step. For instance, even the moderate-resolution single-bifurcation model considered here (model C12) includes more than 17 000 vertices, 88 000 mesh cells and 400,000 degrees of freedom. For a small-scale idealized model such as axisymmetric Model A2, the reduced model uses 2.1% of the number of degrees of freedom but approximately the same amount of memory and longer runtime (0.16 vs. 0.35 s per time step, Table 3). However, the one-dimensional models reduce computational cost substantially for the image-based geometries (Table 3). For the image-based perivascular segment (Model B2), the reduced model uses 0.4% of the number of degrees of freedom, 2.0% of the runtime, and 2.1% of the memory of the full model. For the image-based bifurcating PVS (Model C12), the reduced model uses 0.18% of the number of degrees of freedom, 0.6% of the runtime and 2.0% of the memory of the full model. Overall, the reduced model reduces the computational expense, both in terms of computational time and memory, by several orders of magnitude for image-based PVS segments.

TABLE 3

TABLE 3. The geometrically-reduced models reduce computational cost by orders of magnitude.

Discussion

We have proposed a new mathematical and numerical framework based on topological and geometrical model reduction for computational modelling and simulation of steady and pulsatile fluid flow in deformable perivascular space networks. The reduced model is defined over a perivascular centerline network and predicts the fluid flux and average pressure in each cross-section of each network branch. By numerically comparing direct three-dimensional simulations of the fluid flow with the reduced model results for a range of physiological scenarios, we find that the reduced model accurately captures the important flow characteristics with cross-section peak pressure discrepancies ranging from 0% to 52% and peak flux discrepancies ranging from 0% to 35%. Our findings indicate that reduced model is robust with respect to physiologically relevant spatial and temporal variations in the vascular radius. Moreover and importantly, the computational cost of the reduced model is several orders of magnitude lower than that of the corresponding full model.

While geometrically-reduced network models of pulsatile blood flow have become a standard computational tool [19, 23, 44], network models of perivascular fluid flow have mainly focused either on quantifying flow resistance [7, 10] or predicting steady flow [45]. In the latter, Tithof et al present the results of a network model of glymphatic flow under different parameters, using resistance models to compute flow in idealized domains. For the open channel flow, they compute the flow therein via Darcy’s law v = −(κA/ν)∇p with permeability

$\kappa =\frac{1}{8}\left(R_2^2+R_1^2-\frac{R_2^2-R_1^2}{\ln \left(R_2/R_1\right)}\right).(23)$

This relationship holds under the assumption of Poiseuille flow in the open, annular channel (for which there is an analytic solution) and corresponds to the permeability required for this solution to satisfy Darcy’s law. For steady-state flow (∂_tv = ∂_ssv = 0) driven by a constant pressure difference, the reduced model (6) simplify to the Darcy flow equation with permeability

$\kappa =\frac{1}{\alpha }.(24)$

In the idealized Model A1 scenario, the two definitions of κ (23, 24) agree, with κ = 1.36 × 10^–4 mm², and thus the models coincide within this regime.

Rey and Sarntinoranont [13] also introduced two hydraulic models to predict fluid flow induced by blood pressure wave pulsations, and in particular net flow and transport. Their models also capture the pulsatile flow generated by the volume changes induced by a pulsating inner boundary, but under other modelling assumptions and without considering bifurcations, and thus differ from the one considered here. However, their peak fluid velocities of the order tens of μm/s is of the same order as the fluid velocities predicted in single branches here (Models A2, B2, B3), as are the pressures on the order of up to 0.3 Pa.

Several different bifurcation conditions have been proposed in the literature. In one-dimensional blood flow models, the most common conditions are conservation of flux combined with continuity of pressure [44, 46]. These conditions may be imposed directly on the pressure and flux solution variables [44], or weakly in the variational formulation [46]. Here, we also enforce conservation of flux, but in place of the strong pressure continuity condition, we weakly impose the continuity of the normal stress. This approach gives a natural setting for Stokes flow and allows for a compatible variational formulation using a Lagrange multiplier space.

In terms of limitations, we here focus on models of perivascular flow and the effect of vascular pulsations on perivascular flow, and not on the full interplay between vascular, perivascular and interstitial flow and deformation, nor on the transfer across the blood-brain barrier or the glial limitans. For healthy arterial and venous regions, in which the blood flow dynamics dominate the perivascular flow and pressure, we expect this one-way (vascular-to-perivascular) coupling to capture the leading order dynamics. Moreover, in light of the expected high resistance of the interstitial space [13, 45, 47, 48], we expect the perivascular-interstitial transfer and interstitial flow to be relatively small under physiological conditions. However, in light of the importance of quantifying and characterizing the different potential pathways, coupled fluid dynamics in vascular, perivascular and interstitial spaces will be considered in subsequent work.

We here consider open (in contrast to porous) domains. This is an appropriate modelling choice for surface perivascular spaces surrounding arteries or veins [6, 8]. For parenchymal perivascular spaces, within the pial-glial interface or within the smooth muscle cell basement membranes [49], however, a porous media representation may be more appropriate. In such a case, the Stokes flow (1) are naturally replaced by a Darcy or Brinkman flow model with an additional permeability κ [50]. The analogous reduced model [corresponding to (6)] would include an additional lower order term for the flux $\widehat{q}$ weighted by this permeability. For parenchymal and capillary perivascular spaces, we would also expect the coupled interplay between vascular, perivascular and interstitial spaces to be non-negligible.

Surface PVSs may be of different shapes ranging from annular cylinders with no or some ellipticity to fully separated segments [7], or be defined as more irregular expansions of the subarachnoid space [6, 51]. The image-based vascular geometries used here define high-fidelity representations of inner boundaries of human surface PVSs. However, the representation of surrounding PVSs as annular structures is clearly an approximation, and a response to the lack of appropriate three-dimensional data of human surface PVSs. An interesting point in this regard, and an opportunity for further study, is the quantification of the model error introduced by approximating these non-regular structures by elliptic annular cylinders with a fixed centerline. We would expect more irregular geometries to induce larger differences between the full and reduced models, but the relative importance and role of ellipticity and other geometrical irregularities remain undetermined.

Our simulations rely on a high-order discretization in space to ensure stability of the model, combined with a first-order discretization in time. A temporal sensitivity analysis on key output quantities i.e. pressure and cross-section flux showed expected convergence as the time resolution is reduced, and was used to determine the employed time step. The use of a higher-order discretization in time could also be considered. We also note that we have considered simplified (prescribed traction) boundary conditions at the PVS inlet and outlets. Compliance or resistance-based boundary conditions could of course also be considered, e.g., as in previous work [18], or [37]. We have focused on cardiac pulse wave-induced wall motion and vasomotion, two physiological factors that generate changes in vascular radius of up to 15% [31, 38] and only moderate wall velocities. However, the vascular and perivascular diameters may change more dramatically. For instance, Enger et al [52] report of a nearly 40% increase and 50% decrease in arteriole diameter during cortical spreading depression, and intriguingly the vascular and perivascular wall motions may differ between e.g., sleep states [53]. If these changes lead to significantly higher wall velocities than those considered here, we would expect a further breakdown of the reduced model assumptions, specifically assumption V, which in turn would be expected to impact the accuracy of the reduced models.

While many aspects of brain influx and clearance remain enigmatic, perivascular fluid flow along the cerebral vasculature is widely recognized as a key transport mechanism. The computationally inexpensive yet accurate reduced models presented here give an efficient and flexible framework for computational modelling and simulation of pulsatile flow in idealized or realistic networks including complete representations of e.g., the cerebral arteries or veins and many generations of arterioles/capillaries. This framework thus establishes a foundation for future computational studies of perivascular flow to improve our understanding of brain transport.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.5281/zenodo.5729484.

Author Contributions

CD-C, IG, and MR conceived the experiments. CD-C created the PVS geometries. CD-C, IG, and MR developed the simulation code. CD-C and IG conducted the experiments and analyzed the results. CD-C, IG, and MR created figures. CD-C, IG, and MR wrote the paper. All authors reviewed and approved the final manuscript.

Funding

This study has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement 714892.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Miroslav Kuchta (Simula Research Laboratory) for constructive discussion on topics related to the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.882260/full#supplementary-material

References

1. Rennels ML, Gregory TF, Blaumanis OR, Fujimoto K, Grady PA. 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 (1985) 326:47–63. doi:10.1016/0006-8993(85)91383-6

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Carare RO, Bernardes-Silva M, Newman TA, Page AM, Nicoll JAR, Perry VH, et al. Solutes, but Not Cells, drain from the Brain Parenchyma along Basement Membranes of Capillaries and Arteries: Significance for Cerebral Amyloid Angiopathy and Neuroimmunology. Neuropathol Appl Neurobiol (2008) 34:131–44. doi:10.1111/j.1365-2990.2007.00926.x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Iliff JJ, Wang M, Liao Y, Plogg BA, Peng W, Gundersen GA, et al. A Paravascular Pathway Facilitates CSF Flow through the Brain Parenchyma and the Clearance of Interstitial Solutes, Including Amyloid β. Sci Transl Med (2012) 4:147ra111. doi:10.1126/scitranslmed.3003748

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Wardlaw JM, Benveniste H, Nedergaard M, Zlokovic BV, Mestre H, Lee H, et al. Perivascular Spaces in the Brain: Anatomy, Physiology and Pathology. Nat Rev Neurol (2020) 16:137–53. doi:10.1038/s41582-020-0312-z

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Zhang ET, Inman CB, Weller RO. Interrelationships of the Pia Mater and the Perivascular (Virchow-Robin) Spaces in the Human Cerebrum. J Anat (1990) 170:111–23.

PubMed Abstract | Google Scholar

6. Bedussi B, Almasian M, de Vos J, VanBavel E, Bakker EN. Paravascular Spaces at the Brain Surface: Low Resistance Pathways for Cerebrospinal Fluid Flow. J Cereb Blood Flow Metab (2018) 38:719–26. doi:10.1177/0271678x17737984

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Tithof J, Kelley DH, Mestre H, Nedergaard M, Thomas JH. Hydraulic Resistance of Periarterial Spaces in the Brain. Fluids Barriers CNS (2019) 16:19–3. doi:10.1186/s12987-019-0140-y

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Min Rivas F, Liu J, Martell BC, Du T, Mestre H, Nedergaard M, et al. Surface Periarterial Spaces of the Mouse Brain Are Open, Not Porous. J R Soc Interf (2020) 17:20200593. doi:10.1098/rsif.2020.0593

CrossRef Full Text | Google Scholar

9. Martinac AD, Bilston LE. Computational Modelling of Fluid and Solute Transport in the Brain. Biomech Model Mechanobiol (2019) 19:1–20. doi:10.1007/s10237-019-01253-y

CrossRef Full Text | Google Scholar

10. Faghih MM, Sharp MK. Is Bulk Flow Plausible in Perivascular, Paravascular and Paravenous Channels? Fluids Barriers CNS (2018) 15:17. doi:10.1186/s12987-018-0103-8

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Asgari M, de Zélicourt D, Kurtcuoglu V. Glymphatic Solute Transport Does Not Require Bulk Flow. Sci Rep (2016) 6:38635–11. doi:10.1038/srep38635

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Diem AK, MacGregor Sharp M, Gatherer M, Bressloff NW, Carare RO, Richardson G. Arterial Pulsations Cannot Drive Intramural Periarterial Drainage: Significance for Aβ Drainage. Front Neurosci (2017) 11:475. doi:10.3389/fnins.2017.00475

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Rey J, Sarntinoranont M. Pulsatile Flow Drivers in Brain Parenchyma and Perivascular Spaces: a Resistance Network Model Study. Fluids Barriers CNS (2018) 15:20. doi:10.1186/s12987-018-0105-6

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Keith Sharp M, Carare RO, Martin BA. Dispersion in Porous media in Oscillatory Flow between Flat Plates: Applications to Intrathecal, Periarterial and Paraarterial Solute Transport in the central Nervous System. Fluids Barriers CNS (2019) 16:13. doi:10.1186/s12987-019-0132-y

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lloyd RA, Stoodley MA, Fletcher DF, Bilston LE. The Effects of Variation in the Arterial Pulse Waveform on Perivascular Flow. J Biomech (2019) 90:65–70. doi:10.1016/j.jbiomech.2019.04.030

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kedarasetti RT, Drew PJ, Costanzo F. Arterial Pulsations Drive Oscillatory Flow of CSF but Not Directional Pumping. Sci Rep (2020) 10:10102–12. doi:10.1038/s41598-020-66887-w

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kedarasetti RT, Turner KL, Echagarruga C, Gluckman BJ, Drew PJ, Costanzo F. Functional Hyperemia Drives Fluid Exchange in the Paravascular Space. Fluids Barriers CNS (2020) 17:52–25. doi:10.1186/s12987-020-00214-3

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Daversin-Catty C, Vinje V, Mardal K-A, Rognes ME. The Mechanisms behind Perivascular Fluid Flow. PLOS ONE (2020) 15:e0244442. doi:10.1371/journal.pone.0244442

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Olufsen MS. Structured Tree Outflow Condition for Blood Flow in Larger Systemic Arteries. Am J Physiol Heart Circulatory Physiol (1999) 276:H257–H268. doi:10.1152/ajpheart.1999.276.1.h257

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Sherwin SJ, Franke V, Peiró J, Parker K. One-dimensional Modelling of a Vascular Network in Space-Time Variables. J Eng Math (2003) 47:217–50. doi:10.1023/b:engi.0000007979.32871.e2

CrossRef Full Text | Google Scholar

21. D’Angelo C, Quarteroni A. On the Coupling of 1d and 3d Diffusion-Reaction Equations: Application to Tissue Perfusion Problems. Math Models Methods Appl Sci (2008) 18:1481–504. doi:10.1142/S0218202508003108

CrossRef Full Text | Google Scholar

22. Lesinigo M, D’Angelo C, Quarteroni A. A Multiscale Darcy-Brinkman Model for Fluid Flow in Fractured Porous media. Numer Math (2011) 117:717–52. doi:10.1007/s00211-010-0343-2

CrossRef Full Text | Google Scholar

23. Coccarelli A, Carson JM, Aggarwal A, Pant S. A Framework for Incorporating 3d Hyperelastic Vascular wall Models in 1d Blood Flow Simulations. Biomech Model Mechanobiol (2021) 20:1231. doi:10.1007/s10237-021-01437-5

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Köppl T, Vidotto E, Wohlmuth B. A 3D‐1D Coupled Blood Flow and Oxygen Transport Model to Generate Microvascular Networks. Int J Numer Meth Biomed Engng (2020) 36:e3386. doi:10.1002/cnm.3386

CrossRef Full Text | Google Scholar

25. Koch T, Schneider M, Helmig R, Jenny P. Modeling Tissue Perfusion in Terms of 1d-3d Embedded Mixed-Dimension Coupled Problems with Distributed Sources. J Comput Phys (2020) 410:109370. doi:10.1016/j.jcp.2020.109370

CrossRef Full Text | Google Scholar

26. Vidotto E, Koch T, Köppl T, Helmig R, Wohlmuth B. Hybrid Models for Simulating Blood Flow in Microvascular Networks. Multiscale Model Simul (2019) 17:1076–102. doi:10.1137/18m1228712

CrossRef Full Text | Google Scholar

27. Cattaneo L, Zunino P. A Computational Model of Drug Delivery through Microcirculation to Compare Different Tumor Treatments. Int J Numer Meth Biomed Engng (2014) 30:1347–71. doi:10.1002/cnm.2661

CrossRef Full Text | Google Scholar

28. Possenti L, di Gregorio S, Gerosa FM, Raimondi G, Casagrande G, Costantino ML, et al. A Computational Model for Microcirculation Including Fahraeus-Lindqvist Effect, Plasma Skimming and Fluid Exchange with the Tissue Interstitium. Int J Numer Meth Biomed Engng (2018) 35:e3165. doi:10.1002/cnm.3165

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Possenti L, Cicchetti A, Rosati R, Cerroni D, Costantino ML, Rancati T, et al. A Mesoscale Computational Model for Microvascular Oxygen Transfer. Ann Biomed Eng (2021) 49:3356. doi:10.1007/s10439-021-02807-x

PubMed Abstract | CrossRef Full Text | Google Scholar

30.Aneurisk-Team. AneuriskWeb Project Website. [Dataset] (2012). Available at: http://ecm2.mathcs.emory.edu/aneuriskweb Accessed April 19, 2021.

Google Scholar

31. Mestre H, Tithof J, Du T, Song W, Peng W, Sweeney AM, et al. Flow of Cerebrospinal Fluid Is Driven by Arterial Pulsations and Is Reduced in Hypertension. Nat Commun (2018) 9:4878. doi:10.1038/s41467-018-07318-3

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Daversin-Catty C. PVS Meshing Tools. [Dataset]. Github (2020). Available at: https://github.com/cdaversin/PVS-meshing-tools

Google Scholar

33. Antiga L, Piccinelli M, Botti L, Ene-Iordache B, Remuzzi A, Steinman DA. An Image-Based Modeling Framework for Patient-specific Computational Hemodynamics. Med Biol Eng Comput (2008) 46:1097–112. doi:10.1007/s11517-008-0420-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Schlömer N. Meshio v4.3.10. [Dataset]. Zenodo (2020). doi:10.5281/zenodo.4555995

CrossRef Full Text | Google Scholar

35. Geuzaine C, Remacle J-F. Gmsh: A 3-d Finite Element Mesh Generator with Built-In Pre- and post-processing Facilities. Int J Numer Meth Engng (2009) 79:1309–31. doi:10.1002/nme.2579

CrossRef Full Text | Google Scholar

36. San Martín J, Smaranda L, Takahashi T. Convergence of a Finite element/ALE Method for the Stokes Equations in a Domain Depending on Time. J Comput Appl Math (2009) 230:521–45. doi:10.1016/j.cam.2008.12.021

CrossRef Full Text | Google Scholar

37. Ladrón-de-Guevara A, Shang JK, Nedergaard M, Kelley DH. Perivascular Pumping in the Mouse Brain: Improved Boundary Conditions Reconcile Theory, Simulation, and experiment. J Theor Biol (2022) 542:111103. doi:10.1016/j.jtbi.2022.111103

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Aldea R, Weller RO, Wilcock DM, Carare RO, Richardson G. Cerebrovascular Smooth Muscle Cells as the Drivers of Intramural Periarterial Drainage of the Brain. Front Aging Neurosci (2019) 11:1. doi:10.3389/fnagi.2019.00001

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Alnæs MS, Blechta J, Hake J, Johansson A, Kehlet B, Logg A, et al. The FEniCS Project Version 1.5. Archive Numer Softw (2015) 3:9–23. doi:10.11588/ans.2015.100.20553

CrossRef Full Text | Google Scholar

40. Daversin-Catty C, Vinje V, Mardal K-A, Rognes ME Mechanisms-behind-pvs-flow-v1.0. [Dataset] (2020). doi:10.5281/zenodo.3890133

CrossRef Full Text | Google Scholar

41. Daversin-Catty C, Richardson CN, Ellingsrud AJ, Rognes ME Abstractions and Automated Algorithms for Mixed-Dimensional Finite Element Methods. ACM Trans Math Softw (2021) 47:1–36. doi:10.1145/3471138

CrossRef Full Text | Google Scholar

42. Daversin-Catty C, Gjerde IG, Rognes ME. Geometrically-reduced-pvs-flow-v1.0. [Dataset] (2021). doi:10.5281/zenodo.5729484

CrossRef Full Text | Google Scholar

43. Kuchta M. Assembly of Multiscale Linear PDE Operators. In: Lecture Notes in Computational Science and Engineering. Cham: Springer International Publishing (2020). p. 641–50. doi:10.1007/978-3-030-55874-1_63

CrossRef Full Text | Google Scholar

44. Olufsen M, Nadim A. On Deriving Lumped Models for Blood Flow and Pressure in the Systemic Arteries. Mbe (2004) 1:61–80. doi:10.3934/mbe.2004.1.61

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Tithof J, Boster KAS, Bork PAR, Nedergaard M, Thomas JH, Kelley DH. A Network Model of Glymphatic Flow under Different Experimentally-Motivated Parametric Scenarios. bioRxiv:104258 (2021). doi:10.1101/2021.09.23.461519

CrossRef Full Text | Google Scholar

46. Notaro D, Cattaneo L, Formaggia L, Scotti A, Zunino P. A Mixed Finite Element Method for Modeling the Fluid Exchange between Microcirculation and Tissue Interstitium. In: G Ventura, and E Benvenuti, editors. Advances in Discretization Methods: Discontinuities, Virtual Elements, Fictitious Domain Methods. Cham: Springer International Publishing (2016). p. 3–25. doi:10.1007/978-3-319-41246-7_1

CrossRef Full Text | Google Scholar

47. Holter KE, Kuchta M, Mardal K-A. Sub-voxel Perfusion Modeling in Terms of Coupled 3d-1d Problem. In: FA Radu, K Kumar, I Berre, JM Nordbotten, and IS Pop, editors. Numerical Mathematics and Advanced Applications ENUMATH 2017. Cham: Springer International Publishing (2019). p. 35–47. doi:10.1007/978-3-319-96415-7_2

CrossRef Full Text | Google Scholar

48. Vinje V, Eklund A, Mardal KA, Rognes ME, Støverud KH. Intracranial Pressure Elevation Alters CSF Clearance Pathways. Fluids Barriers CNS (2020) 17:29–19. doi:10.1186/s12987-020-00189-1

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Albargothy NJ, Johnston DA, MacGregor-Sharp M, Weller RO, Verma A, Hawkes CA, et al. Convective Influx/glymphatic System: Tracers Injected into the CSF Enter and Leave the Brain along Separate Periarterial Basement Membrane Pathways. Acta Neuropathol (2018) 136:139–52. doi:10.1007/s00401-018-1862-7

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Brinkman HC. A Calculation of the Viscous Force Exerted by a Flowing Fluid on a Dense Swarm of Particles. Appl Sci Res (1949) 1:27–34. doi:10.1007/bf02120313

CrossRef Full Text | Google Scholar

51. Vinje V, Bakker ENTP, Rognes ME. Brain Solute Transport Is More Rapid in Periarterial Than Perivenous Spaces. Scientific Rep (2021) 11:16085. doi:10.1038/s41598-021-95306-x

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Enger R, Tang W, Vindedal GF, Jensen V, Johannes Helm P, Sprengel R, et al. Dynamics of Ionic Shifts in Cortical Spreading Depression. Cereb Cortex (2015) 25:4469–76. doi:10.1093/cercor/bhv054

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Bojarskaite L, Bjørnstad DM, Pettersen KH, Cunen C, Hermansen GH, Åbjørsbråten KS, et al. Astrocytic Ca2+ Signaling Is Reduced during Sleep and Is Involved in the Regulation of Slow Wave Sleep. Nat Commun (2020) 11:3240–16. doi:10.1038/s41467-020-17062-2

PubMed Abstract | CrossRef Full Text | Google Scholar