==== Front
Sci Rep
Sci Rep
Scientific Reports
2045-2322 Nature Publishing Group UK London 

77787
10.1038/s41598-020-77787-4
Article
Peristaltic flow in the glymphatic system
Romanò Francesco francesco.romano@ensam.eu 1 Suresh Vinod 2 Galie Peter A. 3 Grotberg James B. 4 1 grid.503422.20000 0001 2242 6780Univ. Lille, CNRS, ONERA, Arts et Métiers Institute of Technology, Centrale Lille, UMR 9014 - LMFL - Laboratoire de Mécanique des Fluides de Lille - Kampé de Fériet, 59000 Lille, France 
2 grid.9654.e0000 0004 0372 3343Auckland Bioeng. Inst. and Dept. Eng. Sci., University of Auckland, 70 Symonds Street, Bldg 439, Auckland, 1010 New Zealand 
3 grid.262671.60000 0000 8828 4546Dept. Biomed. Eng., Rowan University, 201 Mullica Hill Rd, Glassboro, NJ 08028 USA 
4 grid.214458.e0000000086837370Dept. Biomed. Eng., University of Michigan, 2123 Carl A. Gerstacker Building, 2200 Bonisteel Boulevard, Ann Arbor, MI 48109-2099 USA 
3 12 2020 
3 12 2020 
2020 
10 2106510 9 2020 12 11 2020 © The Author(s) 2020Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.The flow inside the perivascular space (PVS) is modeled using a first-principles approach in order to investigate how the cerebrospinal fluid (CSF) enters the brain through a permeable layer of glial cells. Lubrication theory is employed to deal with the flow in the thin annular gap of the perivascular space between an impermeable artery and the brain tissue. The artery has an imposed peristaltic deformation and the deformable brain tissue is modeled by means of an elastic Hooke’s law. The perivascular flow model is solved numerically, discovering that the peristaltic wave induces a steady streaming to/from the brain which strongly depends on the rigidity and the permeability of the brain tissue. A detailed quantification of the through flow across the glial boundary is obtained for a large parameter space of physiologically relevant conditions. The parameters include the elasticity and permeability of the brain, the curvature of the artery, its length and the amplitude of the peristaltic wave. A steady streaming component of the through flow due to the peristaltic wave is characterized by an in-depth physical analysis and the velocity across the glial layer is found to flow from and to the PVS, depending on the elasticity and permeability of the brain. The through CSF flow velocity is quantified to be of the order of micrometers per seconds.

Subject terms
Fluid dynamicsComputational modelsBiomedical engineeringUniversity of Auckland Vice Chancellor’s Distinguished Visitor Awardissue-copyright-statement© The Author(s) 2020
==== Body
Introduction
Cerebrospinal fluid serves as a sink for metabolic waste products generated in the brain. The pathway for CSF transport in the brain interstitium has been a puzzle. Recent imaging experiments using in vivo two-photon microscopy have lent support to the hypothesis that CSF enters the brain from the subarachnoid space along the perivascular sheaths surrounding penetrating arteries and ‘leaks’ out into the interstitium across a permeable layer of glial (astrocyte) cells. From there, it is cleared into the perivascular sheaths around veins and the pulsation of the cerebral arteries are identified as an important driver for the transport of CSF into the brain interstitium1,2. Since convective bulk flow of the CSF between these ingress and egress pathways facilitates the clearance of solutes and metabolic waste products from brain tissue, dysfunctions in CSF transport may have implication for a range of neurological conditions such as intracranial hypertension and protein clearance in Alzheimer’s disease and Parkinson’s disease. Empirical studies2–4 indicate that CSF transport is affected by the elastic properties of vessel walls, water permeability of brain tissue and pulsatility of blood flow. However, the difficulty of measuring these parameters in vivo necessitates modeling-based approaches to improve our understanding of fluid transport in the brain. Therefore, the aim of this study is to develop a mathematical model of perivascular transport that provides insight into how these factors alter the direction and magnitude of CSF flow. Since we aim at deriving a leading-order characterization of the CSF flow, the impact of ciliated boundaries and non-Newtonian effects5–7 will be neglected in our model.


The model described here builds upon previous approaches to calculate perivascular fluid flow in idealized geometries. Wang and Olbricht8 studied axial flow and transport in an annulus with impermeable boundaries, but did not address fluid exchange with the interstitium9. Kyrtsos and Baras modeled protein clearance from the interstitium using a compartmental model in which CSF velocity was an input parameter and was assumed to be inversely proportional to vessel stiffness10. Cerebral MRI visualizations of a live rat have been used by Ratner et al.11 to find the direction of the glymphatic flow. Moreover, they made use of a purely diffusive model to estimate the liquid flow through the healthy brain of a rat and reproduced the main experimental features by means of an Optimal Mass Transfer approach which could also estimate the diffusion tensor based on the dynamic flow rate. By means of numerical simulations, Asgari et al.12 claimed that the arterial pulsation due to the peristaltic wave cannot be the tribological driving force responsible of the interstitial solute transport. They address the role of dispersion transport, which is a combination of local mixing and diffusive effects in the para-arterial space. A very different conclusion has been drawn by Aldea et al.13, who proposed a multiscale model of the arteries dealing with the basement membrane as a deformable fluid-filled, poroelastic medium. They rather concluded that the vasomotion-driven intramural periarterial drainage is compatible with experimental observations.

Jin et al.14 modeled the glymphatic system from para-arterial to paravenous cerebrospinal fluid through brain extracellular space. They investigated the glymphatic mechanism for solute clearance in brain by modeling diffusive and convective transport in the cerebral extracellular space, focusing on the short-range transport between para-arterial and paravenous spaces. Based on the numerical simulations of their model, they concluded that the convective transport is not affected by the pressure fluctuations and requires a strong pressure gradient to be significant. Moreover, they found that the convective transport is also fairly insensitive to astrocyte endfoot water permeability and that diffusion transport suffices to explain the experimental data of the transport studies in brain parenchyma. Similarly, Faghih and Sharp15 use a one-dimensional steady, pressure-driven branching flow model to analyze the hydraulic resistance of arterial membranes. They found that the resistance of the periarterial tree is too great to account for physiological estimates of the CSF leakage rate, and that a combined route through the paraarterial and paravenous spaces would also be unlikely based on the magnitude of the transmantle pressure. A similar approach was employed by Rey and Sarntinoranont16, who made use of two resistance network models to study the effect of pulsating flows. They estimated that the peak fluid velocity in the PVS and parenchyma increases with the pulse amplitude and the vessel size, making the convective solute transport less and less relevant.

Our model derives a leading order approximation of the Navier–Stokes equation which is based on lubrication theory and includes the effect of a peristaltic wave in the artery and the deformability of the brain tissue. A further justification of the negligible importance of convective transport compared to diffusion effects in the PVS will be derived from first principles, motivating such a conclusion by dimensional analysis considerations. Thereafter, focusing on CSF exchange between the perivascular space and brain interstitium, we compute the leak velocity using our first-principles approach.

Model
Geometry
The CSF-filled perivascular space was modeled as a thin annular gap between an elastic, impermeable artery and a brain tissue (Fig. 1). An elastic, permeable glial boundary separates the PVS from the brain tissue. The glial boundary is modeled with a linear elastic wall law and we do not solve for the flow in the artery. Instead, the peristaltic wave deformation of the artery is prescribed as a travelling wave. The interstitial pressure was assumed to be constant and was used as the reference pressure. Linear elastic tube laws were used to relate the deformations of the solid boundaries to the pressure difference across them. Governing equations were simplified using lubrication theory.Figure 1 Sketch of the perivascular space between the brain tissue and the artery.



