 Research
 Open Access
 Published:
Hydraulic resistance of periarterial spaces in the brain
Fluids and Barriers of the CNSvolume 16, Article number: 19 (2019)
Abstract
Background
Periarterial spaces (PASs) are annular channels that surround arteries in the brain and contain cerebrospinal fluid (CSF): a flow of CSF in these channels is thought to be an important part of the brain’s system for clearing metabolic wastes. In vivo observations reveal that they are not concentric, circular annuli, however: the outer boundaries are often oblate, and the arteries that form the inner boundaries are often offset from the central axis.
Methods
We model PAS crosssections as circles surrounded by ellipses and vary the radii of the circles, major and minor axes of the ellipses, and twodimensional eccentricities of the circles with respect to the ellipses. For each shape, we solve the governing Navier–Stokes equation to determine the velocity profile for steady laminar flow and then compute the corresponding hydraulic resistance.
Results
We find that the observed shapes of PASs have lower hydraulic resistance than concentric, circular annuli of the same size, and therefore allow faster, more efficient flow of cerebrospinal fluid. We find that the minimum hydraulic resistance (and therefore maximum flow rate) for a given PAS crosssectional area occurs when the ellipse is elongated and intersects the circle, dividing the PAS into two lobes, as is common around pial arteries. We also find that if both the inner and outer boundaries are nearly circular, the minimum hydraulic resistance occurs when the eccentricity is large, as is common around penetrating arteries.
Conclusions
The concentric circular annulus assumed in recent studies is not a good model of the shape of actual PASs observed in vivo, and it greatly overestimates the hydraulic resistance of the PAS. Our parameterization can be used to incorporate more realistic resistances into hydraulic network models of flow of cerebrospinal fluid in the brain. Our results demonstrate that actual shapes observed in vivo are nearly optimal, in the sense of offering the least hydraulic resistance. This optimization may well represent an evolutionary adaptation that maximizes clearance of metabolic waste from the brain.
Background
It has long been thought that flow of cerebrospinal fluid (CSF) in perivascular spaces plays an important role in the clearance of solutes from the brain [1,2,3]. Experiments have shown that tracers injected into the subarachnoid space are transported preferentially into the brain through periarterial spaces at rates much faster than can be explained by diffusion alone [4,5,6]. Recent experimental results from Bedussi et al. [7] and Mestre et al. [8] now show unequivocally that there is pulsatile flow in the perivascular spaces around pial arteries in the mouse brain, with net (bulk) flow in the same direction as the blood flow. The in vivo measurements of Mestre et al. support the hypothesis that this flow is driven primarily by “perivascular pumping” due to motions of the arterial wall synchronized with the cardiac cycle. From the continuity equation (expressing conservation of mass), we know that this net flow must continue in some form through other parts of the system (e.g., along perivascular spaces around penetrating arteries, arterioles, capillaries, venules). This is supported by recent magnetic resonance imaging studies in humans that have demonstrated that CSF tracers are transported deeply into the brain via perivascular spaces [9,10,11].
The in vivo experimental methods of Mestre et al. [8] now enable measurements of the size and shape of the perivascular spaces, the motions of the arterial wall, and the flow velocity field in great detail. With these in vivo measurements, direct simulations can in principle predict the observed fluid flow by solving the Navier–Stokes (momentum) equation. These studies provide important steps in understanding the fluid dynamics of the entire glymphatic system [3, 12], not only in mice but in mammals generally. A handful of numerical [13,14,15,16,17,18] and analytical [19, 20] studies have previously been developed to model CSF flow through PASs. However, these studies have been based on idealized assumptions and have typically simulated fluid transport through only a small portion of the brain. Development of a fullyresolved fluiddynamic model that captures CSF transport through the entire brain is beyond current capabilities for two reasons: (i) the very large computational cost of such a simulation, and (ii) the lack of detailed knowledge of the configuration and mechanical properties of the various flow channels throughout the glymphatic pathway, especially deep within the brain. We note that these limitations and the modest number of publications modeling CSF transport through the brain are in contrast with the much more extensive body of research modeling CSF flow in the spinal canal, which has pursued modeling based on idealized [21,22,23], patientspecific [24, 25], and in vitro [26] geometries (see the recent review articles [27,28,29]).
To simulate CSF transport at a brainwide scale, a tractable first step is to model the flow using a hydraulic network by estimating the hydraulic resistance of the channels that carry the CSF, starting with the PASs. This article is restricted to modeling of CSF flow through PASs in the brain and does not address the question of flow through the brain parenchyma [30, 31], a region where bulk flow phenomena have not been characterized in the same detail as in the PAS. A steady laminar (Poiseuille) flow of fluid down a channel is characterized by a volume flow rate \(\overline{Q}\) that is proportional to the pressure drop \(\Delta p\) along the channel. The inverse of that proportionality constant is the hydraulic resistance \(\overline{\mathcal {R}}\). Higher hydraulic resistance impedes flow, such that fewer mL of CSF are pumped per second by a given pressure drop \(\Delta p\); lower hydraulic resistance promotes flow. Hydraulic resistance is analogous to electrical resistance, which impedes the electrical current driven by a given voltage drop. The hydraulic resistance of a channel for laminar flow can be calculated from the viscosity of the fluid and the length, shape, and crosssectional area of the channel. We note that prior numerical studies have computed the hydraulic resistance of CSF flow in the spinal canal [32, 33], and a few hydraulicnetwork models of periarterial flows have been presented, using a concentric circularannulus configuration of the PAS crosssection (e.g., [16, 34, 35]). As we demonstrate below, the concentric circular annulus is generally not a good model of the crosssection of a PAS. Here we propose a simple but more realistic model that is adjustable and able to approximate the crosssections of PASs actually observed in the brain. We then calculate the velocity profile, volume flow rate, and hydraulic resistance for Poiseuille flow with these crosssections and demonstrate that the shapes of PASs around pial arteries are nearly optimal.
Methods
The basic geometric model of the PAS
In order to estimate the hydraulic resistance of PASs, we need to know the various sizes and shapes of these spaces in vivo. Recent measurements of periarterial flows in the mouse brain by Mestre et al. [8] show that the PAS around the pial arteries is much larger than previously estimated—comparable to the diameter of the artery itself. In vivo experiments using fluorescent dyes show similar results [36]. The size of the PAS is substantially larger than that shown in previous electron microscope measurements of fixed tissue. Mestre et al. demonstrate that the PAS collapses during fixation: they find that the ratio of the crosssectional area of the PAS to that of the artery itself is on average about 1.4 in vivo, whereas after fixation this ratio is only about 0.14.
The in vivo observation of the large size of the PAS around pial arteries is important for hydraulic models because the hydraulic resistance depends strongly on the size of the channel crosssection. For a concentric circular annulus of inner and outer radii \(r_1\) and \(r_2\), respectively, for fixed \(r_1\) the hydraulic resistance scales roughly as \((r_2/r_1)^{4}\), and hence is greatly reduced in a wider annulus. As we demonstrate below, accounting for the actual shapes and eccentricities of the PASs will further reduce the resistance of hydraulic models.
Figure 1 shows images of several different crosssections of arteries and surrounding PASs in the brain, measured in vivo using fluorescent dyes [6, 8, 36, 37] or optical coherence tomography [7]. The PAS around a pial artery generally forms an annular region, elongated in the direction along the skull. For an artery that penetrates into the parenchyma, the PAS is less elongated, assuming a more circular shape, but not necessarily concentric with the artery. Note that similar geometric models have been used to model CSF flow in the cavity (ellipse) around the spinal cord (circle) [21, 22].
We need a simple working model of the configuration of a PAS that is adjustable so that it can be fit to the various shapes that are actually observed, or at least assumed. Here we propose the model shown in Fig. 2. This model consists of an annular channel whose crosssection is bounded by an inner circle, representing the outer wall of the artery, and an outer ellipse, representing the outer wall of the PAS. The radius \(r_1\) of the circular artery and the semimajor axis \(r_2\) (xdirection) and semiminor axis \(r_3\) (ydirection) of the ellipse can be varied to produce different crosssectional shapes of the PAS. With \(r_2 = r_3 > r_1\), we have a circular annulus. Generally, for a pial artery, we have \(r_2 > r_3 \approx r_1\): the PAS is annular but elongated in the direction along the skull. For \(r_3 = r_1 < r_2\), the ellipse is tangent to the circle at the top and bottom, and for \(r_3 \le r_1 < r_2\) the PAS is split into two disconnected regions, one on either side of the artery, a configuration that we often observe for a pial artery in our experiments. We also allow for eccentricity in this model, allowing the circle and ellipse to be nonconcentric, as shown in Fig. 2b. The center of the ellipse is displaced from the center of the circle by distances c and d in the x and y directions, respectively. Using these parameters, we have fit circles and ellipses to the images shown in Fig. 1b–f. Specifically, the fitted circles and ellipses have the same centroids and the same normalized second central moments as the dyed regions in the images. The parameters for the fits are provided in Table 1, and the goodness of these fits can be quantified via the residuals. We define \(A_{out}\) as the image area excluded from the fitted PAS shape even though its color suggests it should be included, and \(A_{in}\) as the image area included in the fitted PAS shape even though its color suggests it should be excluded. Those residuals, normalized by the PAS area, are also listed in Table 1. The model is thus able to match quite well the various observed shapes of PASs. To illustrate the fits, in Fig. 1 we have drawn the inner and outer boundaries (thin and thick white curves, respectively) of the geometric model. We have drawn the full ellipse indicating the outer boundary of the PAS to clearly indicate the fit, but the portion which passes through the artery is plotted with a dotted line to indicate that this does not represent an anatomical structure.
Steady laminar flow in the annular tube
We wish to find the velocity distribution for steady, fully developed, laminar viscous flow in our model tube, driven by a uniform pressure gradient in the axial (z) direction. The velocity u(x, y) is purely in the zdirection and the nonlinear term in the Navier–Stokes equation is identically zero. The basic partial differential equation to be solved is the zcomponent of the Navier–Stokes equation, which reduces to
where \(\mu \) is the dynamic viscosity of the CSF. (Note that the pressure gradient dp/dz is constant and negative, so the constant C we have defined here is positive.) If we introduce the nondimensional variables
then Eq. (1) becomes the nondimensional Poisson’s equation
We want to solve this equation subject to the Dirichlet (noslip) condition \(U=0\) on the inner (circle) and outer (ellipse) boundaries. Analytic solutions are known for simple geometries, and we can calculate numerical solutions for a wide variety of geometries, as described below.
Let \(A_{pas}\) and \(A_{art}\) denote the crosssectional areas of the PAS and the artery, respectively. Now, define the nondimensional parameters
(Note that K is also equal to the volume ratio \(V_{pas}/V_{art}\) of a fixed length of our tube model.) When \(r_1\), \(r_2\), \(r_3\), c, and d have values such that the ellipse surrounds the circle without intersecting it, the crosssectional areas of the PAS and the artery are given simply by
and the area ratio is
In cases where the ellipse intersects the circle, the determination of \(A_{pas}\) is more complicated: in this case, Eqs. (5) and (6) are no longer valid, and instead we compute \(A_{pas}\) numerically, as described in more detail below.
For our computations of velocity profiles in cases with no eccentricity (\(c = d = 0\)), we can choose a value of the area ratio K, which fixes the volume of fluid in the PAS, and then vary \(\alpha \) to change the shape of the ellipse. Thus we generate a twoparameter family of solutions: the value of \(\beta \) is fixed by the values of K and \(\alpha \). In cases where the circle does not protrude past the boundary of the ellipse, the third parameter \(\beta \) varies according to \(\beta = (K+1)/\alpha \). For \(\alpha = 1\) the ellipse and circle are tangent at \(x= \pm r_2\), \(y=0\) and for \(\alpha = K+1\) they are tangent at \(x=0\), \(y = \pm r_3\). Hence, for fixed K, the circle does not protrude beyond the ellipse for \(\alpha \) in the range \(1 \le \alpha \le K+1\). For values of \(\alpha \) outside this range, we have a twolobed PAS, and the relationship among K, \(\alpha \), and \(\beta \) is more complicated.
The dimensional volume flow rate \(\overline{Q}\) is found by integrating the velocityprofile
where \(Q = \overline{Q}/Cr_1^4\) is the dimensionless volume flow rate. The hydraulic resistance \(\overline{\mathcal {R}}\) is given by the relation \(\overline{Q} = \Delta p / \overline{\mathcal {R}}\), where \(\Delta p = (dp/dz) L\) is the pressure drop over a length L of the tube. For our purposes, it is better to define a hydraulic resistance per unit length, \(\mathcal {R} = \overline{\mathcal {R}}/L\), such that
We can use computed values of Q to obtain values of the hydraulic resistance \(\mathcal {R}\). From Eqs. (7) and (8), we have
We can then plot the scaled, dimensionless resistance \(r_1^4 \mathcal {R}/\mu = 1/Q\) as a function of \((\alpha  \beta )/K\) (shape of the ellipse) for different values of K (area ratio). We choose the quantity \((\alpha  \beta )/K\) because it is symmetric with respect to exchange of \(\alpha \) and \(\beta \), larger values of this quantity correspond to a more elongated ellipse, and \((\alpha  \beta )/K=\pm 1\) corresponds to the case in which the ellipse is tangent with the circle.
For viscous flows in ducts of various crosssections, the hydraulic resistance is often scaled using the hydraulic radius \(r_{\text{h}} = 2A/P\), where A is the crosssectional area of the duct and P is the wetted perimeter. In the case of our annular model, however, the hydraulic radius \(r_{\text{h}} = 2A_{pas}/P\) is not a useful quantity: when the inner circle lies entirely within the outer ellipse, both \(A_{pas}\) and P, and hence \(r_{\text{h}}\), are independent of the eccentricity, but (as shown below) the hydraulic resistance varies with eccentricity.
Numerical methods
In order to solve Poisson’s Eq. (3) subject to the Dirichlet condition \(U=0\) on the inner and outer boundaries of the PAS, we employ the Partial Differential Equation (PDE) Toolbox in MATLAB. This PDE solver utilizes finiteelement methods and can solve Poisson’s equation in only a few steps. First, the geometry is constructed by specifying a circle and an ellipse (the ellipse is approximated using a polygon with a high number of vertices, typically 100). Eccentricity may be included by shifting the centers of the circle and ellipse relative to each other. We specify that the equation is to be solved in the PAS domain corresponding to the part of the ellipse that does not overlap with the circle. We next specify the Dirichlet boundary condition \(U=0\) along the boundary of the PAS domain and the coefficients that define the nondimensional Poisson’s Eq. (3). Finally, we generate a fine mesh throughout the PAS domain, with a maximum element size of 0.02 (nondimensionalized by \(r_1\)), and MATLAB computes the solution to Eq. (3) at each mesh point. The volume flow rate is obtained by numerically integrating the velocity profile over the domain. Choosing the maximum element size of 0.02 ensures that the numerical results are converged. Specifically, we compare the numerically obtained value of the flow rate Q for a circular annulus to the analytical values given by Eq. (11) or Eq. (12) below to ensure that the numerical results are accurate to within 1%.
For the case where the circle protrudes beyond the boundary of the ellipse, Eqs. (5) and (6) do not apply. We check for this case numerically by testing whether any points defining the boundary of the circle extend beyond the boundary of the ellipse. If so, we compute the area ratio K numerically by integrating the area of the finite elements in the PAS domain (\(A_{art}\) is known but \(A_{pas}\) is not). In cases where we want to fix K and vary the shape of the ellipse (e.g. Fig. 5a), it is necessary to change the shape of the ellipse iteratively until K converges to the desired value. We do so by choosing \(\alpha \) and varying \(\beta \) until K converges to its desired value within 0.01%.
Analytical solutions
There are two special cases for which there are explicit analytical solutions, and we can use these solutions as checks on the numerical method.
The concentric circular annulus
For a concentric circular annulus we have \(c=d=0\), \(r_2 = r_3 > r_1\), \(\alpha = \beta >1\), and \(K = \alpha ^2 1\). Let r be the radial coordinate, and \(\rho = r/r_1\) be the corresponding dimensionless radial coordinate. The dimensionless velocity profile is axisymmetric, and is given by White [38], p. 114:
and the corresponding dimensionless volume flux rate is given by:
The eccentric circular annulus
There is also an analytical solution for the case of an eccentric circular annulus, in which the centers of the two circles do not coincide [38, 39]. Let c denote the radial distance between the two centers. Then, in cases where the two circles do not intersect, the dimensionless volume flow rate is given by White [38], p. 114:
where \(\epsilon = c/r_1\) is the dimensionless eccentricity and
From this solution, it can be shown that increasing the eccentricity substantially increases the flow rate (see Fig. 310 in [38]). This solution can be used as a check on the computations of the effect of eccentricity in our model PAS in the particular case where the outer boundary is a circle.
Results
The eccentric circular annulus
The eccentric circular annulus is a good model for the PASs around some penetrating arteries (see Fig. 1e, f), so it is useful to show how the volume flow rate and hydraulic resistance vary for this model. This is done in Fig. 3a, where the hydraulic resistance (inverse of the volume flow rate) is plotted as a function of the dimensionless eccentricity \(c/(r_2  r_1) = \epsilon /(\alpha  1)\) for various values of the area ratio \(K = \alpha ^2  1\). The first thing to notice in this plot is how strongly the hydraulic resistance depends on the crosssectional area of the PAS (i.e., on K). For example, in the case of a concentric circular annulus (\(\epsilon = 0\)), the resistance decreases by about a factor of 1700 as the area increases by a factor of 15 (K goes from 0.2 to 3.0).
For fixed K, the hydraulic resistance decreases monotonically with increasing eccentricity (see Fig. 3a). This occurs because the fluid flow concentrates more and more into the wide part of the gap, where it is farther from the walls and thus achieves a higher velocity for a given shear stress (which is fixed by the pressure gradient). (This phenomenon is well known in hydraulics, where needle valves tend to leak badly if the needle is flexible enough to be able to bend to one side of the circular orifice.) The increase of flow rate (decrease of resistance) is well illustrated in Fig. 3c–e, which show numerically computed velocity profiles (as color maps) at three different eccentricities. We refer to the case where the inner circle touches the outer circle (\(\epsilon /(\alpha  1) = 1\)) as the “tangent eccentric circular annulus.”
We have plotted the hydraulic resistance as a function of the area ratio K for the concentric circular annulus and the tangent eccentric circular annulus in Fig. 3b. This plot reveals that across a wide range of area ratios, the tangent eccentric circular annulus (shown in Fig. 3e) has a hydraulic resistance that is approximately 2.5 times lower than the concentric circular annulus (shown in Fig. 3c), for a fixed value of K. Intermediate values of eccentricity (\(0\le \epsilon /(\alpha 1) \le 1\)), where the inner circle does not touch the outer circle (e.g., Fig. 3d) correspond to a reduction in hydraulic resistance that is less than a factor of 2.5. The variation with K of hydraulic resistance of the tangent eccentric annulus fits reasonably well to a power law \(r_1^4 \mathcal {R} / \mu = 8.91 K^{2.78}\) throughout most of the range of observed K values, indicated by the gray shaded region in Fig. 3b.
The concentric elliptical annulus
Now we turn to the results for the elliptical annulus in the case where the ellipse and the inner circle are concentric. Figure 4 shows numerically computed velocity profiles for three different configurations with the same area ratio (\(K=1.4\)): a moderately elongated annulus, the case where the ellipse is tangent to the circle at the top and bottom, and a case with two distinct lobes. A comparison of these three cases with the concentric circular annulus (Fig. 3c) shows quite clearly how the flow is enhanced when the outer ellipse is flattened, leading to spaces on either side of the artery with wide gaps in which much of the fluid is far from the boundaries and the shear is reduced. However, Fig. 4c shows a reduction in the volume flow rate (i.e. less pink in the velocity profile) compared to Fig. 4a, b, showing that elongating the outer ellipse too much makes the gaps narrow again, reducing the volume flow rate (increasing the hydraulic resistance). This results suggests that, for a given value of K (given crosssectional area), there is an optimal value of the elongation \(\alpha \) that maximizes the volume flow rate (minimizes the hydraulic resistance).
To test this hypothesis, we computed the volume flow rate and hydraulic resistance as a function of the shape parameter \((\alpha  \beta )/K\) for several values of the area ratio K. The results are plotted in Fig. 5a. Note that the plot is only shown for \((\alpha  \beta )/K \ge 0\), since the curves are symmetric about \((\alpha  \beta )/K = 0\). The left end of each curve (\((\alpha  \beta )/K = 0\)) corresponds to a circular annulus, and the black circles indicate the value of \(\mathcal {R}\) given by the analytical solution in Eq. (11). These values agree with the corresponding numerical solution to within 1%. The resistance varies smoothly as the outer elliptical boundary becomes more elongated, and our hypothesis is confirmed: for each curve, the hydraulic resistance reaches a minimum value at a value of \((\alpha  \beta )/K\) that varies with K, such that the corresponding shape is optimal for fast, efficient CSF flow. Typically, the resistance drops by at least a factor of two as the outer boundary goes from circular to the tangent ellipse. If we elongate the ellipse even further (beyond the tangent case), thus dividing the PAS into two separate lobes, the resistance continues to decrease but reaches a minimum and then increases. The reason for this increase is that, as the ellipse becomes highly elongated, it forms a narrow gap itself, and the relevant length scale for the shear in velocity is the width of the ellipse, not the distance to the inner circle. For small values of K, we find that the optimal shape parameter \((\alpha  \beta )/K\) tends to be large and the ellipse is highly elongated, while for large values of K the optimal shape parameter is small. The velocity profiles for three optimal configurations (for \(K=0.4\), 1.4, and 2.4) are plotted in Fig. 5c–e.
The hydraulic resistance of shapes with optimal elongation also varies with the area ratio K, as shown in Fig. 5b. As discussed above, the resistance decreases rapidly as K increases and is lower than the resistance of concentric, circular annuli, which are also shown. We find that the optimal elliptical annulus, compared to the concentric circular annulus, provides the greatest reduction in hydraulic resistance for the smallest area ratios K. Although the two curves converge as K grows, they differ substantially throughout most of the range of normalized PAS areas observed in vivo. We find that the variation with K of hydraulic resistance of optimal shapes fits closely to a power law \(r_1^4 \mathcal {R} / \mu = 6.67 K^{1.96}\).
The eccentric elliptical annulus
We have also calculated the hydraulic resistance for cases where the outer boundary is elliptical and the inner and outer boundaries are not concentric (see Fig. 2b). For this purpose, we introduce the nondimensional eccentricities
The hydraulic resistance is plotted in Fig. 6a, b as a function of \(\epsilon _x\) and \(\epsilon _y\), respectively, and clearly demonstrates that adding any eccentricity decreases the hydraulic resistance, similar to the eccentric circular annulus shown in Fig. 3. In the case where the outer boundary is a circle (\(\alpha = \beta > 1\), \(\epsilon = (\epsilon _x^2 + \epsilon _y^2)^{1/2}\)) we employ the analytical solution (12) as a check on the numerical solution: they agree to within 0.4%. Two example velocity profiles are plotted in Fig. 6c, d. Comparing these profiles to the concentric profile plotted in Fig. 4a clearly shows that eccentricity increases the volume flow rate (decreases the hydraulic resistance).
In vivo PASs near pial arteries are nearly optimal in shape
We can compute the velocity profiles for the geometries corresponding to the actual pial PASs shown in Fig. 1b–d (dotted and solid white lines). The parameters corresponding to these fits are provided in Table 1 and are based on the model shown in Fig. 2b, which allows for eccentricity. Figure 7a shows how hydraulic resistance varies with elongation for nonconcentric PASs having the same area ratio K and eccentricities \(\epsilon _x\) and \(\epsilon _y\) as the ones in Fig. 1b–d. The computed values of the hydraulic resistance of the actual observed shapes are plotted as purple triangles. For comparison, velocity profiles for the optimal elongation and the exact fits provided in Table 1 are shown in Fig. 7b–d. Clearly the hydraulic resistances of the shapes observed in vivo are very close to the optimal values, but systematically shifted to slightly more elongated shapes. Even when \((\alpha  \beta )/K\) differs substantially between the observed shapes and the optimal ones, the hydraulic resistance \(\mathcal {R}\), which sets the pumping efficiency and is therefore the biologically important parameter, matches the optimal value quite closely.
Discussion
In order to understand the glymphatic system, and various effects on its operation, it will be very helpful to develop a predictive hydraulic model of CSF flow in the PASs. Such a model must take into account two important recent findings: (i) the PASs, as measured in vivo, are generally much larger than the size determined from postfixation data [7, 8, 36] and hence offer much lower hydraulic resistance; and (ii) (as we demonstrate in this paper) the concentric circular annulus model is not a good geometric representation of an actual PAS, as it overestimates the hydraulic resistance. With these two factors accounted for, we can expect a hydraulicnetwork model to produce results in accordance with the actual bulk flow now observed directly in particle tracking experiments [7, 8].
The relatively simple, adjustable model of a PAS that we present here can be used as a basis for calculating the hydraulic resistance for a wide range of observed PAS shapes, throughout the brain and spinal cord. Our calculations demonstrate that accounting for PAS shape can reduce the hydraulic resistance by a factor as large as 6.45 (see Table 1). We estimate that the pressure gradient required to drive CSF through a murine pial PAS ranges between 0.03 and 0.3 mmHg/cm (this calculation is based on the fit parameters for Fig. 1d, b, respectively, and an average flow speed of 18.7 μm/s [8]). Although CSF pressure gradients have not been measured in PASs, the maximum available pressure to drive such flows arises from arterial pulsations and an upper limit can be estimated based on the arterial pulse pressure, which gives a value on the order of 1 mmHg/cm. We note that our improvements to PAS modeling are also relevant for studies of shearenhanced dispersion of solutes through PASs, a phenomenon that recent numerical works [15, 16, 18] have investigated in the case of an oscillatory, zeromean flow.
We raise the intriguing possibility that the noncircular and eccentric configurations of PASs surrounding pial arteries are an evolutionary adaptation that lowers the hydraulic resistance and permits faster bulk flow of CSF. The in vivo images (e.g., those in Fig. 1b–d) reveal that the crosssection of the PAS around a pial artery is not a concentric circular annulus, but instead is significantly flattened and often consists of two separate lobes positioned symmetrically on each side of the artery. Tracers are mostly moving within these separate tunnels and only to a limited extent passing between them. Our imaging of tens of thousands of microspheres has revealed that crossing is rare, indicating almost total separation between the two tunnels. The arrangement of the two PAS lobes surrounding a pial artery not only reduces the hydraulic resistance but may also enhance the stability of the PAS and prevent collapse of the space during excessive movement of the brain within the skull. Additionally, PASs with wide spaces may facilitate immune response by allowing macrophages to travel through the brain, as suggested by Schain et al. [36]. We note that if CSF flowed through a cylindrical vessel separate from the vasculature (not an annulus), hydraulic resistance would be even lower. However, there are reasons that likely require PASs to be annular and adjacent to the vasculature, including: (i) arterial pulsations drive CSF flow [8], and (ii) astrocyte endfeet, which form the outer boundary of the PAS, regulate molecular transport from both arteries and CSF [40, 41].
The configuration of PASs surrounding penetrating arteries in the cortex and striatum is largely unknown [42]. To our knowledge, all existing models are based on information obtained using measurements from fixed tissue. Our own impression, based on years of in vivo imaging of CSF tracer transport, is that the tracers distribute asymmetrically along the wall of penetrating arteries, suggesting that the PASs here are eccentric. Clearly, we need new in vivo techniques that produce detailed maps of tracer distribution along penetrating arteries. Regional differences may exist, as suggested by the finding that, in the human brain, the striate branches of the middle cerebral artery are surrounded by three layers of fibrous membrane, instead of the two layers that surround cortical penetrating arteries [42]. Accurately characterizing the shapes and sizes of the most distal PASs along the arterial tree is very important, as prior work [35] suggests the hydraulic resistance is largest there. We speculate that the configuration of the PASs at these locations may be optimal as well.
An intriguing possibility for future study is that minor changes in the configuration of PAS spaces may contribute to the sleepwake regulation of the glymphatic system [43]. Also, agedependent changes of the configuration of PASs may increase the resistance to fluid flow, possibly contributing to the increased risk of amyloidbeta accumulation associated with aging [44]. Similarly, reactive remodeling of the PASs in the aftermath of a traumatic brain injury may increase the hydraulic resistance of PASs and thereby increase amyloidbeta accumulation.
There are limitations to the modeling presented here, which can be overcome by straightforward extensions of the calculations we have presented. We have intentionally chosen a relatively simple geometry in order to show clearly the dependence of the hydraulic resistance on the size, shape, and eccentricity of the PAS. However, the fits presented in Fig. 1b–f are imperfect and could be better captured using highorder polygons, which is an easy extension of the numerical method we have employed. Our calculations have been performed assuming that PASs are open channels, which is arguably justified—at least for PASs around pial arteries—by the smooth trajectories observed for 1 μm beads flowing through PASs and the observation that these spaces collapse during the fixation process [8]. However, the implementation of a Darcy–Brinkman model to capture the effect of porosity would simply increase the resistance \(\mathcal {R}\), given a fixed flow rate Q and Darcy number Da, by some multiplicative constant.
The hydraulic resistances we have calculated are for steady laminar flow driven by a constant overall pressure gradient. However, recent quantitative measurements in mice have offered substantial evidence demonstrating that CSF flow in PASs surrounding the middle cerebral artery is pulsatile, driven by peristaltic pumping due to arterial wall motions generated by the heartbeat, with mean (bulk) flow in the same direction as the blood flow [8]. We hypothesize that this “perivascular pumping” occurs mainly in the periarterial spaces around the proximal sections of the main cerebral arteries: at more distal locations the wall motions become increasingly passive, and the flow is driven mainly by the pulsatile pressure gradient generated by the perivascular pumping upstream. Viscous, incompressible duct flows due to oscillating pressure gradients (with either zero or nonzero mean) are well understood: it is a linear problem, and analytical solutions are known for a few simple duct shapes. The nature of the solution depends on the dynamic Reynolds number \(R_d = \omega \ell ^2/\nu \), where \(\omega \) is the angular frequency of the oscillating pressure gradient, \(\nu \) is the kinematic viscosity, and \(\ell \) is the length scale of the duct (e.g., the inner radius of a circular pipe, or the gap width for an annular pipe). (Alternatively, the Womersley number \(W = \sqrt{R_d}\) is often used in biofluid mechanics.) When \(R_d<<1\), as it is in the case of flows in PASs,^{Footnote 1} the velocity profile at any instant of time is very nearly that of a steady laminar flow, and the profile varies in time in phase with the oscillating pressure gradient (see White [38], sec. 34.2). In this case, the average (bulk) volume flow rate will be inversely proportional to exactly the same hydraulic resistance that applies to steady laminar flow. Hence, the hydraulic resistances we have computed here will apply to perivascular spaces throughout the brain, except for proximal sections of main arteries where the perivascular pumping is actually taking place.
In PASs where the perivascular pumping is significant, the picture is somewhat different. Here, the flow is actively driven by traveling wave motions of the arterial wall, or in the context of our model PAS, waves along the inner circular boundary. In the case of an elliptical outer boundary, we expect the flow to be threedimensional, with secondary motions in the azimuthal direction (around the annulus, not down the channel), even if the wave along the inner boundary is axisymmetric. Although we have not yet modeled this flow, we can offer a qualitative description based on an analytical solution for perivascular pumping in the case of concentric circular cylinders [19]. The effectiveness of the pumping scales as \((b/\ell )^2\), where b is the amplitude of the wall wave and \(\ell \) is the width of the gap between the inner and outer boundaries. Although this scaling was derived for an infinite domain, we expect it will also hold for one of finite length. For the case of a concentric circular annulus, the gap width \(\ell \) and hence the pumping effectiveness are axisymmetric, and therefore the resulting flow is also axisymmetric. For an elliptical outer boundary, however, the gap width \(\ell \) varies in the azimuthal direction and so will the pumping effectiveness. Hence, there will be pressure variations in the azimuthal direction that will drive a secondary, oscillatory flow in the azimuthal direction, and as a result the flow will be nonaxisymmetric and the streamlines will wiggle in the azimuthal direction. Increasing the aspect ratio \(r_2/r_3\) of the ellipse for a fixed area ratio will decrease the flow resistance but will also decrease the overall pumping efficiency, not only because more of the fluid is placed farther from the artery wall, but also, in cases where the PAS is split into two lobes, not all of the artery wall is involved in the pumping. Therefore, we expect that there will be an optimal aspect ratio of the outer ellipse that will produce the maximum mean flow rate due to perivascular pumping, and that this optimal ratio will be somewhat different from that which just produces the lowest hydraulic resistance. We speculate that evolutionary adaptation has produced shapes of actual periarterial spaces around proximal sections of main arteries that are nearly optimal in this sense.
Conclusions
Periarterial spaces, which are part of the glymphatic system [6], provide a route for rapid influx of cerebrospinal fluid into the brain and a pathway for the removal of metabolic wastes from the brain. In this study, we have introduced an elliptical annulus model that captures the shape of PASs more accurately than the circular annulus model that has been used in all prior modeling studies. We have demonstrated that for both the circular and elliptical annulus models, nonzero eccentricity (i.e., shifting the inner circular boundary off center) decreases the hydraulic resistance (increases the volume flow rate) for PASs. By adjusting the shape of the elliptical annulus with fixed PAS area and computing the hydraulic resistance, we found that there is an optimal PAS elongation for which the hydraulic resistance is minimized (the volume flow rate is maximized). We find that these optimal shapes closely resemble actual pial PASs observed in vivo, suggesting such shapes may be a result of evolutionary optimization.
The elliptical annulus model introduced here offers an improvement for future hydraulic network models of the glymphatic system, which may help reconcile the discrepancy between the small PAS flow speeds predicted by many models and the relatively large flow speeds recently measured in vivo [7, 8]. Our proposed modeling improvements can be used to obtain simple scaling laws, such as the power laws obtained for the tangent eccentric circular annulus in Fig. 3b or the optimal elliptical annulus in Fig. 5b.
Availability of data and materials
All data generated and analyzed in the course of this study are available from the corresponding author upon reasonable request.
Notes
 1.
