The mechanisms behind perivascular fluid flow

Flow of cerebrospinal fluid (CSF) in perivascular spaces (PVS) is one of the key concepts involved in theories concerning clearance from the brain. Experimental studies have demonstrated both net and oscillatory movement of microspheres in PVS (Mestre et al. (2018), Bedussi et al. (2018)). The oscillatory particle movement has a clear cardiac component, while the mechanisms involved in net movement remain disputed. Using computational fluid dynamics, we computed the CSF velocity and pressure in a PVS surrounding a cerebral artery subject to different forces, representing arterial wall expansion, systemic CSF pressure changes and rigid motions of the artery. The arterial wall expansion generated velocity amplitudes of 60–260 μm/s, which is in the upper range of previously observed values. In the absence of a static pressure gradient, predicted net flow velocities were small (<0.5 μm/s), though reaching up to 7 μm/s for non-physiological PVS lengths. In realistic geometries, a static systemic pressure increase of physiologically plausible magnitude was sufficient to induce net flow velocities of 20–30 μm/s. Moreover, rigid motions of the artery added to the complexity of flow patterns in the PVS. Our study demonstrates that the combination of arterial wall expansion, rigid motions and a static CSF pressure gradient generates net and oscillatory PVS flow, quantitatively comparable with experimental findings. The static CSF pressure gradient required for net flow is small, suggesting that its origin is yet to be determined.