The thickness of the PVS is b, the average radius of the artery and the average inner radius of the brain tissue are r1 and r2, respectively. The deformation about the average radii are h(z,t)=h¯sin2π/λz-ct and d(z, t), where h is the imposed deformation of the travelling peristaltic wave, h¯, λ and c its amplitude, wavelength and velocity, respectively, z is the traveling (axial) direction, t the time, and d the glial boundary deformation. The pressure in the brain tissue is pe, whereas pa and pb<pa are assumed as pressures at the extrema of the PVS.

Governing equations
The Navier–Stokes and continuity equations in dimensional form read  1a ∂u∂t+u∂u∂r+w∂u∂z=-1ρ∂p∂r+μρ1r∂u∂r+∂2u∂r2+∂2u∂z2,  1b ∂w∂t+u∂w∂r+w∂w∂z=-1ρ∂p∂z+μρ1r∂w∂r+∂2w∂r2+∂2w∂z2,  1c 1r∂(ru)∂r+∂w∂z=0,  where u→=(u,w) denotes the velocity field in cylindrical coordinates, p is the pressure, r→=(r,z) and t are the spatial and temporal coordinates, respectively, ρ is the density of the fluid flowing in the perivascular space, μ its dynamic viscosity.

Equation (1) are then scaled with 2 r=bR,z=λZ,t=Tω,u=bωU,w=λωW,p=μωε2P, where ω=c/λ is the frequency of peristaltic wave in the artery, b is the thickness of the perivascular fluid film, and ε=b/λ. The non-dimensional continuity and Navier–Stokes equations read  3a ε3Re∂U∂T+U∂U∂R+W∂U∂Z=-∂P∂R+ε21R∂U∂R+∂2U∂R2+ε2∂2U∂Z2,  3b εRe∂W∂T+U∂W∂R+W∂W∂Z=-∂P∂Z+1R∂W∂R+∂2W∂R2+ε2∂2W∂Z2,  3c 1R∂(RU)∂R+∂W∂Z=0,  where Re=ερωb2/μ is the Reynolds number.

The mathematical problem (3) is closed by the boundary conditions  4a Z=0:P=Pa,  4b Z=L:P=Pb,  4c R=R1+H:U=∂TH,W=0,  4d R=R1+1+D:U→·n→=∂TD+MeP-Pe,  4e U→·t→=0,  4f D=P-Pe/Ee,  where Me=kgμ/bε2, Pe is the non-dimensional pressure at the glial boundary and kg the permeability of the brain tissue, Ee=Egε2/μω(R1+1) and Eg the Young’s modulus of the brain tissue. H=h/b=
H¯sin2πZ-T is the imposed deformation of the peristaltic wave of the artery and H¯=h¯/b its amplitude measured from the middle line R1=r1/b. D=d/b is the deformation of the peristaltic layer at the boundary with the brain tissue; it is measured from the middle line R2=r2/b=R1+1. The normal and tangent unit vector to the glial boundary are denoted by n→ and t→.

Thin film approximation
Assuming that ε≪1, i.e. the wavelenth λ of the peristaltic wave is much larger than the film thickness b, we expand the pressure, the velocity field and the deformation D with the following polynomial series  5a P=P0+εP1+ε2P2+…,  5b U=U0+εU1+ε2U2+…,  5c W=W0+εW1+ε2W2+…,  5d D=D0+εD1+ε2D2+…,  where the subscript 0 denotes the solution at the leading order term, O(ε0), 1 refers to the linear correction in ε, O(ε1), 2 indicates the quadratic correction in ε, O(ε2), and so on.

If Re=O(1) or smaller, the leading order continuity and Navier–Stokes equations read  6a ∂P0∂R=0,  6b ∂P0∂Z=1R∂W0∂R+∂2W0∂R2,  6c 1R∂(RU0)∂R+∂W0∂Z=0.  The system (23) is completed by the boundary conditions at leading order  7a Z=0:P0=Pa,  7b Z=L:P0=Pb,  7c R=R1+H:U0=∂TH,W0=0,  7d R=R1+1+D0:n→=(nr,nz)=(1,0),t→=(tr,tz)=(0,1),  7e U0=MeP0-Pe+1Ee∂P0-Pe∂T,  7f W0=0,  7g D0=P0-PeEe.  More details of the model derivation are given in the section Methods. We however remark that leading order boundary conditions do not include any axial flow at the walls, i.e. W0=0 along the boundaries. For more details about the effect of the axial flow in the glymphatic system, we refer to Albargothy et al.17.

The leading order term of the r-momentum implies that P0 is only function of Z and T and, integrating in r the z-momentum, W0 reads 8 W0=R24∂P0∂Z+C1ln(R)+C2, where C1=C1(Z,T) and C2=C2(Z,T).

Plugging (8) in the continuity equation and integrating in R, we derive the form of U0 9 U0=C3R-R316∂2P0∂Z2-R42ln(R)-1∂C1∂Z-R2∂C2∂Z, where C3=C3(Z,T).

Applying the no-slip boundary conditions at the inner boundary (R=R1+H), i.e. U0=∂TH and W0=0, and the permeable boundary condition at the outer boundary (R=R1+1+D), i.e. U0=Me(P0-Pe)+Ee-1∂T(P0-Pe) and W0=0, the functions C1, C2 and C3 are expressed in terms of P0 and a second-order differential equation is derived for P0 10 ∂P0∂T+A0∂2P0∂Z2+A1∂P0∂Z+A2P0=Ee(R1+H)R1+1+D0∂H∂T+EeMePe+∂Pe∂T. We refer to the section Methods for the definition of A0, A1, A2. The initial condition for P0 is set to be the linear function consistent with the boundary conditions 11 P0(T=0)=Pa(T=0)+Pb(T=0)-Pa(T=0)Z/L. Equation (10) is solved numerically by making use of a collocation spectral method in Z-direction and discretizing in time by implicit Euler method. The time step employed to discretize in time is always set equal to Δt=10-2 and 1000 Chebyshev–Gauss–Lobatto nodes are used in Z. Further details about the numerical method we chose and its implementation in our solver are reported in the section Methods. In order to avoid non-linearities in the discretization algorithm, the explicit outer-wall deformation D0n is employed when computing the solution at time tn+1.

We remark that taking into account the recent study by Mestre et al.4, the annular space around blood vessels in the brain is not uniform in thickness. In this sense, our axisymmetric approach represents a simplification of the perivascular space geometry assuming that the averaged cross-sectional radius is sufficient to capture the leading-order effects of the CSF dynamics. We stress that such an assumption is at the core of the simplified one-dimensional time-dependent partial differential Eq. (10). Including a non-uniform deformation of the gap cross-section would require a non-trivial extension of the model as the pressure, the glial boundary deformation and the flow velocity should depend by all the three spatial coordinates, i.e. R, Z and Φ. The resulting system of PDEs would make the extended model far more complex to solve and we expect it would not provide major improvements in the order-of-magnitude estimate of the steady streaming across the glial boundary.

Results
Physiological parameters
Table 1 Range of the non-dimensional groups for the thin-film problem between an artery and a brain tissue.