For example, for \(\omega = 25.13\ \mathrm{s}^{1}\) (corresponding to a pulse rate of 240 bpm), \(\ell = 20\, \mu \hbox {m}\), and \(\nu = 7.0 \times 10^{7} \mathrm{m}^2\ \mathrm{s}^{1}\), we have \(R_d = 1.4 \times 10^{2}\).
Abbreviations
 CSF:

cerebrospinal fluid
 PAS:

periarterial space
References
 1.
Cserr H, Cooper D, Milhorat T. Flow of cerebral interstitial fluid as indicated by the removal of extracellular markers from rat caudate nucleus. Exp Eye Res. 1977;25:461–73.
 2.
Hladky S, Barrand M. Elimination of substances from the brain parenchyma: efflux via perivascular pathways and via the bloodbrain barrier. Fluids Barriers CNS. 2018;15(1):30.
 3.
Plog B, Nedergaard M. The glymphatic system in central nervous system health and disease: past, present, and future. Annu Rev Pathol. 2018;13:379–94.
 4.
Rennels M, Gregory T, Blaumanis O, Fujimoto K, Grady P. Evidence for a ‘paravascular’ fluid circulation in the mammalian central nervous system, provided by the rapid distribution of tracer protein throughout the brain from the subarachnoid space. Brain Res. 1985;326(1):47–63.
 5.
