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.

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

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 and , respectively. The deformation about the average radii are and d(z, t), where h is the imposed deformation of the travelling peristaltic wave, , 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 , whereas and are assumed as pressures at the extrema of the PVS.

Governing equations

The Navier–Stokes and continuity equations in dimensional form read where denotes the velocity field in cylindrical coordinates, p is the pressure, 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 withwhere is the frequency of peristaltic wave in the artery, b is the thickness of the perivascular fluid film, and . The non-dimensional continuity and Navier–Stokes equations read where is the Reynolds number.

The mathematical problem (3) is closed by the boundary conditions where , is the non-dimensional pressure at the glial boundary and the permeability of the brain tissue, and the Young’s modulus of the brain tissue. is the imposed deformation of the peristaltic wave of the artery and its amplitude measured from the middle line . is the deformation of the peristaltic layer at the boundary with the brain tissue; it is measured from the middle line . The normal and tangent unit vector to the glial boundary are denoted by and .

Thin film approximation

Assuming that , 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, , 1 refers to the linear correction in , , 2 indicates the quadratic correction in , , and so on.

If or smaller, the leading order continuity and Navier–Stokes equations read 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. 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 is only function of Z and T and, integrating in r the z-momentum, readswhere and .

Plugging (8) in the continuity equation and integrating in R, we derive the form of where .