Parameter	Description (definition)	Estimated range	
Re	Reynolds number ερωb2/μ	[5×10-12,2×10-4]	
ε	Perivascular film thickness b/λ	[10-6,3.75×10-2]	
R1	Inner radius of the perivascular layer r1/b	[0.7,103]	
H¯	Amplitude of the peristaltic wave h¯/b	[8×10-4,0.5]	
L	Perivascular length l/λ	[2,20]	
Me	Permeability of the brain tissue kgμ/bε2	[10-2,10]	
Ee	Stiffness of the brain tissue Egε2/μω(R1+1)	[2×10-10,0.2]	


Physiologically relevant parameters for the thin liquid film of interest are derived from the literature. Xie et al.18 estimates the intracranial pressure to be about 2×103  Pa, in accordance with Sakka et al.19, who report pe∈[1300,2000]  Pa. Assuming Pe=Pb=0 as reference pressure in our model, Pa will then be considered in the range Pa∈[5×10-7,5×102]. Wang and Olbricht8 report the values measured in previous studies for: the peristaltic wave frequency ω≈5  Hz3, the inner radius of the perivascular layer r1≈10-5  m20, the outer radius of the perivascular layer r2≈1.1×10-5  m21 (which result in a perivascular film thickness b≈10-6  m), the peristaltic wave amplitude h¯∈[1.25,5]×10-7  m22, the peristaltic wave speed c≈1  m s-120,23,24 (which result in a peristaltic wavelength λ≈0.2  m), and the dynamic viscosity μ=9×10-4  Pa s25. Plus, the fluid density is comparable to the one of water ρ≈ρwater=103  kg m-3. Other values of the artery radius (r1∈[10-4,10-3]  m) are reported in Thorin-Trescases et al.26 Moreover, considering that the artery wall thickness ta is 10 to 100 times smaller than r1, i.e. ta∈[10-6,10-5]  m27, that the elastic modulus of the artery wall is Ea∈[105,106]  N m-2 and that the wall density is ρa≈103  kg m-3, the elastic wave speed is ca=Eata/2r1ρa∈[0.2,7]  m s-1. Considering the relationship given by Atabek28, i.e. c∈[0.1,0.5]ca, an estimated range of the wave speed of the peristaltic wave can be proposed, which includes the estimate of Wang and Olbrich8: c∈[0.02,3.5]  m s-1 resulting in a range for the peristaltic wavelength λ∈[0.004,0.7]  m, which includes the previous estimate λ≈0.2  m. Another estimate of b is given by Iliff et al.1, who reports b=10-5  m, and by Jin et al.14, who reports b=1.5×10-4  m. Finally, The elastic modulus and the permeability of the brain tissue are Eg≈0.01Ea∈[103,104]  Pa and kg≈10-11  m  Pa-1  s-1, respectively29,30. Based on these parameters, the range of the non-dimensional groups of interest for our study is derived and reported in Table 1. From Table 1 it is clear that the leading-order thin-film approach is a very good approximation for our problem since ε≪1 and Re≪1. Higher-order corrections in ε and inertial effects due to Re are, therefore, substantially negligible.

Parametric study
Based on the physiologically relevant parameters, we carried out numerical simulations for the following range of non-dimensional groups 12 R1∈[10,1000],H¯∈[0,0.2],L∈[2,20],Me∈[0.1,5],Ee∈[0.01,1]. We remark that the range of the glial boundary elasticity parameter Ee has been restricted to vary over 3 rather than 10 orders of magnitude, as indicated in Table 1. Indeed, when the Young modulus of the brain tissue is very small, i.e for the softest brain tissue parameters reported in Table 1, Ee≪1, a small pressure difference across the glial boundary is sufficient to induce significant deformations, hence (7e) reduces to P0≈Pe and U0≈0. As a result, in the limit Ee→0, (10) becomes an instantaneous equation that cannot admit any through flow. For these reasons, we limited our parametric studies to Young moduli related to the most interesting CSF dynamics that can admit steady streaming, i.e. Ee∈[0.01,1].

All the simulations are carried out for tfin=100 with Δt=10-4 and the results are interpreted in terms of brain tissue deformation D0 and through-flow velocity U0-∂TD0. Since the solution is time-dependent, the corresponding time averages ⟨D0⟩ and ⟨U0-∂TD0⟩ are analyzed, averaging over t∈[50,100] in order to get rid of the initial transient effects. The boundary conditions in pressure are 13 Pa=10-3,Pe=0,Pb=0. 

It is remarkable that, in our model framework, the pressure distribution P0(Z,T), the deformation of the brain tissue D0(Z,T) and the through-flow velocity Ue(Z,T)=U0|R=R1+1+D0-∂TD0 can be derived from each other taking into account the permeability parameter Me and the elasticity parameter Ee 14 Ue=MeP0,D0=Ee-1P0, hence, analyzing the results in terms of one among P0, Ue or D0 provides information about all three these quantities. A direct implication of it is the pressure boundary conditions play the role of boundary constraints for Ue and D0, too. Hence, given an elasticity parameter Ee, regardless of Me, H¯, L and R1, the brain tissue deformation at Z=0 will always be D0|Z=0=Ee-1Pa and at Z=L
D0|Z=L=0. With the same argument, fixing Me and regardless of Ee, H¯, L and R1, the through-flow velocity on the left will always be Ue|Z=0=MePa and on the right Ue|Z=L=0.

Vanishing peristaltic wave
We first consider all the cases with H¯=0, since the flow reaches a steady state (∂TD0=∂TP0=0) and it can be well understood taking into account the exact solution reported for the flow in an annular pipe with a permeable wall. This represents an asymptotic limit of our problem and can therefore be used as a validation case for our solver. Taking the limit of H¯→0, Ee→∞ and assuming Pe, Pa and Pb constant in time, the problem admits a steady solution and since D0≡0, the pressure P0 becomes an instantaneous field (i.e. ∂TP0≡0), reducing (10) to 15 ∂2P0∂Z2+BP0=BPe, with A1≡0 and the constant B given in the section Methods. Considering that B is always negative, i.e. |B|=-B, the solution of (15) is of the form 16 P0=Pe+γ1eZ|B|+γ2e-Z|B|, where γ1 and γ2 are constants determined by applying the boundary conditions at Z=0 (P0|Z=0=Pa) and Z=L (P0|Z=L=Pb)  17a γ1=(Pa-Pe)e-L|B|-Pb-Pee-L|B|-eL|B|,γ2=(Pb-Pe)-Pa-PeeL|B|e-L|B|-eL|B|. 

For validation purpose, the numerical solution of (10) for Ee=105, Pe=2, Pa=5, Pb=0, Me=1, H¯=0, L=5 and R1=5 at t=1 is compared with the exact solution (16), valid only for Ee→∞ and H¯=0. The very good agreement is depicted in Fig. 2.Figure 2 Comparison between the exact steady solution of (16) (solid line) valid for Ee→∞ and H¯=0 and the corresponding numerical solution of (10) (circles) evaluated at t=1 and computed for H¯=0 and Ee=105. The other parameters of the comparison are Pe=2, Pa=5, Pb=0, Me=1, L=5 and R1=5.



A further confirmation of the derivation of our model is provided in the limit of large inner radius R1. For R1→∞, the curvature effect becomes negligible and (15) tends to the equation for the incompressible flow in a plane shallow channel with a permeable wall 18 ∂2P0∂X2-12MeP0=-12MePe, where the non-dimensional plane coordinates are X→=(X,Y). Equation (18) is derived in the section Methods and it implies that limR1→∞1/B0=-12. The correct asymptotic limit of our model is retrieved, as shown in Fig. 3.Figure 3 Asymptotic limit of 1/B0 to match the plane flow coefficient for R1→∞: limR1→∞1/B0=-12.