Ichimura T, Fraser P, Cserr H. Distribution of extracellular tracers in perivascular spaces of the rat brain. Brain Res. 1991;545(1):103–13.
 6.
Iliff J, Wang M, Liao Y, Plogg B, Peng W, Gundersen G, Benveniste H, Vates G, Deane R, Goldman S, Nagelhus E, Nedergaard M. A paravascular pathway facilitates CSF flow through the brain parenchyma and the clearance of interstitial solutes, including amyloid \(\beta \). Sci Transl Med. 2012;4(147):111–47.
 7.
Bedussi B, Almasian M, de Vos J, VanBavel E, Bakker E. Paravascular spaces at the brain surface: low resistance pathways for cerebrospinal fluid flow. J Cerebr Blood F Met. 2018;38(4):719–26.
 8.
Mestre H, Tithof J, Du T, Song W, Peng W, Sweeney A, Olveda G, Thomas J, Nedergaard M, Kelley D. Flow of cerebrospinal fluid is driven by arterial pulsations and is reduced in hypertension. Nat Commun. 2018;9(1):4878.
 9.
Ringstad G, Vatnehol S, Eide P. Glymphatic MRI in idiopathic normal pressure hydrocephalus. Brain. 2017;140(10):2691–705.
 10.
Ringstad G, Valnes L, Dale A, Pripp A, Vatnehol S, Emblem K, Mardal K, Eide P. Brainwide glymphatic enhancement and clearance in humans assessed with MRI. JCI insight. 2018;3(13):e121537.
 11.