Introduction
The glymphatic theory [1] suggests that the interaction of cerebrospinal fluid (CSF) and interstitial fluid facilitates the brain's clearance of metabolites via perivascular spaces (PVS) in a process faster than diffusion alone. Many experimental findings [2][3][4][5][6][7][8] demonstrate and support that transport is faster than diffusion, while others do not [9]. The glymphatic concept involves an influx of CSF in periarterial spaces, convective flow through the interstitium and finally efflux in perivenous spaces. Convective flow through the interstitium has been challenged [10,11], although even small convective flows may be important for large molecules [12] such as Amyloid-beta and tau. The venous efflux is also not without controversy. Tracers in the SAS have been reported to reach the ventricles without a presence in perivenous spaces [3]. The presence and direction of flow in PVS around arteries is also debated. According to the IPAD hypothesis, fluid is drained out from the brain along the basement membranes of capillaries and arterioles [13][14][15]. In particular, accumulation of amyloid in the walls of cerebral arteries has been seen in cerebral amyloid angiopathy [16]. In conclusion, cerebral fluid flow and transport is still controversial, and the many aspects of the glymphatic hypothesis are still debated almost a decade after its inception [17][18][19]. Initially, the term paravascular spaces was used for a Virchow-Robin type space [1], distinct from the perivascular intramural spaces [15]. The distinction is important in light of the controversy mentioned above. In this paper, we adopt the term perivascular space and consider flow along arteries on the pial surface (surface arteries). We assume a separate compartment within the subarachnoid space (SAS) enclosing the pial arteries, and we consider mechanisms behind flow within this compartment. On the pial surface, some studies have suggested that the perivascular space and SAS form a continuous compartment [3], while others indicate that these define separate spaces [2,20]. Furthermore, several computational studies have questioned whether arterial wall pulsations is a sufficiently effective mechanism for transport in the perivascular spaces [14,[21][22][23][24], while others [2,25] support this concept. Hence, the precise mechanisms and forces involved in perivascular flow have not yet been adequately described.
Perivascular flow appears to originate from forces associated with the cardiac cycle as travelling particles have a distinct cardiac frequency in their motion [2]. The cardiac CSF pulsation is well characterized both in terms of CSF flow and intracranial pressures (ICP) [26]. In humans, ICP is normally 7-15 mmHg [27,28], and pulsates with a temporal peak-to-peak amplitude of around 1-4 mmHg [29,30]. The pulsation is almost synchronous within the whole cranium [31], yet there is a small spatial gradient of 1-5 mmHg/m [29,32]. Less data are available on the values of ICP and in particular pressure gradients in mice. Normal mouse ICP has been reported at 4 mmHg [33] with an approximate peak-to-peak temporal amplitude of 0.5-1 mmHg [34].
Forces inducing PVS flow may originate from local arterial expansions [2], but also from systemic ICP increase and blood pressure oscillations in proximal parts of the vasculature. The forces originating from systemic components are transmitted almost instantaneously to the PVS in terms of a pressure pulsation through the incompressible CSF. Peristalsis driven by the local arterial wall pulsation has received much attention, but computational modeling and theoretical calculations (in idealized geometries) point in different directions as to whether this mechanism is sufficient for net flow [21,22,24,25,[35][36][37].
In this study, we therefore address several forces with the potential to explain both net and oscillatory fluid movement in a realistic PVS geometry. In addition to the pulsatile local arterial expansion, we evaluated systemic CSF gradients of both static and pulsatile nature as well as rigid motions of the artery. We find that all forces combined may induce PVS flow comparable to experimental observations [2], but that the magnitude of the static pressure gradient required for net flow suggests that its origin is still unclear. A small net flow velocity close to the experimental data without the presence of a pressure gradient is only achieved when the PVS geometry is long (close to the wavelength of the arterial pulse).

Methods
To predict flow characteristics and detailed flow patterns in PVS surrounding pial surface arteries, we created several computational models of a CSF-filled PVS surrounding a bifurcating cerebral artery segment (Fig 1A and 1B). This surface PVS was represented as an open (unobstructed) space, deforming in time, and fluid flow in the PVS was modelled via the time-dependent Stokes equations over this moving domain. Flow was induced in the PVS by different combinations of local and systemic effects including pulsatile arterial wall motions, pulsatile arterial rigid motions, and static and/or pulsatile pressure differences between the inlet and outlets (Fig 1C). We computed velocities and pressures in space and time, averaged normal velocities at the inlet and outlets over time, and net flow velocities (Fig 1D). In addition to the realistic geometry models, we considered an idealized PVS to study effects of domain length. A summary of the models considered is presented in Table 1.

PVS geometry and mesh generation
The PVS geometry was generated from image-based models of cerebral arteries (case id C0075) from the Aneurisk dataset repository [38]. The artery model was clipped to define a vessel segment of a healthy middle cerebral artery (MCA M1-M2) including an arterial bifurcation with one inlet vessel and two outlet vessels (Figs 1A and 2). The PVS domain was defined by creating a circular annular cylinder surrounding the artery with the arterial wall as To study the mechanisms behind perivascular fluid flow, we extracted an image-based bifurcating arterial geometry (A) and generated a computational model of a surrounding perivascular space (B) subjected to different forces: arterial wall deformations (red arrows), systemic pressure variations (blue arrows) and rigid motions (black arrows) (C) to predict the induced CSF flow and pressure (D its inner surface. Based on experimental observations of PVS width relative to the adjacent artery [2], we set the width of each PVS proportional to the arterial diameter. We then uniformly scaled the geometry down to the mouse scale: the PVS center line was then of maximal length 1 mm with inlet and outlet branches of comparable lengths (�0.5 mm), PVS widths of 28-42 μm and inner arterial diameters of 32-46 μm. We created a finite element mesh of the PVS with 174924 tetrahedrons and 34265 vertices using VMTK [39]. We also defined a set of idealized PVS domains as annular cylinders of lengths L 2 [1, 5, 10, 50, 100] mm with a annular cross-section width of 40μm. These annular cylinders were represented by one-dimensional axisymmetric finite element meshes with 10L + 1 vertices (mesh size 0.1 mm). Overview of the computational mesh generation. A) The artery geometry was extracted from the Aneurisk dataset repository [38] (case id C0075) and clipped. B) The domain center line was computed using VMTK, and subsequently used to define the extruded PVS. The color indicates the distance from the center line to the vessel wall. A finite element mesh was generated of the full geometry (including both the artery and the PVS) (C), before the outer PVS mesh (D) was extracted for simulations. https://doi.org/10.1371/journal.pone.0244442.g002

PLOS ONE
The mechanisms behind perivascular fluid flow

CSF flow model and parameters
To model the flow of CSF in surface PVS, we consider the CSF as an incompressible, viscous fluid flowing at low Reynolds numbers in a moving domain-represented by the time-dependent Stokes equations over a time-dependent domain. The time-dependent Stokes flow is a valid assumption due to the low Reynolds (Re < 0.01) and Womersley (α < 0.15) numbers. The initial PVS mesh defines the reference domain O 0 for the CSF with spatial coordinates X 2 O 0 . We assume that the PVS domain O t at time t > 0 has spatial coordinates x 2 O t and is defined as a deformation of the reference domain O 0 7 ! O t with x = d(X, t) for a prescribed space-and time-dependent domain deformation d with associated domain velocity w.
The fluid velocity v = v(x, t) for x 2 O t at time t and the fluid pressure p = p(x, t) then solves the following system of time-dependent partial differential equations (PDEs) [40]: The CSF density is set to ρ = 10 3 kg/m 3 and the dynamic viscosity μ = 0.697 × 10 −3 Pa/s [2]. On the inner PVS wall we set the CSF velocity v to match the given domain velocity w. We assume that the outer PVS wall is impermeable and rigid with zero CSF velocity. At the PVS inlet and outlet, we impose given pressures in the form of traction conditions. The system starts at rest and we solve for a number of flow cycles to reach a periodic steady state.