Fixing Pe=Pb=0 and Pa=10-3 the pressure in the annular pipe reads: 19 P0=Pae(Z-L)|B|-e(L-Z)|B|e-L|B|-eL|B|, hence, the solution is nothing but an exponential relaxation from Pa to Pe=Pb=0. This same trend is observed for all the cases with H¯=0, and they are compared in Fig. 4 for Me=1 and Ee=0.01, 0.1 and 1 at t=50. Since the annular channel flow is a limit for Ee→∞ (i.e. D0→0), upon an increase of stiffness of the brain tissue, the pressure distribution tends to (19). It is remarkable that, for the least rigid brain tissue, i.e. Ee=0.01, the exponential relaxation of P0, D0 and Ue blends soon (i.e. approximately for Z>0.5) with a linear trend which holds in most of the thin film.Figure 4 Pressure distribution in a rigid pipe with a permeable wall (solid line) compared to the pressure in the PVS for H¯=0, Ee=0.01 (circles and dashed-line), 0.1 (squares and dashed-line) and 1 (crosses and dashed-line). In all the cases Me=1, L=2, R1=10 (i.e. B0=-0.079557), Pe=Pb=0 and Pa=10-3 at t=50.



Effect of gap length
The effect of the gap length L is investigated, setting H¯=0.1, R1=10, Me=0.5, Ee=0.1 and varying L∈[2,20]. Four PVS lengths are considered and the corresponding average deformation ⟨D0⟩ is depicted in the top panel of Fig. 5: L=2 (dotted), L=5 (dashed–dotted), L=10 (dashed), L=20 (solid). For all the curves it is noticed that the boundary effect which dominates the average deformation distribution is limited to a couple of wavelengths from the boundaries.

The peak at Z=0 is well understood considering the equivalence P0=Ee-1D0, which therefore fixes a steady Dirichlet boundary condition on D0(Z=0)=EePa, hence ⟨D0⟩(Z=0)=EePa. This value is much larger than the average deformation in the bulk, since the flow in the bulk is strongly influenced by the permeability of the glial boundary (see Fig. 4). The second peak near Z=L is typical of non-transparent boundary conditions for wave propagation problems; the steep negative gradient of ⟨D0⟩(Z→L) is a direct results of the Dirichlet boundary condition D0(Z=L)=EePb=⟨D0⟩(Z=L)=0.

In order to get rid of these boundary effects induced by the simplified pressure boundary conditions, we focus on the bulk area where each curve can be well approximated by a straight line, characterized by the only two coefficients A0 and A1 20 ⟨D0⟩|Z∈[5,L-5]≈A0+A1Z. where the coefficient A0 represents the time- and space-averaged brain tissue deformation and the coefficient A1 is the time-averaged axial rate of change of the brain tissue deformation. It is further noticed that the average through flow ⟨Ue⟩ (derived by ⟨D0⟩ multiplying by MeEe) admits a steady streaming since ∫5L-5⟨Ue⟩dZ≠0. The characterization of the steady streaming via A0 and A1 is one of the main aim of our study, as reported in the followings.

Effect of curvature
The effect of the curvature is discussed, setting H¯=0.1, L=20, Me=0.5, Ee=0.1 and varying R1∈[10,1000]. The bottom panel of Fig. 5 compares the average deformation for the three radii of curvature R1=10 (solid line), R1=100 (dashed line) and R1=1000 (dashed–dotted line). The curvature of the annular PVS has relatively small importance in terms of ⟨D0⟩. Increasing the curvature (↓R1) does not have a monotonic trend on the average deformation in the middle of the liquid film, and it tends to preserve the peak near the inflow and outflow boundaries. As also confirmed by Fig. 3, which plots the asymptotic limit of B0-1, the curvature effect becomes negligible when comparing R1=100 and R1=1000.Figure 5 Top: Effect of PVS length for H¯=0.2, R1=10, Me=1, Ee=0.1 is investigated considering four axial lengths: L=2 (dotted), L=5 (dashed–dotted), L=10 (dashed), L=20 (solid). Bottom: Effect of curvature for H¯=0.2, L=20, Me=0.5, Ee=0.1 investigated considering three inner radii: R1=10 (solid line), R1=100 (dashed line) and R1=1000 (dashed–dotted line).



Highest curvature (R1=10) and longest perivascular gap (L=20)
The following results consider R1=10 and L=20. The effect of the peristaltic wave amplitude H¯∈[0,0.2] and of the brain tissue permeability Me∈[0.1,5] is investigated for three cases: soft (Ee=0.01), medium-stiff (Ee=0.1) and rigid (Ee=1) brain tissue.Figure 6 A0 (left panels) and A1 (right panels) coefficients for the average deformation ⟨D0⟩ for R1=10, L=20, Me∈[0.1,5], and Ee=0.01 (top), Ee=0.1 (middle) and Ee=1 (bottom). Six values of Me are considered: Me=0.1 (∙), 0.2 (▪), 0.5 (♦), 1 (▴), 2 (◀), 5 (▶).



Soft brain tissue
The results of A0 (left-top panel) and A1 (right-top panel) for Ee=0.01 are depicted in Fig. 6. Five values of H¯ are considered, i.e. H¯=0, 0.05, 0.1, 0.15, and 0.2, for each six values of Me: Me=0.1 (∙), 0.2 (▪), 0.5 (♦), 1 (▴), 2 (◀), 5 (▶). The case H¯=0 is the only one admitting a steady state for the flow. This reflects on the time-average deformation ⟨D0⟩ (and on ⟨P0⟩ and ⟨Ue⟩, see (14)), which has an exponential trend matching to a linear profile (see Fig. 4). Moreover, all the cases considered, regardless of Me, Ee and H¯, show a decrease of ⟨D0⟩, ⟨P0⟩ and ⟨Ue⟩ as Z increases, i.e. A1 is always negative. Upon an increase of H¯, also the amplitude of the average brain tissue deformation ⟨D0⟩ increases (A1↑). The coefficient A0 is monotonically decreasing with H¯, even if almost negligibly. For soft brain tissues, A0 and A1 strongly depend on the permeability of the brain tissue Me. For the parameter ranges investigated, A1 and A0 show a monotonic trend decreasing as Me increases, if Ee=0.01.

Employing the equivalence Ue=EeMeD0, the coefficients A0 and A1 are used to plot the fitting approximation of ⟨Ue⟩. For soft brain tissues, Fig. 7 reports the time average of the through-flow profile across the glial boundary for Me∈[0.1,5] and H¯∈[0,0.2]. Each panel of Fig. 7 compares the effect of the peristaltic wave amplitude (black: H¯=0; blue: H¯=0.05; red: H¯=0.1; green: H¯=0.15; cyan: H¯=0.2) keeping constant the permeability parameter Me for Ee=0.01. Figure 7 demonstrates that a steady through flow is due to the peristaltic wave amplitude, which increases the overall through flow, ∫0L⟨Ue⟩dZ, for Me≤1 and decreases it for Me>1.Figure 7 ⟨Ue⟩ for R1=10, L=20, Me∈[0.1,5], and Ee=0.01. Five values of H¯ are considered: H¯=0 (black), 0.05 (blue), 0.1 (red), 0.15 (green), 0.2 (cyan).