DeikeHofmann K, Reuter J, Haase R, Paech D, Gnirs R, Bickelhaupt S, Forsting M, Heußel C, Schlemmer H, Radbruch A. Glymphatic pathway of gadoliniumbased contrast agents through the brain: overlooked and misinterpreted. Invest Radiol. 2019;54(4):229–37.
 12.
Jessen N, Munk A, Lundgaard I, Nedergaard M. The glymphatic system: a beginner’s guide. Neurochem Res. 2015;40(12):2583–99.
 13.
Bilston L, Fletcher D, Brodbelt A, Stoodley M. Arterial pulsationdriven cerebrospinal fluid flow in the perivascular space: a computational model. Comput Method Biomech. 2003;6(4):235–41.
 14.
Schley D, CarareNnadi R, Please C, Perry V, Weller R. Mechanisms to explain the reverse perivascular transport of solutes out of the brain. J Theor Biol. 2006;238(4):962–74.
 15.
Asgari M, Zélicourt D, Kurtcuoglu V. Glymphatic solute transport does not require bulk flow. Sci Rep. 2016;6:38635.
 16.
Rey J, Sarntinoranont M. Pulsatile flow drivers in brain parenchyma and perivascular spaces: a resistance network model study. Fluids Barriers CNS. 2018;15(1):20.
 17.