Applying the no-slip boundary conditions at the inner boundary (), i.e. and , and the permeable boundary condition at the outer boundary (), i.e. and , the functions , and are expressed in terms of and a second-order differential equation is derived for We refer to the section Methods for the definition of , , . The initial condition for is set to be the linear function consistent with the boundary conditionsEquation (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 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 is employed when computing the solution at time .

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

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   Pa, in accordance with Sakka et al.[19], who report   Pa. Assuming as reference pressure in our model, will then be considered in the range . Wang and Olbricht[8] report the values measured in previous studies for: the peristaltic wave frequency   Hz[3], the inner radius of the perivascular layer   m[20], the outer radius of the perivascular layer   m[21] (which result in a perivascular film thickness   m), the peristaltic wave amplitude   m[22], the peristaltic wave speed   m s[20,23,24] (which result in a peristaltic wavelength   m), and the dynamic viscosity   Pa s[25]. Plus, the fluid density is comparable to the one of water   kg m. Other values of the artery radius (  m) are reported in Thorin-Trescases et al.[26] Moreover, considering that the artery wall thickness is 10 to 100 times smaller than , i.e.   m[27], that the elastic modulus of the artery wall is   N m and that the wall density is   kg m, the elastic wave speed is   m s. Considering the relationship given by Atabek[28], i.e. , an estimated range of the wave speed of the peristaltic wave can be proposed, which includes the estimate of Wang and Olbrich[8]:   m s resulting in a range for the peristaltic wavelength   m, which includes the previous estimate   m. Another estimate of b is given by Iliff et al.[1], who reports   m, and by Jin et al.[14], who reports   m. Finally, The elastic modulus and the permeability of the brain tissue are   Pa and   m  Pa  s, 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 and . Higher-order corrections in and inertial effects due to are, therefore, substantially negligible.

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 ερωb^2/μ	[5×10^-12,2×10^-4]
ε	Perivascular film thickness b/λ	[10^-6,3.75×10^-2]
R_1	Inner radius of the perivascular layer r_1/b	[0.7,10^3]
H¯	Amplitude of the peristaltic wave h¯/b	[8×10^-4,0.5]
L	Perivascular length l/λ	[2,20]
M_e	Permeability of the brain tissue k_gμ/bε^2	[10^-2,10]
E_e	Stiffness of the brain tissue E_gε^2/μω(R_1+1)	[2×10^-10,0.2]

Parametric study

Based on the physiologically relevant parameters, we carried out numerical simulations for the following range of non-dimensional groupsWe remark that the range of the glial boundary elasticity parameter 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, , a small pressure difference across the glial boundary is sufficient to induce significant deformations, hence (7e) reduces to and . As a result, in the limit , (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. .

All the simulations are carried out for with and the results are interpreted in terms of brain tissue deformation and through-flow velocity . Since the solution is time-dependent, the corresponding time averages and are analyzed, averaging over 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 , the deformation of the brain tissue and the through-flow velocity can be derived from each other taking into account the permeability parameter and the elasticity parameter hence, analyzing the results in terms of one among , or provides information about all three these quantities. A direct implication of it is the pressure boundary conditions play the role of boundary constraints for and , too. Hence, given an elasticity parameter , regardless of , , L and , the brain tissue deformation at will always be and at . With the same argument, fixing and regardless of , , L and , the through-flow velocity on the left will always be and on the right .

Vanishing peristaltic wave

We first consider all the cases with , since the flow reaches a steady state () 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 , and assuming , and constant in time, the problem admits a steady solution and since , the pressure becomes an instantaneous field (i.e. ), reducing (10) towith and the constant B given in the section Methods. Considering that B is always negative, i.e. , the solution of (15) is of the formwhere and are constants determined by applying the boundary conditions at () and ()

For validation purpose, the numerical solution of (10) for , , , , , , and at is compared with the exact solution (16), valid only for and . The very good agreement is depicted in Fig. 2.

Figure 2

Comparison between the exact steady solution of (16) (solid line) valid for and and the corresponding numerical solution of (10) (circles) evaluated at and computed for and . The other parameters of the comparison are , , , , and .

A further confirmation of the derivation of our model is provided in the limit of large inner radius . For , the curvature effect becomes negligible and (15) tends to the equation for the incompressible flow in a plane shallow channel with a permeable wallwhere the non-dimensional plane coordinates are . Equation (18) is derived in the section Methods and it implies that . The correct asymptotic limit of our model is retrieved, as shown in Fig. 3.

Figure 3

Asymptotic limit of to match the plane flow coefficient for : .

Fixing and the pressure in the annular pipe reads:hence, the solution is nothing but an exponential relaxation from to . This same trend is observed for all the cases with , and they are compared in Fig. 4 for and , 0.1 and 1 at . Since the annular channel flow is a limit for (i.e. ), 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. , the exponential relaxation of , and blends soon (i.e. approximately for ) 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 , (circles and dashed-line), 0.1 (squares and dashed-line) and 1 (crosses and dashed-line). In all the cases , , (i.e. ), and at .

Effect of gap length

The effect of the gap length L is investigated, setting , , , and varying . Four PVS lengths are considered and the corresponding average deformation is depicted in the top panel of Fig. 5: (dotted), (dashed–dotted), (dashed), (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.

Figure 5

Top: Effect of PVS length for , , , is investigated considering four axial lengths: (dotted), (dashed–dotted), (dashed), (solid). Bottom: Effect of curvature for , , , investigated considering three inner radii: (solid line), (dashed line) and (dashed–dotted line).

The peak at is well understood considering the equivalence , which therefore fixes a steady Dirichlet boundary condition on , hence . 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 is typical of non-transparent boundary conditions for wave propagation problems; the steep negative gradient of is a direct results of the Dirichlet boundary condition .

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 and where the coefficient represents the time- and space-averaged brain tissue deformation and the coefficient is the time-averaged axial rate of change of the brain tissue deformation. It is further noticed that the average through flow (derived by multiplying by ) admits a steady streaming since . The characterization of the steady streaming via and 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 , , , and varying . The bottom panel of Fig. 5 compares the average deformation for the three radii of curvature (solid line), (dashed line) and (dashed–dotted line). The curvature of the annular PVS has relatively small importance in terms of . Increasing the curvature () 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 , the curvature effect becomes negligible when comparing and .

Highest curvature () and longest perivascular gap ()

The following results consider and . The effect of the peristaltic wave amplitude and of the brain tissue permeability is investigated for three cases: soft (), medium-stiff () and rigid () brain tissue.

Soft brain tissue

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

Figure 6

(left panels) and (right panels) coefficients for the average deformation for , , , and (top), (middle) and (bottom). Six values of are considered: (), 0.2 (), 0.5 (♦), 1 (), 2 (), 5 ().

Employing the equivalence , the coefficients and are used to plot the fitting approximation of . For soft brain tissues, Fig. 7 reports the time average of the through-flow profile across the glial boundary for and . Each panel of Fig. 7 compares the effect of the peristaltic wave amplitude (black: ; blue: ; red: ; green: ; cyan: ) keeping constant the permeability parameter for . Figure 7 demonstrates that a steady through flow is due to the peristaltic wave amplitude, which increases the overall through flow, , for and decreases it for .

Figure 7

for , , , and . Five values of are considered: (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 with highlights the steady pressure component induced by the traveling wave. This is possible only because the brain tissue is deformable , hence . The presence of a time-derivative in (10) allows a phase shift between and H. To better understand it, let us consider the case of with . Since the lubrication approximation considers only linear terms of the momentum equation, if , and (10) becomes an instantaneous equation. As a consequence, for , 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 within a traveling wave period. For , the time derivative 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 are depicted in Fig. 6: (left-middle panel) and (right-middle panel). The same line-style coding is used to denote different , as for the soft tissue case. The first difference with the soft-tissue case is observed in , which is one to two orders of magnitude lower than for soft brain tissues. Once again, this is understood considering the steady case (), which reduced to an almost-exponential relaxation profile (see squares in Fig. 4). Hence, the linear profile inherited by soft tissues from vanished for medium-stiff brain tissues reducing of two orders of magnitudes. The increased rigidity further contributes to this reduction of as . Differently from the soft-tissue case, for , shows a certain dependence on , which grows monotonically for small permeability parameters and decreases monotonically when the permeability of the glial boundary is higher. On the other hand, is always negative and independent (up to the accuracy of our numerical simulation) on , and it is remarkably influenced by up to becoming almost zero if the permeability parameter is high enough (). This is well understood considering that a higher implies a faster relaxation of the average pressure to a constant value, as indicated by the coefficient B of (19). Since , this same consideration applies to and .

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

As a result, using (14), a negative average deformation implies a suction from the brain to the perivascular space . Hence, increasing 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 for and at . Figure 8 reports the time average of the through-flow profile for , and . The same color coding of Fig. 7 is used. For medium-stiff brain tissues, an increase of the peristaltic wave frequency increases if and decreases if , consistently with the trend of for .

Figure 8

for , , , and . Five values of are considered: (black), 0.05 (blue), 0.1 (red), 0.15 (green), 0.2 (cyan).

Rigid brain tissue

The results for are depicted in Fig. 6: (left-bottom panel) and (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 and apply to the rigid brain tissue. Upon an increase of , the amplitude of the average deformation decreases (as expected, see absolute values of and ). Indeed, we remark that must become steady and converge to zero if the rigidity of the brain tissue goes to infinite, i.e. . It is furthermore remarkable that, for , the rigidity of the brain tissue further contributes to creating negative deformation regions resulting in always negative. This has corresponding implications on , 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

for , , , and . Five values of are considered: (black), 0.05 (blue), 0.1 (red), 0.15 (green), 0.2 (cyan).

The trend reported in Fig. 6 for has an interesting minimum at about . Indeed, the integral balance becomes smaller and smaller, in absolute value, upon an increase of , if . 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 (see Fig. 9), hence increasing , and (b) increasing the permeability parameter the brain tissue will oppose less and less resistance to be penetrated, hence . The first effect is dominant for , and the absolute value of increases with ; the second effect is more important for . As a result, the deformation of the brain tissue reduces in amplitude more and more, if the permeability parameter , up to asymptotically leading to an undeformed brain tissues, i.e. . Since is always four orders of magnitude smaller than , the linear component of , and 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, and , 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 , and of the perivascular length if . A very comprehensive characterization of the through flow across the glial boundary is provided within our model framework, quantifying the leading order pressure , deformation and through flow 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 , and for whatever perivascular space with and .

Among the major outcomes of our study, we estimate the average leak flow velocity for a large physiologically relevant parameter space, finding that ranges between . Considering that typical peristaltic wave frequencies are   Hz, the dimensional average through flow is between   s  s, where   m is the thickness of the perivascular space. Hence, our model estimates that   m/s   m/s. We remark that this result is consistent with experimental measurements and other model results, since is typically some orders of magnitude smaller than , which is supposed to be in the range of 1 m/s  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. (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 , 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 and . 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

Expanding the boundary conditions (4) leads

If or smaller, the leading order continuity and Navier–Stokes equation read The system (23) is completed by the boundary conditions at leading order

Equation (23)b can be recast in the formkeeping in mind that is just function of Z and T, and integrating in R, it yieldswhere is just function of Z and T. Dividing by R and integrating once again in radial direction, it yieldswhich corresponds to (8), where is just function of Z and T. The leading-order boundary conditions in are and . Substituting them in (27) yields Subtracting the two equations, we eliminate , and determine By substitution of in (28a), is determinedwhere and are functions of Z and T.

Equation (23)c can be recast in the formSubstituting (27) into (31), it readsand integrating yieldswhere is a function of Z and T. Dividing by R, (9) is retrievedApplying the leading-order boundary conditions on , yields Eliminating by combining (35a) and (35b) leads towhich is equivalent to (10). The coefficient is then computed by substituting the solution 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 , and , 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 , and using the scaling (2), the non-dimensional channel flow problem at leading order reads where denotes the non-dimensional plane coordinates, and are the pressure and the velocity field at leading order. The system of Eqs. (39) is completed by the boundary conditions Considering that is only function of X, and integrating the X-momentum twice in Y, it yieldsApplying the boundary conditions in Y direction, we find and . Plugging in the continuity equation and integrating once in Y-direction, it yieldsApplying the boundary conditions in Y direction, we find and the following relation for 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 formwhere , , and are known functions, is the space variable (Z in our thin-film flow) and 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 and the first- and second-order discrete derivation matrices in space constructed using N Chebyshev–Gauss–Lobatto nodes, (44) discretized in s readswhere , , , and are 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 the current time and with the next instant such that , the time-discrete version of (45) readswhere the superscripts n and denote the times and , respectively, is the identity matrix and is the diagonal matrix resulting from distributing the array along the diagonal of an matrix.

To test the numerical implementation of our code, we assume such that the exact solution of (44) isDirichlet boundary conditions are derived from (48) and set at and in (46), together with the initial condition . We stress that the arbitrary choices made in (47) are representative of the problem of interest in our study.

The numerical solution is then compared to the exact solution at each time step by means of the infinite norm of the error function computed at each time point. The simulations are carried out for setting and varying . 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 , depicting it in a log-log plot against the 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 is plotted as function of time for the largest and the smallest time step (i.e. and 0.0005, respectively) in the two insets of fig. 10.

Figure 10

Maximum in time of the infinite norm of the error function (, bullets), slope-1 line assumed as reference for the solver accuracy (dashed line). The two insets depict the infinite norm of the numerical error as function of t for and 0.0005.