The understanding of a steady streaming component in the through flow is an interesting outcome of our model. In fact, considering that H is a zero-mean deformation, the increase of ⟨P0⟩=Ee⟨D0⟩ with H¯ highlights the steady pressure component induced by the traveling wave. This is possible only because the brain tissue is deformable Ee↛∞, hence ∂TD0=Ee-1∂TP0≠0. The presence of a time-derivative in (10) allows a phase shift between D0 and H. To better understand it, let us consider the case of Ee→∞ with H¯≠0. Since the lubrication approximation considers only linear terms of the momentum equation, if Ee→∞, ∂TD0=Ee-1∂TP0=0 and (10) becomes an instantaneous equation. As a consequence, for Ee→∞, the fluid flow becomes fully reversible in time and a symmetric zero-mean deformation H, as the one we consider, would produce a zero-mean streaming ⟨Ue⟩≡0 within a traveling wave period. For Ee↛∞, the time derivative ∂TP0 carries the memory of the previous states and makes the flow non-reversible in time, which allows for steady streaming.

Medium-stiff brain tissue
The results for Ee=0.1 are depicted in Fig. 6: A0 (left-middle panel) and A1 (right-middle panel). The same line-style coding is used to denote different H¯, as for the soft tissue case. The first difference with the soft-tissue case is observed in A0, which is one to two orders of magnitude lower than for soft brain tissues. Once again, this is understood considering the steady case (H¯=0), which reduced to an almost-exponential relaxation profile (see squares in Fig. 4). Hence, the linear profile inherited by soft tissues from H¯=0 vanished for medium-stiff brain tissues reducing A0 of two orders of magnitudes. The increased rigidity Ee further contributes to this reduction of A0 as D0=P0Ee-1. Differently from the soft-tissue case, for Ee=0.1, A0 shows a certain dependence on H¯, which grows monotonically for small permeability parameters Me=0.1 and decreases monotonically when the permeability of the glial boundary is higher. On the other hand, A1 is always negative and independent (up to the accuracy of our numerical simulation) on H¯, and it is remarkably influenced by Me up to becoming almost zero if the permeability parameter is high enough (Me≳1). This is well understood considering that a higher Me implies a faster relaxation of the average pressure to a constant value, as indicated by the coefficient B of (19). Since P0=UeMe-1=D0Ee, this same consideration applies to ⟨D0⟩ and ⟨Ue⟩.Figure 8 ⟨Ue⟩ for R1=10, L=20, Me∈[0.1,5], and Ee=0.1. Five values of H¯ are considered: H¯=0 (black), 0.05 (blue), 0.1 (red), 0.15 (green), 0.2 (cyan).



The hallmark of the steady exponential relaxation due to the pressure gradient is hardly visible when comparing the time-dependent profiles of D0 for H¯=0 and H¯=0.4. This is the direct consequence of the stiffness parameter, since increasing Ee reduces the deformation at Z=0 for a given Pa, i.e. D0|Z=0=Ee-1Pa. The flow is then dominated by the peristaltic wave deformation H which gives rise to an interesting phenomenon: increasing the permeability parameter Me, for very permeable brain tissues Me≳1, the average brain tissue deformation ⟨D0⟩ becomes negative.

As a result, using (14), a negative average deformation ⟨D0⟩<0 implies a suction from the brain to the perivascular space ⟨Ue⟩<0. Hence, increasing Me gives rise to an opposite direction of the steady streaming, which now flows from the brain to the PVS. Based on Fig. 6, the sign change occurs at Me≈0.5 for H¯≥0.15 and at Me≈1
∀H¯. Figure 8 reports the time average of the through-flow profile for Ee=0.1, Me∈[0.1,5] and H¯∈[0,0.2]. The same color coding of Fig. 7 is used. For medium-stiff brain tissues, an increase of the peristaltic wave frequency increases ⟨Ue⟩ if Me=0.1 and decreases ⟨Ue⟩ if Me≥0.2, consistently with the trend of A0 for Ee=0.1.

Rigid brain tissue
The results for Ee=1 are depicted in Fig. 6: A0 (left-bottom panel) and A1 (right-bottom panel) using the same line-style coding of the previous cases. Very similar qualitative considerations done for the medium-stiff brain tissue about ⟨D0⟩ and ⟨Ue⟩ apply to the rigid brain tissue. Upon an increase of Ee, the amplitude of the average deformation ⟨D0⟩ decreases (as expected, see absolute values of A0 and A1). Indeed, we remark that ⟨D0⟩ must become steady and converge to zero if the rigidity of the brain tissue goes to infinite, i.e. limEe→∞⟨D0⟩=limEe→∞D0=0. It is furthermore remarkable that, for Ee=1, the rigidity of the brain tissue further contributes to creating negative deformation regions resulting in ⟨D0⟩ always negative. This has corresponding implications on ⟨Ue⟩, which admits more and more extended suction regions, making permeable stiff brain tissues streaming fluid, in average, exclusively from the glial boundary to the perivascular space. This is clearly demonstrated by Fig. 9, where the average through flow across the glial boundary is depicted using the same template of Figs. 7 and 8.Figure 9 ⟨Ue⟩ for R1=10, L=20, Me∈[0.1,5], and Ee=1. Five values of H¯ are considered: H¯=0 (black), 0.05 (blue), 0.1 (red), 0.15 (green), 0.2 (cyan).



The trend reported in Fig. 6 for A0 has an interesting minimum at about Me≈0.5. Indeed, the integral balance A0≈∫0LD0dZ becomes smaller and smaller, in absolute value, upon an increase of Me, if Me≥0.5. This behavior is understood considering the competition between two opposite effects: (a) increasing the permeability of the glial boundary, more fluid can go through the brain tissue |⟨Ue⟩|↑ (see Fig. 9), hence increasing |⟨D0⟩|, and (b) increasing the permeability parameter Me the brain tissue will oppose less and less resistance to be penetrated, hence |⟨D0⟩|↓. The first effect is dominant for Me≤0.5, and the absolute value of A0 increases with Me; the second effect is more important for Me≥1. As a result, the deformation of the brain tissue reduces in amplitude more and more, if the permeability parameter Me≥1, up to asymptotically leading to an undeformed brain tissues, i.e. limMe→∞⟨D0⟩=limMe→∞D0=0. Since A1 is always four orders of magnitude smaller than A0, the linear component of ⟨P0⟩, ⟨D0⟩ and ⟨Ue⟩ can safely be neglected for rigid brain tissues.

Discussion
The cerebrospinal fluid flow across the glial boundary of the brain tissue has been investigated by means of a tribological model derived from first principles. We demonstrate that the phase shift between the arterial peristaltic wave and the glial boundary deformation is a necessary conditions to break the flow symmetry and have a steady streaming. Depending on the elasticity and permeability parameters of the glial boundary, Ee and Me, the steady streaming either enters or exits the brain. For physiologically relevant parameters, we proved that such flow is almost insensitive to curvature effects of the annular perivascular gap for R1>10, and of the perivascular length if L>5. A very comprehensive characterization of the through flow across the glial boundary is provided within our model framework, quantifying the leading order pressure ⟨P0⟩, deformation ⟨D0⟩ and through flow ⟨Ue⟩ across the glial boundary, averaged in time. A reduced order model can be readily derived for such quantities from our model, implementing the fitting functions ⟨P0⟩≈EeA0+A1Z, ⟨D0⟩≈A0+A1Z and ⟨Ue⟩≈EeMeA0+A1Z for whatever perivascular space with R1>10 and L>5.