Aldea R, Weller R, Wilcock D, Carare R, Richardson G. Cerebrovascular smooth muscle cells as the drivers of intramural periarterial drainage of the brain. Front Aging Neurosci. 2019;11:1.
 18.
Sharp M, Carare R, Martin B. Dispersion in porous media in oscillatory flow between flat plates: applications to intrathecal, periarterial and paraarterial solute transport in the central nervous system. Fluids Barriers CNS. 2019;16(1):13.
 19.
Wang P, Olbricht W. Fluid mechanics in the perivascular space. J Theor Biol. 2011;274(1):52–7.
 20.
Sharp M, Diem A, Weller R, Carare R. Peristalsis with oscillating flow resistance: a mechanism for periarterial clearance of amyloid beta from the brain. Ann Biomed Eng. 2016;44(5):1553–65.
 21.
Loth F, Yardimci M, Alperin N. Hydrodynamic modeling of cerebrospinal fluid motion within the spinal cavity. J Biomed Eng. 2001;123:71–9.
 22.
Linge S, Haughton V, Løvgren A, Mardal K, Langtangen H. CSF flow dynamics at the craniovertebral junction studied with an idealized model of the subarachnoid space and computational flow analysis. Am J Neuroradiol. 2010;31(1):185–92.
 23.
Heil M, Bertram C. A poroelastic fluidstructure interaction model of syringomyelia. J Fluid Mech. 2016;809:360–89.
 24.
