ORIGINAL RESEARCH article
Sec. Statistical and Computational Physics
Geometrically Reduced Modelling of Pulsatile Flow in Perivascular Networks
- 1Simula Research Laboratory, Oslo, Norway
- 2Department of Mathematics, University of Bergen, Bergen, Norway
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.
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 . 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 . Theoretical models have quantified the resistance in PVS networks , while detailed numerical simulations can predict perivascular fluid velocities and pressures in idealized [11–17] and image-based geometries . 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 , 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 Ωi, with Ω ⊆∪i∈IΩi, spatial coordinates x ∈ Ω and time t ≥ 0. The boundary is denoted ∂Ω, with boundary normal n. We assume that each generalized annular cylinder Ωi has a well-defined and oriented, topologically one-dimensional centerline Λi with coordinate s. We set Λ = ∪i∈IΛi. Along s, we define the cross-sections Ci = Ci(s, t) of Λi with area Ai = Ai (s, t). We denote the inner radius of Ωi by
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 , 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  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. 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.
The 3D PVS construction and the 1D centerline extraction are performed using PVS-meshing-tools , largely based on VMTK . 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  using meshio  and GMSH . The centerline meshes consist of topologically one-dimensional intervals embedded in three dimensions. The bifurcation points
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 (102–104) for 1 μm spherical particles transported at low Reynolds number , a mean flow speed of up to 60 μm/s, and parabolic flow profiles . 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]:
where ρ is the fluid density and μ is the dynamic fluid viscosity. To model CSF at body temperature, we set the fluid density to ρ = 103 kg/m3 and the dynamic viscosity to μ = 0.697 × 10–3 Pa s. As in our previous full models of perivascular flow , 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
2.3 Boundary Conditions, Initial Conditions and Periodicity
At the PVS ends, we prescribe a traction condition corresponding to a known, applied pressure
We either prescribe (i) zero pressure at both ends
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:
with reference to the initial (fixed) mesh with coordinates X and prescribe
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 Ωi(t) with centerline Λi 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 vs, i.e., the velocity component in the axial direction can be decomposed in the form
where vvp is a given velocity profile varying radially only,
For the velocity profile vvp, we here choose a normalized annular Poiseuille flow:
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
Moreover, Ai = Ai (s, t) denotes the cross-section area, while
We also define the (one-dimensional) normal stress induced by
which corresponds to an average of the axial (s-)component of the normal stress in (2) over each cross-section
At the bifurcation points
where Λp and
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
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 . This solver builds on the standard FEniCS finite element software suite , and is openly available .
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
• The flux space Vh is the space of continuous piecewise quadratics over
• The (average) pressure space Qh is the space of continuous piecewise linears on
• The Lagrange multiplier space
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 tk, given
for all finite element test functions ψ ∈ Vh, ϕ ∈ Qh, and ξ ∈ Rh. The left-hand side bilinear form a is defined by:
where λb (or ξb) is simply the entry of the vector λ (or ξ) corresponding to bifurcation point b, and we define the natural jump:
The right-hand side linear form L is:
where the superscript iI (iO) 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 . The solver, and in particular the definition of the partially continuous flux space, builds on mixed-domain features  and relies on the latest development version of FEniCS. All data and source code are available via Zenodo .
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. Overview of computational models parameterized by domain, prescribed pressure gradient
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
for a quadrature scheme with points xk and weights wk defined over C and an approximation |C| of the cross-section area. Here nC 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 .
With this in hand we define the percentagewise relative model discrepancy Eq(t) in the flux by
and similarly for the pressure Ep. 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 Eq(T) and Ep(T), where T denotes the final time.
For pulsatile flow, we typically compare the percentagewise relative error in peak pressure
and the corresponding quantity associated with the flux
where the integration in time is over one period [0, T].
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:
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/m2, 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
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
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 R1 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 ) 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
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
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. 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. 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
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
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  (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
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−3 μL/s, and 1.2 × 10−3 μL/s and 7.9, ×, 10−4 μ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 (sout) with a relative difference
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−5, and 2.9 × 10−5 μ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.
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 . 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
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
Rey and Sarntinoranont  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 , or weakly in the variational formulation . 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 , 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 κ . The analogous reduced model [corresponding to (6)] would include an additional lower order term for the flux
Surface PVSs may be of different shapes ranging from annular cylinders with no or some ellipticity to fully separated segments , 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 , or . 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  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 . 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.
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.
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.
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.
We thank Miroslav Kuchta (Simula Research Laboratory) for constructive discussion on topics related to the manuscript.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.882260/full#supplementary-material
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
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
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
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
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
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
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
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
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
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
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
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
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
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
30.Aneurisk-Team. AneuriskWeb Project Website. [Dataset] (2012). Available at: http://ecm2.mathcs.emory.edu/aneuriskweb Accessed April 19, 2021.
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
32. Daversin-Catty C. PVS Meshing Tools. [Dataset]. Github (2020). Available at: https://github.com/cdaversin/PVS-meshing-tools
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
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
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
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
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
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
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
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
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
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
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
Keywords: biomedical flows, low-dimensional models, bifurcation, variational methods, computational methods
Citation: Daversin-Catty C, Gjerde IG and Rognes ME (2022) Geometrically Reduced Modelling of Pulsatile Flow in Perivascular Networks. Front. Phys. 10:882260. doi: 10.3389/fphy.2022.882260
Received: 23 February 2022; Accepted: 06 April 2022;
Published: 11 May 2022.
Edited by:André H. Erhardt, Weierstrass Institute for Applied Analysis and Stochastics (LG), Germany
Reviewed by:Kartik Jain, University of Twente, Netherlands
Jeff Heys, Montana State University, United States
Copyright © 2022 Daversin-Catty , Gjerde and Rognes . This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Cécile Daversin-Catty , firstname.lastname@example.org