Among the major outcomes of our study, we estimate the average leak flow velocity for a large physiologically relevant parameter space, finding that ⟨Ue⟩ ranges between -0.0027≤⟨Ue⟩≤0.0005. Considering that typical peristaltic wave frequencies are ω≈5  Hz, the dimensional average through flow is between -0.0135b  s-1≤⟨ωbUe⟩≤0.0025b  s-1, where 10-6≤b≤1.5×10-4  m is the thickness of the perivascular space. Hence, our model estimates that -2.25  μm/s ≤⟨ωbUe⟩≤0.4  μm/s. We remark that this result is consistent with experimental measurements and other model results, since ⟨ωbUe⟩ is typically some orders of magnitude smaller than maxtωbUe, which is supposed to be in the range of 1 μm/s ≤maxtωbUe≤100 μm/s. In particular, considering CSF transport in the perivascular space, Faghih and Sharp15 also mention that arterial pulsations can account for the physiological flow rates through these high flow-resistant spaces. Overall, our model elucidates the dependence of CSF transport on the factors listed in Table 1, and therefore provides a framework to better understand the effect of physiological parameters on perivascular transport. For example, the model can be used to predict how pathologies known to modify parameters like extracellular matrix stiffness (e.g. glial scarring following central nervous system injury) alter the magnitude and direction of CSF flow. Therefore, in addition to calculating specific flow rates, the model described here improves our conceptual understanding of perivascular transport in the brain.

A few concluding remarks about the model robustness and its possible extensions. Owing to the very small values of the non-dimensional film thickness, i.e. 0.000001<ε<0.0375 (see Table 1), the thin-film approximation represents the most insightful and numerically robust leading-order model for Newtonian creeping flows with a permeable boundary. If we consider the complete axisymmetric creeping flow model, the pressure would depend on both coordinates, R and Z. Still, as ε≪1, the pressure would be a very weak function of R, and passing from the thin-film to the complete creeping flow model would mean a significant increment of the model complexity with negligible advantages at leading-order. On the other hand, assuming that the P does not depend on R, as the simplification (6a) does, does not lead to remarkable model inaccuracies. On top of it, owing to the small ε, solving numerically the creeping flow equations is a much more challenging task than solving the thin-film equations because the creeping flow system becomes stiffer and stiffer the smaller ε is. In a recent paper, Ladron-de-Guevara et al.31 point out that a correct modeling of the outflow boundary condition is important when one wants to model perivascular pumping. We further stress that our model does not include any restrictive assumption on the kind of boundary conditions that can be considered. In fact, the extension of the model to pulsatile boundary conditions is straightforwardly achieved by replacing (7a) and (7b) by Z=0:P0=Pa(t) and Z=L:P0=Pb(t). We further remark that including the pulsatile nature of the boundary conditions can induce an improvement of the model accuracy, and we propose it as a very relevant objective for future studies.

Methods
Analytic details
Plugging the asymptotic expansion (5) into (3), it yields  21a ε3Re∂U0∂T+U0∂U0∂R+W0∂U0∂Z+ε∂U1∂T+U1∂U0∂R+U0∂U1∂R+W1∂U0∂Z+W0∂U1∂Z+O(ε2)=-∂P0∂R-ε∂P1∂R+O(ϵ2),  21b εRe∂W0∂T+U0∂W0∂R+W0∂W0∂Z+ε∂W1∂T+U1∂W0∂R+U0∂W1∂R+W1∂W0∂Z+W0∂W1∂Z+O(ε2)=-∂P0∂Z-ε∂P1∂Z+1R∂W0∂R+εR∂W1∂R+∂2W0∂R2+ε∂2W1∂R2+O(ε2),  21c 1R∂(RU0)∂R+ε∂(RU1)∂R+∂W0∂Z+ε∂W1∂Z+O(ε2)=0. 

Expanding the boundary conditions (4) leads  22a Z=0:P0+εP1+O(ε2)=Pa,  22b Z=L:P0+εP1+O(ε2)=Pb,  22c R=R1+H:U0+εU1+O(ε2)=∂TH,  22d W0+εW1+O(ε2)=0,  22e R=R1+1+D:n→=(nz,nr)=-ε∂ZD0+O(ε2),11+ε2∂ZD02+O(ε3),t→=(tz,tr)=1,ε∂ZD0+O(ε2)1+ε2∂ZD02+O(ε3),U0+εU1-ε∂ZD0W0+O(ε2)1+ε2∂ZD02+O(ε3)=MeP0+εP1-Pe+O(ε2)+∂TP0+εP1-Pe+O(ε2)/Ee,W0+εW1+εU0∂ZD0+O(ε2)=0,D0+εD1+O(ε2)=P0+εP1+O(ε2)-Pe/Ee, 

If Re=O(1) or smaller, the leading order continuity and Navier–Stokes equation read  23a ∂P0∂R=0,  23b ∂P0∂Z=1R∂W0∂R+∂2W0∂R2,  23c 1R∂(RU0)∂R+∂W0∂Z=0.  The system (23) is completed by the boundary conditions at leading order  24a Z=0:P0=Pa,  24b Z=L:P0=Pb,  24c R=R1+H:U0=∂TH,W0=0,  24d R=R1+1+D:n→=(nz,nr)=(0,1),t→=(tz,tr)=(1,0),  24e U0=MeP0-Pe+1Ee∂P0-Pe∂T,  24f W0=0,  24g D0=P0-PeEe. 

Equation (23)b can be recast in the form 25 R∂P0∂Z=∂∂RR∂W0∂R, keeping in mind that P0 is just function of Z and T, and integrating in R, it yields 26 R22∂P0∂Z+C1=R∂W0∂R, where C1 is just function of Z and T. Dividing by R and integrating once again in radial direction, it yields 27 R24∂P0∂Z+C1ln(R)+C2=W0, which corresponds to (8), where C2 is just function of Z and T. The leading-order boundary conditions in W0 are W0|R=R1+H=0 and W0|R=R1+1+D0=0. Substituting them in (27) yields  28a W0|R=R1+H=R1+H24∂P0∂Z+C1lnR1+H+C2=0,  28b W0|R=R1+1+D0=R1+1+D024∂P0∂Z+C1lnR1+1+D0+C2=0.  Subtracting the two equations, we eliminate C2, and determine C1 29 C1=α∂P0∂Z,α=(R1+H)2-(R1+1+D0)24ln(R1+1+D0)-ln(R1+H). By substitution of C1 in (28a), C2 is determined 30 C2=β∂P0∂Z,β=-(R1+H)24-αln(R1+H), where α and β are functions of Z and T.