Rutkowska G, Haughton V, Linge S, Mardal K. Patientspecific 3D simulation of cyclic CSF flow at the craniocervical region. Am J Neuroradiol. 2012;33(9):1756–62.
 25.
Støverud K, Langtangen H, Ringstad G, Eide P, Mardal K. Computational investigation of cerebrospinal fluid dynamics in the posterior cranial fossa and cervical subarachnoid space in patients with Chiari I malformation. PloS ONE. 2016;11(10):0162938.
 26.
Pahlavian S, Bunck A, Thyagaraj S, Giese D, Loth F, Hedderich D, Kröger J, Martin B. Accuracy of 4D flow measurement of cerebrospinal fluid dynamics in the cervical spine: an in vitro verification against numerical simulation. Ann Biomed Eng. 2016;44(11):3202–14.
 27.
Elliott N, Bertram C, Martin B, Brodbelt A. Syringomyelia: a review of the biomechanics. J Fluid Struct. 2013;40:1–24.
 28.
Haughton V, Mardal K. Spinal fluid biomechanics and imaging: an update for neuroradiologists. Am J Neuroradiol. 2014;35(10):1864–9.
 29.
Linninger A, Tangen K, Hsu C, Frim D. Cerebrospinal fluid mechanics and its coupling to cerebrovascular dynamics. Annu Rev Fluid Mech. 2016;48:219–57.
 30.