Pulsatile wall motion and velocity
We stipulate that arterial blood flow pulsations induce a pulsatile movement of the inner PVS boundary Λ. We let this boundary deform in the direction of the boundary normal with a spatially and temporally varying amplitude A: where n denotes the outward pointing boundary normal. For the amplitude A, we consider the combination of (i) the wall motions reported by Mestre [41]. To represent the spatial variation, we assume that the arterial pulse wave takes the form of a periodic travelling wave with wave speed c = 1 m/s [2] (and corresponding wave length λ = c/f for a given cardiac frequency f). We then set AðX; tÞ ¼ À 0:5 � dðt À kX À X 0 k=cÞ � R PVS ; dðtÞ ¼ 0:01Dd=dðsÞ; where s is the fraction of the cardiac cycle: s = (t � f)mod 1, with an average PVS width R PVS = 4.4 × 10 −2 mm, and X 0 is a fixed point close to the center of the inlet. We let the frequency vary between the different models (see Table 1).

Static and pulsatile pressure gradients
Pressure gradients in the arterial PVS can also be a consequence of the systemic phase-shift in pressure pulsations between larger proximal arteries, CSF, and the venous system [42], a general pressure increase in the CSF due to infusion [43], or other factors affecting the relative timing of the arterial and CSF pulse pressure [44]. To examine these systemic effects, we considered different static and pulsatile pressure gradients.
We associated static pressure gradients with the third circulation [45] (0.01 mmHg/m), the cardiac cycle (1.46 mmHg/m), and respiration (0.52 mmHg/m) [32]. The latter two values correspond to the peak amplitude of the pulsatile pressure gradients associated with these cycles [31,32], and should be considered as upper estimates of any associated static pressure gradients. During an infusion a change of at least 0.03 mmHg may occur between the CSF and the PVS [43]. The pressure drop occurs over at least a cortex thickness of 2.5 mm [46], and an upper estimate of the static pressure gradient due to infusion can thus be computed as dp ¼ 0:03 mmHg=2:5 mm ¼ 12 mmHg=m: ð3Þ Inspired by Bilston et al [44], we also considered a cardiac-induced pulsatile pressure gradient between the inlet and outlet with a phase shift relative to the pulsatile arterial wall movement of the form where f is the pulsatile frequency, and θ is the relative phase shift ranging from 0 to 1. The pressure gradients where weighted by the lengths of each branch to ensure that average pressure gradients from the inlet to the two different outlets are equal, and then applied as pressure differences between the inlet and outlets as traction boundary conditions (c 1 , c 2 , a 1 and a 2 in Table 1).

Modeling resistance and compliance
The PVS is part of a larger CSF system, and resistance and compliance far away from the local boundaries of the PVS model may affect flow in the PVS [23]. These global effects may be modeled as a resistance and compliance boundary condition according to a Windkessel model [47]. At the boundary of Model F (see Table 1) we thus solve to include the compliance C and resistance R of the brain. Here, p is the pressure at the boundary, and Q = R Γ u � nds is the volumetric flow rate (outflow) over the boundary Γ. Initially we set C = 1.798μL/mmHg and R = 1.097 mmHg/μL/min [23]. However, since the volume and volumetric flow rate of a single PVS (which we model) is much lower than the entire brain and CSF compartment, we also tested a compliance C = 0.001798μL/mmHg, i.e. 3 orders of magnitude lower than for the entire brain.