Equation (23)c can be recast in the form 31 ∂RU0∂R=-R∂W0∂Z. Substituting (27) into (31), it reads 32 ∂RU0∂R=-R34∂2P0∂Z2-Rln(R)∂C1∂Z-R∂C2∂Z, and integrating yields 33 RU0=C3-R416∂2P0∂Z2-R22ln(R)-14∂C1∂Z-R22∂C2∂Z, where C3 is a function of Z and T. Dividing by R, (9) is retrieved 34 U0=C3R-R316∂2P0∂Z2-R2ln(R)-14∂C1∂Z-R2∂C2∂Z. Applying the leading-order boundary conditions on U0, yields  35a U0|R=R1+H=C3R1+H-(R1+H)316∂2P0∂Z2-(R1+H)2ln(R1+H)-14∂C1∂Z-(R1+H)2∂C2∂Z=∂H∂T,  35b U0|R=R1+1+D0=C3R1+1+D0-(R1+1+D0)316∂2P0∂Z2-(R1+1+D0)2ln(R1+1+D0)-14∂C1∂Z-(R1+1+D0)2∂C2∂Z=MeP0-Pe+∂D0∂T=MeP0-Pe+1Ee∂P0-Pe∂T.  Eliminating C3 by combining (35a) and (35b) leads to 36 R1+H∂H∂T+(R1+H)316∂2P0∂Z2+(R1+H)2ln(R1+H)-14∂C1∂Z+(R1+H)2∂C2∂Z=R1+1+D0MeP0-Pe+1Ee∂P0-Pe∂T+(R1+1+D0)316∂2P0∂Z2+(R1+1+D0)2∂C2∂Z+(R1+1+D0)2ln(R1+1+D0)-14∂C1∂Z, which is equivalent to (10). The coefficient C3 is then computed by substituting the solution P0 and its derivatives in (35a). The coefficients in (10) are  37a A0=Eeα4R1+1+D02lnR1+1+D0-1-(R1+H)2R1+1+D02ln(R1+H)-1+Eeβ2R1+1+D0-(R1+H)2R1+1+D0+Ee16R1+1+D03-(R1+H)4R1+1+D0,  37b A1=Ee4∂α∂ZR1+1+D02lnR1+1+D0-1-(R1+H)2R1+1+D02ln(R1+H)-1+Ee2∂β∂ZR1+1+D0-(R1+H)2R1+1+D0,  37c A2=EeMe,  37d α=(R1+H)2-(R1+1+D0)24ln(R1+1+D0)-ln(R1+H),  37e β=-(R1+H)24-αln(R1+H),  37f C1=α∂P0∂Z,  37g C2=β∂P0∂Z,  37h C3=(R1+H)∂H∂T+(R1+H)316∂2P0∂Z2+R1+H42ln(R1+H)-1∂C1∂Z+R1+H2∂C2∂Z. 

The coefficient B of (15) is defined by  38a B0=A0|H=0,D0=0Ee=α|H=0,D0=04R1+12lnR1+1-1-R12R1+12ln(R1)-1+β|H=0,D0=02R1+1-R12R1+1+116R1+13-R14R1+1,  38b B2=A2|H=0,D0=0Ee=Me,  38c α|H=0,D0=0=R12-(R1+1)24ln(R1+1)-ln(R1),  38d β|H=0,D0=0=-R124-αln(R1),  38e B=B2B0. 

Shallow channel with a permeable wall
If Ee→∞, H¯=0 and R1→∞, the flow in a shallow channel with a permeable wall represents an asymptotic limit of our thin film problem. Denoting the channel height with b and the channel length with L, if L≫b, and using the scaling (2), the non-dimensional channel flow problem at leading order reads  39a ∂P0∂Y=0,  39b ∂P0∂X=∂2U0∂Y2,  39c ∂U0∂X+∂V0∂Y=0.  where X→=(X,Y) denotes the non-dimensional plane coordinates, P0 and U→0=(U0,V0) are the pressure and the velocity field at leading order. The system of Eqs. (39) is completed by the boundary conditions  40a V0|Y=0=0,  40b V0|Y=1=Me(P0-Pe),  40c U0|Y=0=0,  40d U0|Y=1=0,  40e P0|X=0=Pa,  40f P0|X=L=Pb.  Considering that P0 is only function of X, and integrating the X-momentum twice in Y, it yields 41 U0=∂P0∂XY22+C~1Y+C~2. Applying the boundary conditions in Y direction, we find C~1=-1/2∂XP0 and C~2=0. Plugging ∂XU0 in the continuity equation and integrating once in Y-direction, it yields 42 V0=-12∂2P0∂X2Y33-Y22+C~3. Applying the boundary conditions in Y direction, we find C~3=0 and the following relation for P0 holds 43 ∂2P0∂X2-12MeP0=-12MePe. 

Convergence test
The Navier–Stokes and continuity equation of an incompressible flow in a perivascular thin film have been reduced to the solution of an equation in the form 44 ∂f∂t+α∂f∂s2+β∂f∂s+γf=σ, where α=α(s,t), β=β(s,t), γ=γ(s,t) and σ=σ(s,t) are known functions, s∈[0,Λ] is the space variable (Z in our thin-film flow) and t∈[0,tfin] denotes the time variable.

We discretize (44) in space making use of a spectral collocation method which employs Gauss–Lobatto nodes based on Chebyshev polynomials. Denoting [DN]∈RN×N and [DN2]∈RN×N the first- and second-order discrete derivation matrices in space constructed using N Chebyshev–Gauss–Lobatto nodes, (44) discretized in s reads 45 ∂f→N∂t+α→N[DN2]f→N+β→N[DN]f→N+γ→Nf→N=σ→N, where f→N, α→N, β→N, γ→N and σ→N are N×1 arrays which gather the values of f, α, β, γ and σ at the location of the N nodes at each instant of time t. The time discretization is carried out using the implicit Euler scheme. Denoting with tn the current time and with tn+1 the next instant such that Δt=tn+1-tn, the time-discrete version of (45) reads 46 [IN]/Δt+diagα→Nn+1[DN2]+diagβ→Nn+1[DN]+diagγ→Nn+1f→Nn+1=σ→Nn+1+f→Nn/Δt, where the superscripts n and n+1 denote the times tn and tn+1, respectively, [IN] is the N×N identity matrix and diag∗ is the diagonal matrix resulting from distributing the N×1 array ∗ along the diagonal of an N×N matrix.Figure 10 Maximum in time of the infinite norm of the error function (||Err||∞, bullets), slope-1 line assumed as reference for the solver accuracy (dashed line). The two insets depict the infinite norm of the numerical error ||Err||∞ as function of t for Δt=0.1 and 0.0005.



To test the numerical implementation of our code, we assume  47a α=-sin2π(s-t)-1.05,  47b β=sinπ(s-2t)/6,  47c γ=cos(2s),  47d σ=πcosπ(s+t)-απ2sinπ(s+t)+βπcosπ(s+t)+γsinπ(s+t),  such that the exact solution of (44) is 48 f(s,t)=sinπ(s+t). Dirichlet boundary conditions are derived from (48) and set at s=0 and s=Λ=4 in (46), together with the initial condition f→N0=f(s→N,t=0)=sinπs→N. We stress that the arbitrary choices made in (47) are representative of the problem of interest in our study.

The numerical solution f→N is then compared to the exact solution at each time step by means of the infinite norm of the error function Err(tn)=f(s→N,t=tn)-f→Nn computed at each time point. The simulations are carried out for t=tfin=100 setting N=100 and varying Δt. Figure 10 depicts the convergence curve of the error function, which demonstrate the correctness of our numerical code. The bullets denote the maximum in time of ||Err||∞, depicting it in a log-log plot against the Δt to demonstrate that the solver is first-order accurate in time (see dashed line with slope 1), as expected. The infinite norm of the numerical error ||Err||∞ is plotted as function of time for the largest and the smallest time step (i.e. Δt=0.1 and 0.0005, respectively) in the two insets of fig. 10.

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements
This work was initiated during JBG’s visit to Auckland supported by the University of Auckland Vice Chancellor’s Distinguished Visitor Award.

Author contributions
F.R. derived the current lubrication model based on a previous version proposed by V.S. and J.B.G. All the authors contributed to the design of the research, starting from a preliminary study of V.S. and J.B.G. F.R. programmed the numerical solver, carried out all simulations, and analyzed the data. All the authors edited the article.

Code availability
The code used in this paper is an in-house developed software that will be made available upon request to the corresponding author.