Smith A, Yao X, Dix J, Jin B, Verkman A. Test of the ‘glymphatic’ hypothesis demonstrates diffusive and aquaporin4independent solute transport in rodent brain parenchyma. eLife. 2017;6:e27679.
 31.
Mestre H, Hablitz L, Xavier A, Feng W, Zou W, Pu T, Monai H, Murlidharan G, Castellanos Rivera R, Simon M, Pike M, Plà V, Du T, Kress B, Wang X, Plog B, Thrane A, Lundgaard I, Abe Y, Yasui M, Thomas J, Xiao M, Hirase H, Asokan A, Iliff J, Nedergaard M. Aquaporin4dependent glymphatic solute transport in the rodent brain. eLife. 2018;7:40070.
 32.
Mardal K, Rutkowska G, Linge S, Haughton V. Estimation of CSF flow resistance in the upper cervical spine. Neuroradiol J. 2013;26(1):106–10.
 33.
Tangen K, Hsu Y, Zhu D, Linninger A. CNS wide simulation of flow resistance and drug transport due to spinal microanatomy. J Biomech. 2015;48(10):2144–54.
 34.
Asgari M, De Zélicourt D, Kurtcuoglu V. How astrocyte networks may contribute to cerebral metabolite clearance. Sci Rep. 2015;5:15024.
 35.