Image analysis of rigid motion and particle positions
We defined the arterial rigid motions by juxtaposing 28 screenshots, extracted at a fixed frequency, from [2, S2 Movie]. Comparing the artery outlines and motion, we estimated the peak amplitude γ of the motion to be no more than 6 μm ( Fig 5A) and identified a center point for the rigid motion X c close to the center of the bifurcation. We then defined the signed amplitude of the rigid motion as where δ min (resp. δ max ) is the minimum (resp. maximum) relative diameter variation, and kX c − X 0 k � 0.40 mm. We defined the (normalized) direction of the rigid motion r as normal to the main axis between X 0 and X c . We then investigated the impact of rigid motion on perivascular flow by imposing the following boundary displacement (in place of (2): The resulting rigid motions were used in two models as additional movement of the arterial wall (Model D and E listed in Table 1). From the screenshots, we also tracked the position of a number of sample microspheres over time using GIMP-Python [48].

Numerical solution
To compute numerical solutions of (1), we consider the Arbitrary Lagrange-Eulerian (ALE) formulation [40] with a first order implicit Euler scheme in time and a second-order finite element scheme in space. For each discrete time t k (k = 1, 2, . . .), we evaluate the boundary deformation d| Λ given by (2) and extend the deformation to the entire mesh by solving an auxiliary elliptic PDE. The computational mesh is deformed accordingly and thus represents O t k. We also evaluate the first order piecewise linear discrete mesh velocity w kÀ 1;k h associated with this mesh deformation.
Next, at each discrete time t k , we solve for the approximate CSF velocity v k h and pressure p k for finite element test functions ϕ h and q h defined on O t n, Δt is the time step size,p is the prescribed boundary pressure at the inlet and/or outlet Γ t n, and v nÀ 1 h and � nÀ 1 h are the approximate velocity and test function respectively at the previous discrete time t n−1 , defined over the previous domain O t n−1. We set Δt = 0.001 s.

Computation of output quantities
With the computed velocity field v, we define the flow rate at the inlet QðtÞ ¼ R G in vðx; tÞ � n dx, and the average normal velocity (see e.g. Fig 3G) was computed as v avg (t) = Q(t)/A in , where A in is the area of the inlet. From the average normal velocity, the position of a particle at time t was computed as Finally, the net flow velocity was computed as the slope between the peaks of x(t), using the two last cardiac cycles.

Computational verification
All numerical results were computed using the FEniCS finite element software suite [49]. Key output quantities were compared for a series of mesh resolutions and time step sizes to confirm the convergence of the computed solutions (S3 and S4 Figs). The simulation code, meshes and associated data are openly available [50].

Vascular wall pulsations induce oscillatory bi-directional flow patterns in the PVS
When inducing flow in the PVS by pulsatile arterial wall motions (Model A), the fluid in the PVS oscillated with the same frequency (10 Hz) and in-phase with the wall. During systole flow was bi-directional: the arterial radius rapidly increased, pushing fluid out of the domain at both the inlet (top) and the outlets (bottom) (Fig 3A). During diastole, flow was reversed (S1 Movie). Peak velocity magnitude occurred close to the inlet of the PVS model ( Fig 3A). From the inlet, the velocity magnitude decreased along the PVS until reaching a minimum close to the bifurcation. The velocity magnitude then increased towards the outlets but did not reach the same magnitude as at the inlet. The velocity profile in axial cross-sections of the PVS followed a Poiseuille-type flow pattern with high magnitude in central regions and low magnitude close to the walls (Fig 3A). The (average normal) velocity at the inlet was negative (downwards into the PVS) during diastole reaching close to -45 μm/s and positive (upwards, out of the PVS) during systole, reaching nearly 220 μm/s, giving a peak-to-peak amplitude of 265 μm/s (Fig 3F). While the pressure difference between the PVS inlet and outlets was set at zero, the pressure within the PVS again oscillated with the cardiac frequency and varied almost piecewise linearly throughout the PVS (Fig 3A). I.e. from the inlet, the pressure increased linearly to the point of peak pressure before it decreased linearly towards the outlet from there. The peak pressure was 0.38 Pa (0.0029 mmHg), and occurred in the smallest daughter vessel after the bifurcation (� 0.36 mm from the left outlet). The time of peak pressure coincided with the time of peak velocity (t = 0.048 s). At time of peak pressure, the pressure gradient magnitude was nearly uniform throughout the domain, with an average gradient magnitude of 6.31 mmHg/m. High pressure gradients were observed locally in a small narrow region of the inlet vessel reaching a maximum gradient magnitude of 93.0 mmHg/m.

Flow is slow and laminar around the bifurcation
In all models without a static pressure gradient, flow around the bifurcation followed a similar pattern as in the rest of the domain. The flow was slow and laminar, with reversal of flow direction going from systole to diastole (Fig 4). PVS fluid velocities in the bifurcation region typically reached 25 μm/s during systole and 5 μm/s during diastole (Fig 4B and 4D). The primary reason for lower velocities in the bifurcation region was the central placement of the bifurcation within the domain. No recirculation regions or circular flow were observed around the bifurcation despite the fact that bidirectional flow occurred in these regions (Fig 4).

Static pressure gradients induce net PVS flow
Static CSF pressure gradients may occur as a direct consequence of tracer infusions [43]. However, such gradients also occur naturally due to e.g. the third circulation as a pressure gradient driving slow steady flow from the choroid plexus to the arachnoid granulations [32,45]. To assess the static systemic effect and pulsatile local effect, we simulated flow and pressure under a static pressure difference between the inlet and outlets of the PVS model in combination with the pulsatile wall motion (Model B). Several pressure gradients were examined, representing forces involved in the cardiac and respiratory cycle, the third circulation, and in infusion tests [32,43].
In general, the additional static pressure gradient induced net flow in the downwards direction, with the presence of oscillatory flow including backflow depending on the magnitude of the applied gradient and location. With a static pressure gradient of 1.46 mmHg/m (Fig 3B), which is representative of the pulsatile pressure gradient induced by the cardiac cycle [32], net flow velocities varied between 20 and 30 μm/s depending on the location of measurement: at the inlet, the net flow velocity was 28 μm/s. The peak-to-peak amplitude of the average normal velocity was unchanged at 265 μm/s (Fig 3G). A particle suspended at the PVS inlet would then experience a pulsatile back-and-forth motion with a net movement downstream: 1-2 μm upstream during systole and 5 μm downstream during diastole (Fig 3H). At the outlets, the average normal velocities were lower (in absolute value) due to the larger area and flow was nearly stagnant during diastole (backflow was negligible) (S1 Fig). The static pressure gradient did otherwise not change the shape of the velocity pulse. We note that the pulsatile motion of particles was not easily visible, except for particles close to the inlet or outlets (Fig 3B).
The gradient induced by steady production of CSF i.e. corresponding to the third circulation (0.01 mmHg/m [32]) generated negligible net flow, while a gradient equal to the respiratory gradient (0.52 mmHg/m) gave about threefold lower net flow velocity than the cardiac gradient described in detail above. An upper estimate of a gradient induced by infusion (� 12 mmHg/m) resulted in net flow velocity of more than 100 μm/s.

Flow induced by blood and CSF asynchrony
Flow in cranial or spinal PVS may be influenced by the relative timing of pulsatile blood and CSF pressures [44]. With physiological pressure gradients, we investigated (Model C) to what extent differences in phase-between the pulsatile arterial wall motion and some pulsatile systemic pressure gradient-would induce fluid velocities and net flow in the PVS [2,3].
A pulsatile cardiac pressure gradient with peak amplitude 1.46 mmHg/m (=0.195 Pa/mm), frequency 10 Hz, combined with the pulsatile arterial wall motion at a phase shift θ, again induced laminar flow in the PVS with streamlines taking the shortest path between the inlet and outlets (Fig 3C). Further, the velocity profiles were comparable to those induced by the arterial wall motion alone (Fig 3F-3H, S1 Fig). The velocity profile resulting from a phase shift of 10% of the cardiac cycle, while still similar, differed the most from the previous experiments. Qualitative differences could be observed during diastole, and the peak-to-peak amplitude of the velocity was slightly reduced at 262 μm/s. However, no substantial net movement of fluid in any direction was observed: a particle suspended at the inlet of this systemic model would oscillate back-and-forth with a peak-to-peak amplitude of 2 μm (Fig 3H). Regardless of the phase shift applied in the systemic pressure gradient, net flow velocity did not exceed 0.5 μm/s (S2 Fig).

Arterial rigid motions induce complex flow patterns
In addition to its pulsatile expansions and contractions, an artery can undergo pulsatile rigid motions i.e. rotations or local translations, possibly independent from movement of the rest of the body. The potential effect of such arterial rigid motions on PVS flow as well as the underlying causes of such movements are poorly understood. To investigate, we extracted experimentally observed arterial rigid motions [2, S2 Movie] and simulated how this additional movement could affect flow in the PVS (Model D, Fig 5A).
The rigid motions extracted were synchronous with the arterial wall pulsations. When combined with these pulsations, the rigid motion of the artery increased fluid motion within the PVS (Fig 3D). As the artery shifted, the displaced fluid tended to move to the other side of the artery, thus yielding more complex streamline patterns and swirls. The peak-to-peak velocity amplitude was 251 μm/s, which is slightly lower than for the other models. However, the rigid motion did not affect the overall shape of the velocity pulse at the inlet or outlets (Fig 3F and  3G, S1 Fig) or the net flow velocity (Fig 3H). The rigid motion resulted in more complex flow, which will enhance local mixing.

Arterial pulsation frequency modulates PVS flow velocity
The typical duration of the mouse cardiac cycle has been reported as 80-110 ms [51], corresponding to a cardiac frequency of 9-12.5 Hz. However, experimental studies of perivascular flow also also reveal reduced heart rates as low as 2.2 Hz [2, S2 Movie]. Reducing the frequency of the arterial wall pulsations from 10 to 2.2 Hz (Model E) reduced the peak velocity by a similar factor: the peak-to-peak amplitude of the average normal velocity at the inlet was reduced from 260 to 60 μm/s. Adding a static gradient of 1.46 mmHg/m to the 2.2 Hz pulsations again induced net flow velocities of 20-30 μm, but the pulsatile motion of particles was small compared to net flow (S2 Movie). Combining rigid motions with the experimentally observed arterial wall pulsation frequency of 2.2 Hz and a static pressure gradient of 1.46 mmHg/m, induced oscillatory PVS flow with non-trivial flow patterns, backflow, and net flow (S3 Movie, Fig 5B). The pulsatile rigid motions induced oscillatory particle movement normal to the arterial wall, while the pulsatile arterial expansion induced back-and-forth movement along the PVS. Superimposed on the steady downwards flow induced by the static gradient, the movement of particles were thus similar to existing experimental observations [2] both in terms of net flow velocities and peak-to-peak pulsation amplitude.

Resistance and compliance do not increase net flow
Adding compliance and resistance at the outlets (Model F) did not change the pressure, and the outlet pressure thus remained close to 0 during the entire cardiac cycle. Consequently, the flow rate did not change. As in the previous models without a pressure drop, the pulsations pushed fluid out of each end of the PVS during systole, and during diastole the PVS refilled from both sides (data not shown). Reducing the compliance by three orders of magnitude allowed for less fluid to leave the PVS before the pressure increased in response. The pressure at the two outlets were similar, but differed slightly due to different outflow rates (Fig 6). The resistance prevents most flow at the outlets during diastole (Fig 6). In the initial phase of systole, the arterial expansion results in a large pressure in the central PVS, causing fluid to leave the domain at both ends, but much slower at the outlets (Fig 6). As more fluid leaves the PVS at the outlets, the compliance comes into play due to the accumulated volume going out. Therefore, the flow direction changes earlier at the outlets than at the inlet. Peak pressure The pressure is similar at each outlet, but differ in time of peak and peak value. At the inlet, the pressure is always close to 0. C) Fluid velocity at the inlet and outlets. Compliance and resistance at the outlets restrict flow over these boundaries. The inflow velocity is more than four times larger than the outflow velocities. Peak velocity occurs earlier at the outlets than at the inlet, as the pressure at the outlets increase as fluid starts to move out. Net flow is still negligible (net flow velocity <0.2μm/s).
https://doi.org/10.1371/journal.pone.0244442.g006 occurs at slightly different time points for each outlet. At the inlet, the pressure is always close to 0 (zero traction condition). Net flow did not change by adding resistance and compliance to the PVS model and were still negligible (net flow velocity <0.2μm/s).

