Peristaltic flow in the glymphatic system

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.

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.
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 interstitium 1,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 studies [2][3][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 effects [5][6][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 Olbricht 8 studied axial flow and transport in an annulus with impermeable boundaries, but did not address fluid exchange with the interstitium 9 . 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 stiffness 10 . 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 Sharp 15 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 Sarntinoranont 16 , 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.
The thickness of the PVS is b, the average radius of the artery and the average inner radius of the brain tissue are r 1 and r 2 , respectively. The deformation about the average radii are h(z, t) =h sin [2π/ (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 p e , whereas p a and p b < p a are assumed as pressures at the extrema of the PVS.
Governing equations. The Navier-Stokes and continuity equations in dimensional form read www.nature.com/scientificreports/ 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 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 where Re = ερωb 2 /µ is the Reynolds number.

The mathematical problem (3) is closed by the boundary conditions
where M e = k g µ/bε 2 , P e is the non-dimensional pressure at the glial boundary and k g the permeability of the brain tissue, E e = E g ε 2 /µω(R 1 + 1) and E g the Young's modulus of the brain tissue.
is the imposed deformation of the peristaltic wave of the artery and H =h/b its amplitude measured from the middle line R 1 = r 1 /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 R 2 = r 2 /b = R 1 + 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 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 www.nature.com/scientificreports/

The system (23) is completed by the boundary conditions at leading order
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. W 0 = 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 P 0 is only function of Z and T and, integrating in r the z-momentum, W 0 reads where C 1 = C 1 (Z, T) and C 2 = C 2 (Z, T).
Plugging (8) in the continuity equation and integrating in R, we derive the form of U 0 Applying the no-slip boundary conditions at the inner boundary ( R = R 1 + H ), i.e. U 0 = ∂ T H and W 0 = 0 , and the permeable boundary condition at the outer boundary ( R = R 1 + 1 + D ), i.e. U 0 = M e (P 0 − P e ) + E −1 e ∂ T (P 0 − P e ) and W 0 = 0 , the functions C 1 , C 2 and C 3 are expressed in terms of P 0 and a second-order differential equation is derived for P 0 We refer to the section Methods for the definition of A 0 , A 1 , A 2 . The initial condition for P 0 is set to be the linear function consistent with the boundary conditions 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 D n 0 is employed when computing the solution at time t n+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.

Physiological parameters.
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 × 10 3 Pa, in accordance with Sakka et al. 19 , who report p e ∈ [1300, 2000] Pa. Assuming P e = P b = 0 as reference pressure in our model, P a will then be considered in the range P a ∈ [5 × 10 −7 , 5 ×  26 Moreover, considering that the artery wall thickness t a is 10 to 100 times smaller than r 1 , i.e. t a ∈ [10 −6 , 10 −5 ] m 27 , that the elastic modulus of the artery wall is E a ∈ [10 5 , 10 6 ] N m −2 and that the wall density is ρ a ≈ 10 3 kg m −3 , the elastic wave speed is c a = √ E a t a /2r 1 ρ a ∈ [0.2, 7] m s −1 . Considering the relationship given by Atabek 28 , i.e. c ∈ [0.1, 0.5]c a , an estimated range of the wave speed of the peristaltic wave can be proposed, which includes the estimate of Wang and Olbrich 8 : 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 E g ≈ 0.01E a ∈ [10 3 , 10 4 ] Pa and k g ≈ 10 −11 m Pa −1 s −1 , respectively 29,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 We remark that the range of the glial boundary elasticity parameter E e 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, E e ≪ 1 , a small pressure difference across the glial boundary is sufficient to induce significant deformations, hence (7e) reduces to P 0 ≈ P e and U 0 ≈ 0 . As a result, in the limit E e → 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. E e ∈ [0.01, 1].
All the simulations are carried out for t fin = 100 with t = 10 −4 and the results are interpreted in terms of brain tissue deformation D 0 and through-flow velocity U 0 − ∂ T D 0 . Since the solution is time-dependent, the corresponding time averages D 0 and �U 0 − ∂ T D 0 � are analyzed, averaging over t ∈ [50, 100] in order to get rid of the initial transient effects. The boundary conditions in pressure are It is remarkable that, in our model framework, the pressure distribution P 0 (Z, T) , the deformation of the brain tissue D 0 (Z, T) and the through-flow velocity U e (Z, T) = U 0 | R=R 1 +1+D 0 − ∂ T D 0 can be derived from each other taking into account the permeability parameter M e and the elasticity parameter E e hence, analyzing the results in terms of one among P 0 , U e or D 0 provides information about all three these quantities. A direct implication of it is the pressure boundary conditions play the role of boundary constraints for U e and D 0 , too. Hence, given an elasticity parameter E e , regardless of M e , H , L and R 1 , the brain tissue deformation at Z = 0 will always be D 0 | Z=0 = E −1 e P a and at Z = L D 0 | Z=L = 0 . With the same argument, fixing M e and regardless of E e , H , L and R 1 , the through-flow velocity on the left will always be U e | Z=0 = M e P a and on the right U e | Z=L = 0. www.nature.com/scientificreports/ Vanishing peristaltic wave. We first consider all the cases with H = 0 , since the flow reaches a steady state ( ∂ T D 0 = ∂ T P 0 = 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 , E e → ∞ and assuming P e , P a and P b constant in time, the problem admits a steady solution and since D 0 ≡ 0 , the pressure P 0 becomes an instantaneous field (i.e. ∂ T P 0 ≡ 0 ), reducing (10) to with A 1 ≡ 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 where γ 1 and γ 2 are constants determined by applying the boundary conditions at Z = 0 ( P 0 | Z=0 = P a ) and For validation purpose, the numerical solution of (10) for E e = 10 5 , P e = 2 , P a = 5 , P b = 0 , M e = 1 , H = 0 , L = 5 and R 1 = 5 at t = 1 is compared with the exact solution (16), valid only for E e → ∞ and H = 0 . The very good agreement is depicted in Fig. 2.
A further confirmation of the derivation of our model is provided in the limit of large inner radius R 1 . For R 1 → ∞ , the curvature effect becomes negligible and (15) tends to the equation for the incompressible flow in a plane shallow channel with a permeable wall where the non-dimensional plane coordinates are � X = (X, Y ) . Equation (18) is derived in the section Methods and it implies that lim R 1 →∞ 1/B 0 = −12 . The correct asymptotic limit of our model is retrieved, as shown in Fig. 3.
Fixing P e = P b = 0 and P a = 10 −3 the pressure in the annular pipe reads: hence, the solution is nothing but an exponential relaxation from P a to P e = P b = 0 . This same trend is observed for all the cases with H = 0 , and they are compared in Fig. 4 for M e = 1 and E e = 0.01 , 0.1 and 1 at t = 50 . Since the annular channel flow is a limit for E e → ∞ (i.e. D 0 → 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. E e = 0.01 , the exponential relaxation of P 0 , D 0 and U e blends soon (i.e. approximately for Z > 0.5 ) with a linear trend which holds in most of the thin film.    www.nature.com/scientificreports/ The peak at Z = 0 is well understood considering the equivalence P 0 = E −1 e D 0 , which therefore fixes a steady Dirichlet boundary condition on D 0 (Z = 0) = E e P a , hence �D 0 �(Z = 0) = E e P a . 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 �D 0 �(Z → L) is a direct results of the Dirichlet boundary condition D 0 (Z = L) = E e P b = �D 0 �(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 A 0 and A 1 where the coefficient A 0 represents the time-and space-averaged brain tissue deformation and the coefficient A 1 is the time-averaged axial rate of change of the brain tissue deformation. It is further noticed that the average through flow U e (derived by D 0 multiplying by M e E e ) admits a steady streaming since  (14)), which has an exponential trend matching to a linear profile (see Fig. 4). Moreover, all the cases considered, regardless of M e , E e and H , show a decrease of D 0 , P 0 and U e as Z increases, i.e. A 1 is always negative. Upon an increase of H , also the amplitude of the   Figure 7 demonstrates that a steady through flow is due to the peristaltic wave amplitude, which increases the overall through flow, L 0 U e dZ , for M e ≤ 1 and decreases it for M e > 1. 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 �P 0 � = E e �D 0 � with H highlights the steady pressure component induced by the traveling wave. This is possible only because the brain tissue is deformable E e → ∞ , hence ∂ T D 0 = E −1 e ∂ T P 0 � = 0 . The presence of a time-derivative in (10) allows a phase shift between D 0 and H. To better understand it, let us consider the case of E e → ∞ with H � = 0 . Since the lubrication approximation considers only linear terms of the momentum equation, if E e → ∞ , ∂ T D 0 = E −1 e ∂ T P 0 = 0 and (10) becomes an instantaneous equation. As a consequence, for E e → ∞ , 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 �U e � ≡ 0 within a traveling wave period. For E e → ∞ , the time derivative ∂ T P 0 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 E e = 0.1 are depicted in Fig. 6: A 0 (left-middle panel) and A 1 (rightmiddle 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 A 0 , 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 almostexponential 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 A 0 of two orders of magnitudes. The increased rigidity E e further contributes to this reduction of A 0 as D 0 = P 0 E −1 e . Differently from the soft-tissue case, for E e = 0.1 , A 0 shows a certain dependence on H , which grows monotonically for small permeability parameters M e = 0.1 and decreases monotonically when the permeability of the glial boundary is higher. On the other hand, A 1 is always negative and independent (up to the accuracy of our numerical simulation) on H , and it is remarkably influenced by M e up to becoming almost zero if the permeability parameter is high enough ( M e 1 ). This is well understood considering that a higher M e implies a faster relaxation of the average pressure to a constant value, as indicated by the coefficient B of (19). Since P 0 = U e M −1 e = D 0 E e , this same consideration applies to D 0 and U e .
The hallmark of the steady exponential relaxation due to the pressure gradient is hardly visible when comparing the time-dependent profiles of D 0 for H = 0 and H = 0.4 . This is the direct consequence of the stiffness parameter, since increasing E e reduces the deformation at Z = 0 for a given P a , i.e. D 0 | Z=0 = E −1 e P a . The flow is then dominated by the peristaltic wave deformation H which gives rise to an interesting phenomenon: increasing the permeability parameter M e , for very permeable brain tissues M e 1 , the average brain tissue deformation D 0 becomes negative.
As a result, using (14), a negative average deformation D 0 < 0 implies a suction from the brain to the perivascular space U e < 0 . Hence, increasing M e gives rise to an opposite direction of the steady streaming, Rigid brain tissue. The results for E e = 1 are depicted in Fig. 6: A 0 (left-bottom panel) and A 1 (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 D 0 and U e apply to the rigid brain tissue. Upon an increase of E e , the amplitude of the average deformation D 0 decreases (as expected, see absolute values of A 0 and A 1 ). Indeed, we remark that D 0 must become steady and converge to zero if the rigidity of the brain tissue goes to infinite, i.e. lim E e →∞ �D 0 � = lim E e →∞ D 0 = 0 . It is furthermore remarkable that, for E e = 1 , the rigidity of the brain tissue further contributes to creating negative deformation regions resulting in D 0 always negative. This has corresponding implications on U e , 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.  www.nature.com/scientificreports/ The trend reported in Fig. 6 for A 0 has an interesting minimum at about M e ≈ 0.5 . Indeed, the integral balance A 0 ≈ L 0 D 0 dZ becomes smaller and smaller, in absolute value, upon an increase of M e , if M e ≥ 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 |�U e �| ↑ (see Fig. 9), hence increasing |�D 0 �| , and (b) increasing the permeability parameter M e the brain tissue will oppose less and less resistance to be penetrated, hence |�D 0 �| ↓ . The first effect is dominant for M e ≤ 0.5 , and the absolute value of A 0 increases with M e ; the second effect is more important for M e ≥ 1 . As a result, the deformation of the brain tissue reduces in amplitude more and more, if the permeability parameter M e ≥ 1 , up to asymptotically leading to an undeformed brain tissues, i.e. lim M e →∞ �D 0 � = lim M e →∞ D 0 = 0 . Since A 1 is always four orders of magnitude smaller than A 0 , the linear component of P 0 , D 0 and U e 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, E e and M e , 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 R 1 > 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 P 0 , deformation D 0 and through flow U e 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 �P 0 � ≈ E e (A 0 + A 1 Z) , �D 0 � ≈ A 0 + A 1 Z and �U e � ≈ E e M e (A 0 + A 1 Z) for whatever perivascular space with R 1 > 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 U e ranges between −0.0027 ≤ �U e � ≤ 0.0005 . Considering that typical peristaltic wave frequencies are ω ≈ 5 Hz, the dimensional average through flow is between −0.0135b s −1 ≤ �ωbU e � ≤ 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 ≤ �ωbU e � ≤ 0.4 µm/s. We remark that this result is consistent with experimental measurements and other model results, since ωbU e is typically some orders of magnitude smaller than max t ωbU e , which is supposed to be in the range of 1 µm/s ≤ max t ωbU e ≤ 100 µm/s. In particular, considering CSF transport in the perivascular space, Faghih and Sharp 15 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 thinfilm 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 : P 0 = P a (t) and Z = L : P 0 = P b (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)  www.nature.com/scientificreports/ keeping in mind that P 0 is just function of Z and T, and integrating in R, it yields where C 1 is just function of Z and T. Dividing by R and integrating once again in radial direction, it yields which corresponds to (8), where C 2 is just function of Z and T. The leading-order boundary conditions in W 0 are W 0 | R=R 1 +H = 0 and W 0 | R=R 1 +1+D 0 = 0 . Substituting them in (27) yields Subtracting the two equations, we eliminate C 2 , and determine C 1 By substitution of C 1 in (28a), C 2 is determined where α and β are functions of Z and T. Equation (23) .
(30) 3 16 www.nature.com/scientificreports/ which is equivalent to (10). The coefficient C 3 is then computed by substituting the solution P 0 and its derivatives in (35a). The coefficients in (10) are The coefficient B of (15) is defined by Shallow channel with a permeable wall. If E e → ∞ , H = 0 and R 1 → ∞ , 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 , www.nature.com/scientificreports/ where � X = (X, Y ) denotes the non-dimensional plane coordinates, P 0 and � U 0 = (U 0 , V 0 ) are the pressure and the velocity field at leading order. The system of Eqs. (39) is completed by the boundary conditions Considering that P 0 is only function of X, and integrating the X-momentum twice in Y, it yields Applying the boundary conditions in Y direction, we find C 1 = −1/2 ∂ X P 0 and C 2 = 0 . Plugging ∂ X U 0 in the continuity equation and integrating once in Y-direction, it yields Applying the boundary conditions in Y direction, we find C 3 = 0 and the following relation for P 0 holds 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 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, t fin ] 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 [D N ] ∈ R N×N and [D 2 N ] ∈ R N×N the first-and second-order discrete derivation matrices in space constructed using N Chebyshev-Gauss-Lobatto nodes, (44) discretized in s reads 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 t n the current time and with t n+1 the next instant such that t = t n+1 − t n , the time-discrete version of (45) reads where the superscripts n and n + 1 denote the times t n and t n+1 , respectively, [I N ] 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.
To test the numerical implementation of our code, we assume (43) ∂ 2 P 0 ∂X 2 − 12M e P 0 = −12M e P e .
(44) ∂f ∂t + α ∂f ∂s 2 + β ∂f ∂s + γ f = σ , 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(t n ) = f (� s N , t = t n ) − � f n N computed at each time point. The simulations are carried out for t = t fin = 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.

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.