{'fraction_non_alphanumeric': 0.03580810965162764, 'fraction_numerical': 0.015276984580239863, 'mean_word_length': 4.806831371248549, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.


Introduction
Flow of cerebrospinal fluid (CSF) in perivascular spaces (PVSs) is a key transport mechanism in and around the brain [1][2][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][5][6][7][8] , PVSs are typically represented as (elliptic) annular structures or pipes 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][12][13][14][15][16][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][20][21][22][23][24][25][26][27][28][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.  Table 1. Geometrical or numerical PVS domain characteristics for domains A, B, C. denotes an approximate domain length, = 2 1 is the range of the arterial diameters, pvs indicates the range of widths of the perivascular space ( pvs = 2 − 1 , so that 2 = 2.95 1 ) cells and vertices indicate the number of mesh cells and mesh vertices respectively for the full (2D or 3D) model and reduced models, and ℎ max denotes the maximal mesh cell size for each mesh. The vertices for the one-dimensional geometries are uniformly spaced in the interior of the domain.
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.

PVS geometries (3D and 1D)
In general, we consider a perivascular tree-like domain Ω consisting of a network of branching generalized annular cylinders Ω , with Ω ⊆ ∪ ∈ Ω , spatial coordinates ∈ Ω and time ≥ 0. The boundary is denoted Ω, with boundary normal . We assume that each generalized annular cylinder Ω has a well-defined and oriented, topologically one-dimensional centerline Λ with coordinate . We set Λ = ∪ ∈ Λ . Along , we define the cross-sections = ( , ) of Λ with area = ( , ). We denote the inner radius of Ω by 1 and the outer radius of Ω by 2 ; these radii will in practice vary with , and the angular coordinate . We denote the set of bifurcation points i.e. the points at which the centerlines of branches meet by 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). Three-dimensional PVS flow in geometries A and C have been studied previously 18 and will be used for comparison. 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,30 . 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 .
The 3D PVS construction and the 1D centerline extraction are performed using PVS-meshing-tools 31 , largely based on VMTK 32 . 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 31 using meshio 33 and GMSH 34 . The centerline meshes consist of topologically one-dimensional intervals embedded in three dimensions. The bifurcation points ∈ B ⊂ Ω 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/19
A 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 a realistic 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 a realistic perivascular space including a bifurcation.

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 2 − 10 4 ), a mean flow speed of up to 60 m/s, and parabolic flow profiles 30 . 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 Ω = Ω( ) representing the PVS. The fluid velocity = ( , ) for ∈ Ω( ) at time and the CSF pressure = ( , ) then solve the following system of time-dependent partial differential equations (PDEs) 18,35 : where is the fluid density and is the dynamic fluid viscosity. To model CSF at body temperature, we set the fluid density to = 10 3 kg/m 3 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 Ω( ) at time > 0 is given by a deformation of the reference domain: We denote the domain velocity associated with by (thus = ).

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˜= 0, or (ii) a constant-in-time pressure gradient Δ˜> 0 by setting in = out Δ˜at the inlet, letting˜o ut = 0 at the outlet furthest from the inlet with distance out , and setting˜o ut at any other outlets such that the average pressure gradient over each branch path (˜−˜o ut )/ out is constant and equal to the prescribed pressure gradient Δ˜mmHg/m. This static pressure difference can represent e.g. a hydrostatic pressure difference, a venous pressure differential, or some other systemic pressure difference.
On the inner and outer PVS walls (along the length of the PVS), we set the fluid velocity to match a known, prescribed domain velocity = ( , ). For the inner PVS wall, we either (i) consider a rigid wall and set = = 0, or (ii) impose a pulsating wall displacement: with reference to the initial (fixed) mesh with coordinates and prescribe = = . To represent wall motion induced by the cardiac pulse wave, we let the amplitude be defined by the juxtaposition of an experimentally-observed wall motion time series 30 either applied uniformly along the length of the PVS or as a travelling wave along the PVS length with wave speed = 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 36 , and an amplitude of 7.5% of the initial inner radius 1 . We note that for all models, the wall moves in the normal (radial) direction only. For the outer PVS wall Ω outer , we set = = 0.
The system starts at rest with = = 0 at = 0. The system reaches the periodic steady state nearly immediately, and we report results starting from the first cycle.

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 37 . For each branch Ω ( ) with centerline Λ and local coordinate system ( , , ), where represents the path length (or axial coordinate), 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;