Model length modulates PVS velocity and net flow
Mathematical modeling of PVS flow in idealized geometries has demonstrated that, under certain conditions, peristaltic motion of the arterial walls could induce substantial net flow velocities [25,36]. However, these findings have not been supported by computational models [52]. The mathematical model [36] represents the PVS as an infinitely long annular cylinder, while in-vivo and in relevant computational models, the PVS is considerably shorter than the wavelength of the arterial pulse wave [52]. To examine the effect of model length on PVS velocities and net flow, we considered a idealized axisymmetric model of an annular cylinder (Model G) of different lengths L (1, 5, 10, 50 and 100 mm) for a fixed frequency (10 Hz) and arterial pulse wavelength λ (100 mm).
When the model length is shorter than the wavelength, velocities are highly dependent on the length of the PVS (Fig 7). For L � λ, the wall displacement is close to uniform along the PVS, and more fluid will leave the domain through the inlet and outlet in a longer artery ( Fig  7C). Thus, for a given relative wall displacement and model lengths smaller than half the arterial pulse wavelength, the velocity at the inlet (or outlet) will increase with increasing PVS model length (Fig 7A). The shape of the velocity pulse also change: for longer models, at the inlet, we observe a longer period of upwards flow (out of the domain) and a corresponding shorter period of downwards flow (into the domain). When the domain length is equal to an integer multiple of the wavelength, using a zero pressure drop or a symmetry boundary condition will model an infinitely long cylinder. For this case, the velocity will not increase further with increasing PVS length. To obtain a peak-to-peak velocity amplitude of �20μm/s [2], a frequency of 10 Hz required a PVS length of 0.10 mm. Changing the frequency to 2.2 Hz, required a PVS length of 0.47 mm to reach the same amplitude.
For the same set of geometries, net flow also increased with model lengths up to the wavelength, and in the long idealized models of lengths 50 and 100 mm, net flow velocities of 4.7 and 7.0 μm/s were observed (Fig 7B). For the other PVS lengths tested, the net flow velocity was lower than 1 μm/s. Net flow velocities were small compared to the large average normal velocity amplitudes: a particle suspended at the PVS inlet could experience a change in position of up to 60 μm over one cardiac cycle (10 Hz, model length 50 mm) (Fig 7B).