Competing interests
The authors declare no competing interests.
==== Refs
References
1. Iliff, J. J., Wang, M., Liao, Y., Plogg, B. A., Peng, W., Gundersen, G. A., Benveniste, H., Vates, G. E., Deane, R., Goldman, S. A. & Nagelhus, E. A. A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid β. Sci. Trans. Med.4, 147–111 (2012).
2. Iliff JJ  Wang M  Zeppenfeld DM  Venkataraman A  Plog BA  Liao Y  Deane R  Nedergaard M   Cerebral arterial pulsation drives paravascular CSF-interstitial fluid exchange in the murine brain J. Neurosci. 2013 33 18190 18199 10.1523/JNEUROSCI.1592-13.2013 24227727 
3. Hadaczek, P. et al. The, “perivascular pump” driven by arterial pulsation is a powerful mechanism for the distribution of therapeutic molecules within the brain. Mol. Ther.14, 69–78 (2006).
4. Mestre H  Tithof J  Du T  Song W  Peng W  Sweeney AM  Olveda G  Thomas JH  Nedergaard M  Kelley DH   Flow of cerebrospinal fluid is driven by arterial pulsations and is reduced in hypertension Nat. Commun. 2018 9 4878 10.1038/s41467-018-07318-3 30451853 
5. Saleem, A., Qaiser, A., Nadeem, S., Ghalambaz, M. & Issakhov, A. Physiological flow of non-Newtonian fluid with variable density inside a ciliated symmetric channel having compliant wall. Arab. J. Sci. Eng.1–12, (2020).
6. Saleem A  Akhtar S  Alharbi FM  Nadeem S  Ghalambaz M  Issakhov A   Physical aspects of peristaltic flow of hybrid nano fluid inside a curved tube having ciliated wall Res. Phys. 2020 19 103431 
7. Saleem A  Akhtar S  Nadeem S  Alharbi FM  Ghalambaz M  Issakhov A   Mathematical computations for peristaltic flow of heated non-Newtonian fluid inside a sinusoidal elliptic duct Phys. Scr. 2020 95 105009 10.1088/1402-4896/abbaa3 
8. Wang P  Olbricht WL   Fluid mechanics in the perivascular space J. Theor. Biol. 2011 274 52 57 10.1016/j.jtbi.2011.01.014 21241713 
9. Schley D  Carare-Nnadi R  Please CP  Perry VH  Weller RO   Mechanisms to explain the reverse perivascular transport of solutes out of the brain J. Theor. Biol. 2006 238 962 974 10.1016/j.jtbi.2005.07.005 16112683 
10. Kyrtsos, C. R. & Baras, J. S. Modeling the role of the glymphatic pathway and cerebral blood vessel properties in Alzheimer’s disease pathogenesis. PLoS ONE10, e0139574 (2015).
11. Ratner V  Gao Y  Lee H  Nedergaard M  Benveniste H  Tannenbaum A   Cerebrospinal fluid and interstitial fluid motion via the glymphatic pathway modelled by optimal mass transport NeuroImage 2017 152 530 537 10.1016/j.neuroimage.2017.03.021 28323163 
12. Asgari M  de Zélicourt D  Kurtcuoglu V   Glymphatic solute transport does not require bulk flow Sci. Rep. 2016 6 38635 10.1038/srep38635 27929105 
13. 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 17 10.3389/fnagi.2019.00001 30740048 
14. Jin, B.-J., Smith, A. J. & Verkman, A. S. Spatial model of convective solute transport in brain extracellular space does not support a “glymphatic” mechanism. J. Gen. Physiol.148, 489–501 (2016).
15. Faghih MM  Keith Sharp M   Is bulk flow plausible in perivascular, paravascular and paravenous channels? Fluids Barriers CNS 2018 15 17 10.1186/s12987-018-0103-8 29903035 
16. Rey J  Sarntinoranont M   Pulsatile flow drivers in brain parenchyma and perivascular spaces: a resistance network model study Fluids Barriers CNS 2018 15 20 10.1186/s12987-018-0105-6 30012159 
17. Albargothy NJ  Johnston DA  MacGregor-Sharp M  Weller RO  Verma A  Hawkes CA  Carare RO   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 152 10.1007/s00401-018-1862-7 29754206 
18. Xie X  Zhang X  Fu J  Wang H  Jonas JB  Peng X  Tian G  Xian J  Ritch R  Li L  Kang Z  Zhang S  Yang D  Wang N   Noninvasive intracranial pressure estimation by orbital subarachnoid space measurement: the Beijing Intracranial and Intraocular Pressure (iCOP) study Critical Care 2013 17 R162 10.1186/cc12841 23883736 
19. Sakka L  Coll G  Chazal J   Anatomy and physiology of cerebrospinal fluid Eur. Ann. Otorhinolary. 2011 128 309 316 10.1016/j.anorl.2011.03.002 
20. Lightfoot EN   Transport Phenomena and Living Systems 1974 New York Wiley 
21. Ichimura T  Fraser PA  Cserr HF   Distribution of extracellular tracers in perivascular spaces of the rat brain Brain Res. 1991 545 103 113 10.1016/0006-8993(91)91275-6 1713524 
22. Bilston LE  Fletcher DF  Brodbelt AR  Stoodley MA   Arterial pulsation-driven cerebrospinal fluid flow in the perivascular space: a computational model Comput. Method. Biomech. Biomed. Eng. 2003 6 235 241 10.1080/10255840310001606116 
23. Gladdish S  Manawadu D  Banya W  Cameron J  Bulpitt CJ  Rajkumar C   Repeatability of non-invasive measurement of intracerebral pulse wave velocity using transcranial Doppler Clin. Sci. 2005 108 433 439 10.1042/CS20040251 
24. Eng J  Hoeks APG  Brands PJ  Willigers JM  Reneman RS   Non-invasive measurement of mechanical properties of arteries in health and disease Proc. Inst. Mech. Eng. H 1999 213 195 202 10.1243/0954411991534924 10490292 
25. Bloomfield IG  Johnston IH  Bilston LE   Effects of proteins, blood cells and glucose on the viscosity of cerebrospinal fluid Pediatr. Neurosurg. 1998 28 246 251 10.1159/000028659 9732257 
26. Thorin-Trescases N  Bartolotta T  Hyman N  Penar PL  Walters CL  Bevan RD  Bevan JA   Diameter dependence of myogenic tone of human pial arteries: possible relation to distensibility Stroke 1997 28 2486 2492 10.1161/01.STR.28.12.2486 9412638 
27. Baumbach GL  Heistad DD  Siems JE   Effect of sympathetic nerves on composition and distensibility of cerebral arterioles in rats J. Physiol. 1989 416 123 140 10.1113/jphysiol.1989.sp017753 2607446 
28. Atabek HB   Wave propagation through a viscous fluid contained in a tethered, initially stressed, orthotropic elastic tube Biophys. J. 1968 8 626 649 10.1016/S0006-3495(68)86512-9 5699800 
29. Helton ES  Palladino S  Ubogu EE   A novel method for measuring hydraulic conductivity at the human blood-nerve barrier in vitro Microvasc. Res. 2017 109 1 6 10.1016/j.mvr.2016.08.005 27592219 
30. Li G  Yuan W  Fu BM   A model for the blood-brain barrier permeability to water and small solutes J. Biomech. 2010 43 2133 2140 10.1016/j.jbiomech.2010.03.047 20434157 
31. Ladron-de-Guevara, A., Shang, J. K., Nedergaard, M. & Kelley, D. H. Perivascular pumping in the mouse brain: realistic boundary conditions reconcile theory, simulation, and experiment, bioRxiv (2020).