Faghih M, Sharp M. Is bulk flow plausible in perivascular, paravascular and paravenous channels? Fluids Barriers CNS. 2018;15:17.
 36.
Schain A, MeloCarrillo A, Strassman A, Burstein R. Cortical spreading depression closes paravascular space and impairs glymphatic flow: implications for migraine headache. J Neurosci. 2017;37(11):2904–15.
 37.
Achariyar T, Li B, Peng W, Verghese P, Shi Y, McConnell E, Benraiss A, Kasper T, Song W, Takano T, Holtzman D, Nedergaard M, Deane R. Glymphatic distribution of CSFderived apoE into brain is isoform specific and suppressed during sleep deprivation. Mol Neurodegener. 2017;1–20.
 38.
White F. Viscous fluid flow. 3rd ed. New York: McGrawHill; 2006.
 39.
Piercy N, Hooper M, Winny H. LIII. Viscous flow through pipes with cores. Philos Mag. 1933;15(99):647–76.
 40.
Rouach N, Koulakoff A, Abudara V, Willecke K, Giaume C. Astroglial metabolic networks sustain hippocampal synaptic transmission. Science. 2008;322(5907):1551–5.
 41.
Pellerin L, Magistretti P. Glutamate uptake into astrocytes stimulates aerobic glycolysis: a mechanism coupling neuronal activity to glucose utilization. Proc Natl Acad Sci. 1994;91(22):10625–9.
 42.
Zhang E, Inman C, Weller R. Interrelationships of the pia mater and the perivascular (VirchowRobin) spaces in the human cerebrum. J Anat. 1990;170:111.
 43.
Xie L, Kang H, Xu Q, Chen M, Liao Y, Thiyagarajan M, O’Donnell J, Christensen D, Nicholson C, Iliff J, Takano T, Deane R, Nedergaard M. Sleep drives metabolite clearance from the adult brain. Science. 2013;342(6156):373–7.
 44.
Kress B, Iliff J, Xia M, Wang M, Wei H, Zeppenfeld D, Xie L, Kang H, Xu Q, Liew J, Plog B, Ding F, Deane R, Nedergaard M. Impairment of paravascular clearance pathways in the aging brain. Ann Neurol. 2014;76(6):845–61.
Acknowlegements
We thank Dan Xue for assistance with illustrations.
Funding
This work was supported by a grant from the NIH/National Institute of Aging (RF1 AG05757501 to MN, JHT, and DHK) and a grant from the Army Research Office (MURI W911NF1910280 to MN, DHK, and JHT).
Author information
Affiliations
Contributions
JHT developed the theoretical ideas and the geometric model and outlined the calculations. JT and DHK carried out the calculations. HM and MN provided information on actual PAS shapes and flows. JHT, JT, and DHK analyzed the results and wrote the paper. All authors read and approved the final manuscript.
Corresponding author
Correspondence to John H. Thomas.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Perivascular flow
 Cerebrospinal fluid
 Bulk flow
 Hydraulic resistance
 Fluid mechanics
 Glymphatic system