4/19
(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. = ( , ); (V) Axial velocity profile The axial velocity , i.e. the velocity component in the axial direction can be decomposed in the form where vp is a given velocity profile varying radially only,ˆis to be determined.
For the velocity profile vp , we here choose a normalized annular Poiseuille flow: This velocity profile is parabolic in (as for Poiseuille flow in a cylinder) with a logarithmic correction that accounts for the annulus.
In particular, the domain velocity 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.

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 37 : find the cross-section fluxˆ=ˆ( , ) and the cross-section average pressureˆ=ˆ( , ) such that for each centerline Λ (denotingˆ| Λ =ˆandˆ| Λ =ˆ): hold.ˆ( Moreover, = ( , ) denotes the cross-section area, whileˆ=ˆ( , ) is a lumped flow parameter that depends on the domain geometry and the choice of velocity profile vp : and where¯v p is the velocity profile integrated over each cross-section: We also define the (one-dimensional) normal stress induced byˆandˆ: which corresponds to an average of the axial ( -)component of the normal stress in (2) over each cross-section.

5/19
At the bifurcation points ∈ B ⊂ Ω, we impose the following two conditions representing conservation of flux and continuity of normal stress, respectively: where Λ and Λ 1 , Λ 2 represent the centerlines of the parent and two daughter branches, respectively, associated with the bifurcation point and · = · ( ) where · denotes the map from three-dimensional bifurcation point to the one-dimensional centerline coordinate for each branch Ω · . The system (6) defines a set of equations for each branch centerline Λ 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ˆcf. (10). In this manner, the (one-dimensional) solutionsˆandˆof the reduced model (6) define approximations of the (three-dimensional) axial flux and pressure solving (1) integrated or averaged over each cross-section:ˆ( The factor originates from integrating in cylindrical coordinates. We note that the wall velocity , which defines a boundary condition for the full PVS model (1), enters as a body force in the reduced model (6).

Numerical solution and software
We solve the full PVS equations (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 38 , and is openly available 39 .
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 37 . The finite element mesh T of the centerline Λ is composed of mesh segments T , one for each centerline branch Λ . Each mesh segment is a mesh consisting of intervals embedded in R 3 . We label the set of bifurcation points B, inlet points I and outlet points O, and define the following finite element spaces: • The flux space ℎ is the space of continuous piecewise quadratics over T for each .
• The (average) pressure space ℎ is the space of continuous piecewise linears on T .
• The Lagrange multiplier space ℎ = R where 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.

6/19
where (or ) is simply the entry of the vector (or ) corresponding to bifurcation point , and we define the natural jump: The right-hand side linear form is: where the superscript ( ) 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 38 . The solver, and in particular the definition of the partially continuous flux space, builds on mixed-domain features 40 and relies on the latest development version of FEniCS.

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.
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 ( = 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 = 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.
For each cross section , we compare a numerical approximation¯ℎ of the average of the full pressure solution: for a quadrature scheme with points and weights defined over and an approximation | | of the cross-section area 37 . 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 41 .
Similarly, we compute a numerical flux approximation¯ℎ by the same numerical integration of the axial velocity over the cross-section . We define the total relative model discrepancy in the flux by and similarly for the pressure ( ). Finally, we define the net flux associated with the velocity = ( , ) as: and the corresponding quantity associated with the fluxˆ=ˆ( ): where the integration in time is over one period [0, ].

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 equations (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.

Reduced model exactly predicts pressure-driven axisymmetric flow characteristics
Flow in an axisymmetric annular cylinder of length ℓ driven by a constant pressure difference Δ (Model A1) is described by the analytic expression:ˆ( where is the lumped flow parameter given by (8) ( = 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.

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, Figure 2b Pa for the reduced model. The relative model discrepancy in the peak cross-section flux (difference between the full and reduced peak flux) at the inlet is 4.1% and in the peak cross-section (average) pressure is 5.3%. There is thus a small discrepancy between the two models, as expected by the violation of the reduced model assumptions.

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 ) 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 ( = 0.5), the reduced pressure approximationˆvaries nearly linearly along the length of the domain as expected, and the reduced flux approximationˆis essentially constant along the centerline with valueˆ= 4.28 × 10 −4 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 total relative model discrepancy (18)

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 (Figure 4c, Figure 5), but with angularly varying velocity profiles (Figure 4b, Figure 5).
For the cardiac wall motion, the cross-section average of the full pressure¯ℎ ranges from −0.05 to 0.26 Pa, while the full cross-section flux¯ℎ 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 (Figure 4d-Figure 4g). For the reduced model, the cross-section pressureˆℎ ranges from −0.06 to 0.29 Pa, while the cross-section fluxˆℎ 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, we observe that the reduced model slightly overestimates the peak pressure and flux when compared to the full model ( Figure 4e, Figure 4f). The relative difference in peak positive pressure between the two models at this cross-section is 19% and 30% in peak negative pressure. For the flux, the corresponding model discrepancies are 1.2% and 11%.
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.001m/s 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). 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ˆℎ matches pulsatile behaviour of the average cross-section pressure in the full model¯ℎ (Figure 5b) but that the peak amplitude is higher. The largest model differences in pressure at lower cross-section is at the peak pressure; there the relative difference in peak pressure is 52%. The similar observations hold for the flux, but the model discrepancies are lower: the relative difference in peak flux is 15%. 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).

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 S2). 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-Womersleytype velocity profiles. Comparing the full and reduced average pressure and flux at the time of peak velocity (Figure 6f, Figure 6j), we note that the reduced model captures the qualitative and quantitative flow and pressure characteristics. The bifurcation conditions are satisfied at the bifurcation point (Figure 6c, Figure 6g 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 peak flux model discrepancies at the inlet and outlets, we note that the discrepancy is largest at larger daughter outlet with a relative difference of 12%. Comparing the full and reduced peak pressures at the upper, middle and lower cross-sections, we find relative model differences of 1%, 4%,   Table 3. The geometrically-reduced models reduce computational cost by orders of magnitude. Number of degrees of freedom . . . , computational time (average time for a single time step) and memory usage (peak memory usage throughout the simulation) for the full models (2D/3D) and reduced models (1D). and 16%. The analogous numbers for the fluxes are 1%, 34%, and 5%. Thus, the model discrepancies for the flux are larger near the bifurcation region (Figure 6f, Figure 6j). 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 L, and 2.9 × 10 −5 L for the reduced model, corresponding to a relative difference of 17%.

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)

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,42 , network models of perivascular fluid flow have mainly focused either on quantifying flow resistance 7,10 or predicting steady flow 43 In the idealized Model A1 scenario, the two definitions of ( (21) and (22)) agree, with = 1.36 × 10 −4 mm 2 , 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 42,44 . These conditions may be imposed directly on the pressure and flux solution variables 42 , or weakly in the variational formulation 44 . 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,43,45,46 , 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 47 , however, a porous media representation may be more appropriate. In such a case, the Stokes flow equations (1) are naturally replaced by a Darcy or Brinkman flow model with an additional permeability 48 . The analogous reduced model (corresponding to (6)) would include an additional lower order term for the fluxˆ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.
Furthermore, we have approximated the PVS as an (elliptic) annular structure, while surface PVSs may be of different shapes 6,7,49 . An interesting point is the quantification of the model error introduced by approximating these non-regular structures by elliptic annular cylinders with a fixed centerline. Gjerde et al 37 addresses this point numerically and via theoretical analysis, including the balance between numerical and model errors. Finally, 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 16/19 work 18 . 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% 30,36 and only moderate wall velocities. However, the vascular and perivascular diameters may change more dramatically. For instance, Enger et al 50 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 51 . 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 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.