Discussion
Experimental studies of perivascular flow have found substantial velocities and net particle movement, predominantly in a uni-directional pattern. Mestre et al [2] reported a peak-topeak velocity amplitude of �20μm/s and a typical net flow velocity of 18.7 μm/s. With a period of 0.45 s (2.2 Hz), a particle could then be expected to move no more than 9 μm back-andforth per cycle. Similarly, Bedussi et al [3] reported an average net flow velocity of 17 μm/s and a mean amplitude of movement of 14 μm. However, the shorter cardiac period of 0.15 s points at much higher velocity amplitudes (at least 100-200 μm/s) in the latter study. Our peak-topeak velocity amplitudes (251-265 μm/s for 10 Hz, 60 μm/s for 2.2 Hz) are thus at the upper range of experimental values reported. Our observations further point at the impact of cardiac frequency on velocity amplitude, which may explain the difference in estimated velocity amplitudes between these two experimental studies. On the other hand, a particle suspended at the PVS inlet in our model would oscillate back-and-forth with a peak-to-peak amplitude of 2 μm, with an additional net downstream movement only if a static pressure gradient is imposed. These changes in position are thus at the lower end of the experimental observations, pointing at the likely presence of a static CSF pressure gradient in the experimental configurations.
Waste products such as amyloid beta have been reported to concentrate distal to bifurcations [16]. In the bifurcation region in our model, the pulsatile motion of particles was substantially slower than in the rest of the PVS. However, in the presence of a static pressure gradient, net flow velocities responsible for particle movement were of similar magnitude all along the PVS. The difference in pulsatile motion between the bifurcation and the rest of the PVS is caused by the central placement of the bifurcation in our domain. As such, low pulsatile velocities near the bifurcation found in our model can not be associated with accumulation of particles in this region.
A static CSF pressure gradient-of magnitude corresponding to the pulsatile gradient induced by the cardiac cycle-was sufficient to create net flow velocities of 30-40 μm/s in the PVS. A pressure gradient associated with the third circulation [32] was not sufficient to drive net fluid movement in the PVS. The respiratory gradient is approximately one third of the cardiac gradient [32], and was sufficient to drive some PVS flow when applied as a static pressure difference. Thus, longer waves (such as those induced by respiration and vasomotion) may also play a role in net fluid movement in the PVS. Our upper estimate of the pressure gradient induced between the CSF and the PVS during an infusion test (in humans) [43] resulted in net flow velocities of several hundred μm/s, much higher than what has been observed in mice [2,3]. Overall, our observations indicate that the static pressure gradient sufficient to drive net flow is small (�1.5 mmHg/m, i.e. a pressure difference of 0.015 mmHg per cm) compared to the intracranial pressure increase of 1.4-3 mmHg observed in mice during tracer infusion [33,53].
Arterial rigid motions were of greater amplitude than the arterial wall expansions, but had minimal effect on average and net flow velocities. Indeed, the rigid motions of the artery did not force fluid to leave the PVS domain, but rather displaced fluid within. As such, arterial rigid motions can create oscillatory movement of cardiac frequency within the PVS without adding to the net movement of particles. In our model, rigid motions and arterial expansion combined can explain oscillatory motion as seen by Mestre et al. [2], but were not sufficient to generate net flow. However, the rigid motion introduced complex swirling that would significantly enhance local mixing and potentially contribute to increased clearance.
A systemic CSF pulsation of physiological amplitude [32,43]-out of synchrony with the arterial wall pulsation-did not induce net fluid movement in the brain PVS. The relative timing of arterial and CSF pulse waves, a possible factor for net fluid movement in spinal cord PVS [44], is thus not a likely factor for explaining higher average or net flow velocities in brain PVS. Adding resistance and compliance at the PVS outlets suppressed pulsatile flow amplitudes and resulted in a phase shift of the outlet flow, but did not affect net fluid movement. These effects of resistance and compliance in our model are thus consistent with a recent report by Ladrón-de-Guevara et al. [23] The length of the PVS segment is important for the observed fluid dynamics in the domain, and is also an important modeling parameter. When the PVS model is much shorter than the wavelength of the arterial pulse wave, the velocity amplitudes at the inlet and outlets of the PVS are directly linked to the PVS length. Similar effects have previously been noted by Asgari et al. [21] as an increase in dispersion effects with length, and by Rey and Sarntinoranont [37] as a variable Péclet number throughout the PVS domain. Considerable net flow in our models were seen only for very long geometries.
Several modelling studies have now tried to explain net movement of fluid in surface and parenchymal PVS driven by local arterial pulsations. While theoretical considerations have explained net flow by arterial wall pulsations alone [25,36], most computational studies [21,35,37,52] suggest that the local effect of arterial wall pulsations is not sufficient to drive net flow in the PVS of magnitude comparable to experimental observations at the pial surface. Theoretical considerations have assumed an infinitely long cylinder, which we here show overestimates net flow and velocity amplitudes compared to PVS models of physiologically relevant lengths. In our idealized computational models however, net flow velocity is higher than the theoretical model by Wang and Olbricht [36] predicts: with a 0.7% half-amplitude of the arterial wall expansion, inner and outer radii of 20 and 60 μm, and an arterial wave speed of c = 1 m/s, their model predicts an average net flow velocity of 1.53 μm/s, which is lower than our estimates of 6.7-7 μm/s by a factor of four. It should be emphasized that while the Wang and Olbricht net flow model [36] does not differentiate between PVS lengths, it is sensitive to PVS width and half-amplitude; thus specific such parameters could yield a higher net flow velocity. It should also be noted that the model used by Wang and Olbricht has a greater hydraulic resistance due to the porous media assumption, which may explain the lower flow magnitude compared to our model.
In terms of limitations, PVS flow was modelled as incompressible viscous fluid flow ignoring potential barriers to flow (reduced permeability) in contrast to e.g. [24,36]. For pial surface PVS, this may be a reasonable assumption. The addition of a finite permeability would be expected to yield lower velocities. We also ignored nonlinear (turbulent) effects, which, given the low Reynolds numbers involved, seems a reasonable approximation. The PVS domain was assumed to have a constant cross-sectional width, while other studies have suggested that the PVS cross-section is elliptic [2,54]. An elliptic PVS surrounding a circular artery may increase PVS flow up to a factor of two compared to the circular annulus used in our study [54]. For the representation of the rigid motion, we estimated its magnitude (�6μm) without isolating the �1 μm wall pulsations. The findings reported here thus represents an upper estimate of the impact of arterial rigid motions. Finally, all pressure gradients used to drive flow in our models originated from human measurements, while both the PVS size and previously reported velocities were obtained in mice. Compared to humans, mice have smaller CSF volumes, lower ICP and ICP amplitudes and a shorter CSF turnover time [55]. The latter suggests that the third circulation gradient is larger in mice than in humans. However, an increase of a factor � 50 from human to mouse would be needed for the third circulation gradient to drive any substantial net flow. We mainly considered static pressure gradients, but other non-zero cardiac cycle averaged pressure gradients could yield similar results.
In conclusion, our simulations indicate that the combination of arterial wall pulsations and a systemic static pressure gradient larger than that associated with the third circulation can explain experimental findings on pulsatile perivascular flow. The required static gradient need not necessarily be caused mainly by an infusion, as such a gradient is at the order of physiological pressure gradients in the brain. Without a pressure gradient, net flow was only achieved for very long PVS geometries (on the order of the wavelength, here 100 mm), explaining why theoretical considerations of infinitely long cylinders yield net flow. Finally, rigid arterial motion can induce complex flow patterns in the PVS. Cardiac frequency: 10 Hz. The outward pointing normal at the inlet defines the positive direction; negative values thus refer to flow downwards (into the PVS at the inlet, and out of the PVS at the outlet). All models predict bi-directional flow during systole, with fluid leaving the domain at both inlet and outlets. The peak velocity amplitude is slightly higher at the inlet during systole mainly due to the smaller area for flow compared to the combined area at the outlets. The velocity during diastole differ more between the models, and more (in absolute value) between at the inlet and outlet. . Key output quantities converge as the time resolution is reduced as expected. For a PVS length of 100 mm, the time step required for convergence was lower than for the shorter models. C) Peaks of the position over time for different time resolutions (L = 100 mm). We note that the computed net flow velocity strongly depends on the time resolution, but that a time step of 1 ms is sufficient. D) Position over time for the bifurcating arterial geometry (C0075) for different time resolutions. The time step of 1 ms is again sufficient. (TIF)

S4 Fig. Numerical verification and convergence analysis of the computational models. A)
Flow rate in the 2D model with different mesh resolutions. B) Close-up for a very short time period in A) to emphasize the small differences between the different mesh resolutions in 2D. C) Very small differences were also observed in 3D where two meshes were tested. Peak velocity obtained in the two meshes differed by � 2.5%. (TIF)