A network model of glymphatic flow under different experimentally-motivated parametric scenarios

Summary Flow of cerebrospinal fluid (CSF) through perivascular spaces (PVSs) in the brain delivers nutrients, clears metabolic waste, and causes edema formation. Brain-wide imaging cannot resolve PVSs, and high-resolution methods cannot access deep tissue. However, theoretical models provide valuable insight. We model the CSF pathway as a network of hydraulic resistances, using published parameter values. A few parameters (permeability of PVSs and the parenchyma, and dimensions of PVSs and astrocyte endfoot gaps) have wide uncertainties, so we focus on the limits of their ranges by analyzing different parametric scenarios. We identify low-resistance PVSs and high-resistance parenchyma as the only scenario that satisfies three essential criteria: that the flow be driven by a small pressure drop, exhibit good CSF perfusion throughout the cortex, and exhibit a substantial increase in flow during sleep. Our results point to the most important parameters, such as astrocyte endfoot gap dimensions, to be measured in future experiments.


INTRODUCTION
The brain lacks lymph vessels, so scientists have questioned whether a flow of cerebrospinal fluid (CSF) might play a pseudo-lymphatic role in transporting metabolic waste products (Milhorat, 1975). Early speculation was motivated by studies that found that tracers injected into the CSF were transported at rates faster than is possible by diffusion alone (Cserr et al., 1981;Rennels et al., 1985). Now, renewed interest has followed the in vivo observations of Iliff et al. (2012), who reported bulk flow of CSF through perivascular spaces (PVSs; annular channels around brain vasculature) of the murine brain, which aids clearance of amyloid-b, a peptide linked to Alzheimer's disease; they named this clearance pathway the ''glymphatic'' (glial-lymphatic) system. Soon thereafter, Xie et al. (2013) demonstrated that this system is active primarily during sleep. Growing evidence suggests that glymphatic dysfunction may contribute to the progression of dementia (Nedergaard and Goldman, 2020) and worsened outcomes following stroke , brain trauma (Sullan et al., 2018), and many other neurological disorders (Rasmussen et al., 2018).
The glymphatic pathway is hypothesized to consist of an influx of CSF along periarterial spaces which subsequently exchanges with extracellular fluid via bulk flow, facilitated by aquaporin-4 channels on the astrocyte endfeet lining the outer wall of PVSs, followed by an efflux along perivenous spaces and nerve sheaths . Recent studies in humans have confirmed many of the key features of the glymphatic hypothesis (Ringstad et al., 2018;Fultz et al., 2019;Eide et al., 2021). Several experimental methods have been used to probe various parts of the glymphatic system. Two-photon microscopy offers excellent temporal and spatial resolution for in vivo measurements, but typically requires invasive surgery to place a cranial window and is limited to regions near the surface of the brain (Iliff et al., 2012;Schain et al., 2017;Mestre et al., 2018bMestre et al., , 2020. Magnetic resonance imaging (MRI) provides noninvasive brain-wide measurements, but temporal and spatial resolution are orders of magnitude lower, rendering PVSs smaller than the spatial resolution (Ringstad et al., 2018;Fultz et al., 2019;Taoka and Naganawa, 2020). Although ex vivo analysis of brain tissue offers high resolution throughout the brain, recent studies have revealed abnormal CSF flow immediately following cardiac arrest (Ma et al., 2019;Du et al., 2021) and collapse of PVSs during tissue fixation (Mestre et al., 2018b), casting doubt on such measurements. Hence, there remains much uncertainty regarding the precise CSF flow pathway and transport rates, including glymphatic efflux routes. Resolving such details may lead to novel strategies for prevention, diagnosis, and treatment of neurological disorders (Rasmussen et al., 2018).
Numerical modeling offers a powerful tool in which governing equations and physical constraints can fill voids where experimental measurements are not feasible. Indeed, much insight into the glymphatic system has already resulted from such studies (see the review articles by ; Thomas (2019); Benveniste et al. (2019); Martinac and Bilston (2020); Rasmussen et al. (2021)). Here we develop numerical models of CSF flow through a substantial portion of the glymphatic system and use this model to make predictions under different scenarios that account for uncertainties in important geometric and material parameters. Because a fully-resolved fluid-dynamic model is not computationally feasible, our approach employs a hydraulic network model, as in prior work (Asgari et al., 2015;Faghih and Sharp, 2018;Rey and Sarntinoranont, 2018;Vinje et al., 2020). We investigate whether most CSF flows through the parenchyma or PVSs surrounding precapillaries, which we model as parallel pathways. Our attention to precapillary PVSs is motivated by (1) early experimental evidence of tracer transport through capillary PVSs (Rennels et al., 1985), (2) recent characterization of molecular markers suggesting PVSs are continuous from arterioles to capillaries to veins , and (3) recent theoretical arguments that diffusive transport in the parenchyma coupled with advective transport in precapillary PVSs might provide an effective clearance mechanism (Thomas, 2019).
In order to improve on prior idealizations of the glymphatic pathway (Faghih and Sharp, 2018;Wang et al., 2021), we have developed a model of CSF flow in the murine brain based on measurements of the vascular connectivity performed by Blinder et al. (2010Blinder et al. ( , 2013. We use the connectivity between different vessels in this model ( Figures 1A-1C) to separately simulate either blood flow (for characterizing the influence of our (B) Circuit schematic of the pial vasculature (black), which has several penetrating arterioles (red) branching from it. (C) Circuit schematic of a penetrating arteriole (red) which has a total of 11 precapillaries (green) branching from it (only three are shown). When we use a similar model to predict glymphatic CSF flow, we also include an equal number of parenchymal channels (purple). The gray circuit elements in B-C are not shown in A. (D-F) Pressure, volume flow rate, and speed for blood flow; in all three cases, the shaded regions indicate the range of values for a real vascular topology reported by Blinder et al. (2013), while the symbols and error bars indicate the mean and range of values, respectively, computed using the idealized geometry shown in panel A. See also Figure S1. iScience Article idealized geometry) or CSF flow. The model includes flow associated with one of the major arteries branching from the circle of Willis, e.g., the middle cerebral artery (MCA), and thus includes flow in approximately one-fifth of the cortex. MRI studies Stanton et al., 2021) show that CSF enters pial PVSs at the circle of Willis, which is represented by the inlet node in our model, labeled in Fig. 1A-1B. The location of the outlet of our model (labeled ''grounded node'' in Figures 1B-1C) is ambiguous, as glymphatic efflux is not yet well-characterized; this could, for example, correspond to the subarachnoid space (SAS) or meningeal lymphatic vessels. The model geometry for the pial vasculature ( Figure 1B) is based on a branching hexagonal model proposed by Blinder et al. (2010), with nine pial generations amounting to 45 hexagonal units and a total of 324 penetrating arterioles. This latter value approximately matches the number of penetrating arterioles in the vicinity of the MCA, 320, which we obtained by inspecting the pial arterial reconstructions available in the Supplemental Material of Adams et al. (2018). From data reported by Blinder et al. (2013), we determined that, on average, 11 precapillaries branch from each of the penetrating arterioles, which we assumed to be uniformly spaced ( Figure 1C). Our hydraulic network model relates flow to the pressure differences that drive the flow and the hydraulic resistances that oppose the flow (pressure and resistance being analogous to voltage and electrical resistance in circuits). Note that this approach describes the time-averaged (net) volume flow rate and therefore neglects the oscillatory component of CSF flow, which is a reasonable approach since the Womersley number for PVS flow is small (Thomas, 2019). For blood flow (or CSF flow), the resistance through the capillary bed (or capillary PVSs) and venous circulation (or venous PVSs) is modeled using single parallel resistors, shown in gray in Figure 1C, with resistance 2:25310 7 mmHg,min/mL (or one mmHg,min/mL); see STAR Methods for details. Parenchymal flow (implemented only for CSF flow) is modeled using hydraulic resistances based on an analytical expression provided by Holter et al. (2017) (see STAR Methods). A full list of the parameters for the model is given in Table 1.

Characterizing the effects of network geometry idealization via blood flow simulations
In order to investigate whether the idealizations of our vascular model (e.g., hexagonal connectivity, homogeneity of pial artery diameter) significantly alter the distribution of flow, we compared blood flow in our model with blood flow predicted for the realistic network measured by Blinder et al. (2013). The idealized network was adjusted to cover an extent of vascular territory similar to that of the Blinder et al. study by matching the number of penetrating arterioles, resulting in a network with two pial generations (three hexagonal units), in contrast to the network shown in Figure 1A, which consists of nine pial generations or 45 hexagonal units. In Blinder et al. (2013), the authors measured the location and radius of all of the vessels within a section of the cortex, noting the connectivity between the vessels, and assigned a resistance to each segment based on a modified Hagen-Poiseuille law, where r is the vessel radius, L is the vessel length, R is the resistance of that segment of vessel, and m is the dynamic viscosity of water. They then applied a constant pressure difference of 50 mmHg between the arterioles and venules at the surface of the cortex and solved for the flow in each vessel. The resulting ranges of pressures, volumetric flow rates, and velocities for one mouse are indicated by the shaded regions shown in Figures 1D-1F (see Figure S1 for results for two more mice). Based on Equation (1) and with a pressure difference of 50 mmHg between the inlet and outlet, we also predicted pressures, volume flow rates, and velocities for the idealized vascular geometry, which are plotted in Figures 1D-1F with solid symbols; the error bars indicate the range of values. The good agreement between the results for the realistic geometry (Blinder et al., 2013) and for the idealized geometry indicates that the idealization does not substantially alter the salient features of blood flow through the network. The smaller range of values observed for the idealized geometry is a result of the homogeneity of the idealization. These insights suggest that the idealized vascular geometry, which provides a framework for modeling glymphatic flow, is reasonable. Though it does not address the geometry of CSF circulation, we can infer that our results predicting glymphatic flow based on this idealized vascular geometry will likely also exhibit a narrower variation in pressure, volume flow rate, and flow speed than the actual network which has much greater heterogeneity.

Dependence of glymphatic flow on permeability and PVS size
To model CSF flow through the glymphatic network, we enabled parenchymal flow (purple stars in Figures 1A and 1C), modeled three different types of PVSs -pial, penetrating, and precapillary -and assumed homogeneity in the shapes, sizes, and porosity of each of these different PVS types (see STAR Methods for a description of how the hydraulic resistance was computed for each pathway). Several variables needed to model fluid flow through the PVSs and parenchyma are unknown or have substantial uncertainty in their estimates. To overcome this challenge, we performed multiple simulations by bracketing the uncertain quantities (i.e., using the highest and lowest estimates of the uncertain quantities), based on a wide survey of the literature. We emphasize that in most cases, these bounds do not represent strict limits on feasible parameter ranges, but rather correspond to the extrema of values that have been reported in the literature or can be inferred from experimental data.
Bracketed parameters are indicated with an asterisk in Table 1. We considered four scenarios that lead to an overall resistance for the glymphatic network that is either maximal (R max ), minimal (R min ), or intermediate (Intermediate scenario 1, 2; i.e., a combination of one maximal and one minimal parameter set). For all these simulations, we matched the median pial PVS velocity to experimental measurements of 18.7 mm/s (Mestre et al., 2018b;Bedussi et al., 2017;Raghunandan et al., 2021); to obtain this match in flow speed, a different effective pressure drop Dp eff was required for each different scenario. We modeled the pial PVSs as open (i.e., not porous) (Min Rivas et al., 2020) with a realistic, oblate shape (Tithof et al., 2019) and a PVS-to-artery cross-sectional area ratio of G pial = 1:4 (Mestre et al., 2018b). In vivo imaging studies suggest that pial PVSs are demarcated from the SAS (Mestre et al., 2018b;Schain et al., 2017), although some fluid may flow between the two compartments through stomata, which are pores up to a few microns in diameter ; our model does not account for stomata. For penetrating PVSs, we used either an approximate upper bound on the area ratio G pial = 1:4 (i.e., that of pial PVSs) or an approximate lower bound on the area ratio G pial = 0:36 (i.e., the upper bound for the precapillary PVSs, discussed below). We modeled flow through the parenchyma, as well as porous penetrating and precapillary PVSs, using Darcy's law; open (non-porous) penetrating PVSs were modeled as a tangent eccentric annulus (Tithof The four different scenarios we modeled arise from combining either the highest or lowest estimate of (1) the total parenchymal resistance and (2) the penetrating and precapillary PVS permeability, as detailed in Table 2. Minimum/maximum estimates of the total parenchymal resistance were obtained by lumping together the resistance from the gaps between astrocyte endfeet and the extracellular space (ECS; Figures 2A and 2B; see STAR Methods). Note that prior studies (Asgari et al., 2015;Jin et al., 2016;Wang et al., 2021) suggest CSF from penetrating PVSs primarily enters the ECS via gaps between astrocyte endfeet. We modeled flow through the gaps in the endfeet as flow between infinite parallel plates, with a fixed endfoot thickness of 0.45 mm based on measurements of the ''intercellular cleft length'' (Mathiisen et al., 2010) and variable width and cavity fraction for the gaps (Mathiisen et al., 2010;Korogod et al., 2015) (see Table 1 and Figure S2). The upper and lower bounds that we set on the parenchymal permeability k par come from two commonly cited studies (Basser, 1992;Holter et al., 2017); multiple other studies Prabhu et al., 1998;Bobo et al., 1994;Neeves et al., 2006;Smith and Humphrey, 2007) have reported k par values within these bounds. Basser (1992) performed experimental measurements that estimated k par = 4:5310 À 15 m 2 . However, Holter et al. (2017) performed a numerical reconstruction of the neuropil, estimating k par = 1:2310 À 17 m 2 , and speculated that the discrepancy with the earlier findings of Basser and other experimental studies may be because of fluid escaping to high-permeability pathways such as PVSs in those experiments. We therefore used this hypothesis as the basis for our R max scenario, with k par = 1:2310 À 17 m 2 and k PVS = 4:5310 À 15 m 2 . For the R min scenario, we supposed that measurements of k par = 4:5310 À 15 m 2 from Basser (1992)  For each of the four scenarios we considered, we varied G precap ( Figure 2A) from 0.07 to 0.36 (i.e., the precapillary PVS gap ranged from 0.1 to 0.5 mm). These values come from estimates of the size of a basement membrane (Yurchenco, 2011) or the endothelial glycocalyx (Reitsma et al., 2007), which are the respective smallest and largest anatomical structures likely to form the contiguous portion of the PVS network at the precapillary level. In general, the anatomical details of which spaces are contiguous with penetrating PVSs are not well-understood; for a more in-depth discussion of potential PVS routes at the level of microvessels, see Hladky and Barrand (2018). For the R min and Intermediate 2 scenarios (with k PVS = k Basser = 4:5310 À 15 m 2 ), we found that an open precapillary PVS would result in an effective permeability k open < k Basser for G precap < 0:16 ( Figure 2E). By definition, k open provides the upper limit on permeability, and because these two scenarios assume k Basser provides the lower limit of PVS permeability, we exclude G precap < 0:16 (i.e., precapillary PVS gap widths below 0.23 mm) from further analysis in these two scenarios.
The effective pressure drop Dp eff required to drive flow through the glymphatic network is plotted in Figure 2F for all four scenarios. By ''effective'' pressure drop, we mean that we have driven flow through the iScience Article network using a single pressure source ( Figure 1B), but the actual pressure gradients driving glymphatic flow -the source of which is actively debated -may be much more complex. The effective pressure drop may be thought of as Q total =R total , where Q total and R total are the total volume flow rate and resistance for the entire glymphatic network, respectively, even if an external pressure drop of that magnitude does not exist. Potential sources of the pressure gradients that drive the observed flows include arterial pulsations (Mestre et al., 2018b), functional hyperemia (Kedarasetti et al., 2020), and osmotic effects Halnes et al., 2019;Rasmussen et al., 2021). The largest effective pressure drop (43 mmHg) is required for the R max case with G precap = 0:17, whereas the R min case requires a drop of only 0.21 mmHg (which does not vary appreciably with G precap ). CSF pressure gradients have never been measured in mice, but a prior study (Penn and Linninger, 2009) numerically estimated that the largest feasible transmantle pressure difference in humans is approximately 1 mmHg, and hence the R min scenario is possible and Intermediate Scenario 1 is of marginal feasibility, whereas the R max and Intermediate 2 scenarios are very unlikely. The total volumetric flow rate through the entire network ( Figure 2G), which is approximately one-fifth of the full cortical glymphatic network, varies from 0.063 to 0.089 mL/min for all cases considered here; this relatively narrow range of values is a consequence of our requirement that the median flow speed in the pial PVSs match experimental measurements (Mestre et al., 2018b). With negligible variation in Q total for a given scenario, the total hydraulic resistance of the network ( Figure 2H) is linearly proportional to the effective pressure drop, resulting in a similar functional dependence for each scenario.
We next investigated the percentage of flow that passes through the parenchyma versus the precapillary PVSs and the associated flow speed for each case ( Figures 2I-2P). We found that when k par = 1:2310 À 17 m 2 (R max and Intermediate 1 scenarios) there is a comparable fraction of total flow through the parenchyma and precapillary PVSs ( Figures 2I and 2M). However, if k par = 4:5310 À 15 m 2 (R min and Intermediate 2 scenarios), virtually all of the flow passes through the parenchyma with a negligible amount passing through the precapillary PVSs ( Figures 2K and 2O). Consequently, only the former two cases show a substantial dependence on G precap , with the percentage of flow through precapillary PVSs varying from 20% to 35% as G precap is varied from 0.17 to 0.36 for the R max scenario, or from 1.8% to 68% as G precap is varied from 0.07 to 0.36 for Intermediate scenario 1. The average flow speeds are plotted in Figures 2J, 2L, 2N, and 2P, with error bars indicating the full range of the data. The mean values for the flow speed through the parenchyma are quite similar for all four scenarios, with the average speed varying from 0.053 mm/s to 0.065 mm/s for the R max scenario, or from 0.060 mm/s to 0.019 mm/s for Intermediate scenario 1, as G precap is increased. For the R min scenario (Intermediate scenario 2), the mean speed is 0.086 (0.081) mm/s and does not vary appreciably with G precap . We caution that the plotted parenchymal flow speeds are not mean values across the parenchyma; they are computed at the outer wall of the PVS, so they should be interpreted as upper bounds on the parenchymal flow speed, which varies spatially. The mean precapillary flow speeds, in contrast to parenchymal speeds, show substantial variation throughout the four scenarios. The average speed varies from 13 mm/s to 10 mm/s for the R max scenario, or from 2.7 mm/s to 20 mm/s for Intermediate scenario 1, as G precap is increased. For the R min scenario, the mean speed varies from 0.0058 to 0.13 mm/s as G precap is increased, but for Intermediate scenario two the mean speed is 0.021 mm/s and does not vary with G precap . Figures S3 and S4 show how the speed and pressure vary throughout the network for the four scenarios each with maximum or minimum G precap .

Quantifying CSF perfusion of tissue for different scenarios
Numerous studies in both humans and mice have reported that tracers injected into CSF penetrate below the brain's surface over relatively short time scales (Gaberel et al., 2014;Ringstad et al., 2018;Eide et al., 2021;Taoka and Naganawa, 2020). Furthermore, there is growing evidence that CSF flow through the Figure 2. Simulations of CSF flow through the glymphatic network for different scenarios (A) Schematic illustrating the geometry of a penetrating PVS segment below the cortical surface (the same segment depicted in Figure 1C), with flow continuing through precapillary PVSs and/or the parenchyma. (B) Circuit schematic for the geometry shown in A (a greater portion of the network is shown in Figures 1B and 1C). Throughout this article, CSF flows through the precapillary PVSs or parenchyma are consistently plotted with green or purple arrows/symbols, respectively. iScience Article glymphatic pathway is important for the removal of metabolic waste (Koundal et al., 2020;Gu et al., 2020;Eide et al., 2021), including amyloid-b (Iliff et al., 2012;Xie et al., 2013;Roberts et al., 2014;Shokri-Kojori et al., 2018), which is produced throughout the brain. Hence, one may reasonably expect a uniform perfusion of CSF throughout the depth of the cortex to explain observations in tracer experiments and the physiological necessity of adequate waste removal. Consequently, we next computed the volume flow rate through pial PVSs, penetrating PVSs, precapillary PVSs, and the parenchyma for each of eight cases (the four scenarios introduced previously, each with either the maximum or minimum value of G precap ), as shown in the left columns of each scenario in Figure 3. It is immediately clear that when penetrating PVS resistance is minimal (R min and Intermediate 1 scenarios), a significant volume of CSF penetrates into the deep cortex ( Figures 3E, 3G, 3I, and 3K). However, if penetrating PVS resistance is high (R max and Intermediate 2 scenarios), the volume flow rate drops off more rapidly with depth ( Figures 3A, 3C, 3M, and 3O). Figure S5 provides a visualization of how the volume flow rate varies throughout the network.
To characterize the CSF perfusion, we plotted the cumulative flow fraction (i.e., the fraction of the total volume flow rate perfused from the surface of the brain to a given depth) in the right columns for each scenario in Figure 3. The R max scenario has fairly poor CSF perfusion, with 81% (G precap = 0:17) to 84% (G precap = 0:36) of the total CSF exiting each penetrating PVS within 270 mm of the surface. Comparing the flows for small versus large precapillary PVSs ( Figures 3B and 3D), it is clear that more flow reroutes through the PVSs in the latter case, consistent with Figure 2I. Among all cases, Intermediate scenario 2 exhibits the worst CSF perfusion, with 100% of the total CSF perfusing within 270 mm of the surface ( Figures 3N  and 3P); this scenario exhibits negligible dependence on G precap . In contrast, the R min scenario exhibits moderately uniform CSF perfusion, with 63% of the total CSF perfusing within 270 mm of the surface ( Figures 3F and 3H); this scenario also exhibits weak dependence on G precap . By far the best CSF perfusion is observed in Intermediate scenario 1, for which 27 and 28% of the CSF is perfused within 270 mm of the surface for G precap = 0:07 and G precap = 0:36, respectively (Figures 3J and 3L; perfectly uniform CSF perfusion corresponds to 27% at 270 mm). Although the total CSF perfusion remains approximately constant for different precapillary PVS sizes, as G precap is increased a greater fraction of the flow reroutes from the parenchyma to the precapillary PVSs (compare Figures 3I-3J with 3K-3L), consistent with the flow fractions plotted in Figure 2M.
The variations in CSF perfusion through the depth of the cortex for these different scenarios can be understood by comparing the hydraulic resistance of individual segments of the network, as shown in Figures 3Q-3X; several of these values are also provided in Table S1. When R pen is substantially smaller than both R par and R precap (Intermediate scenario 1; Figures 3U-3V), excellent, uniform CSF perfusion occurs ( Figures 3J and 3L). However, a lesser separation in resistance values (R max and R min scenarios; Figures 3Q-3T) leads to less uniformity in the CSF perfusion ( Figures 3B, 3D, 3F, and 3H). When R pen is much greater than R par (Intermediate scenario 2; Figures 3W-3X), virtually all fluid exits through the parenchymal nodes closest to the surface of the brain and CSF perfusion is negligible at deeper nodes ( Figures 3N and 3P). The relative flow through the parenchyma versus precapillary PVSs can also be understood by comparing R precap and R par . For cases where there is substantial CSF perfusion, if the value of R precap and R par are comparable ( Figures 3Q-3R and 3V), then a comparable fraction of fluid will flow through each route ( Figures 3B, 3D, and 3L), with greater flow through the path of lower resistance. Two additional points are notable. The value of R pial is much less than R pen in every scenario, which ensures uniform perfusion of CSF across the pial PVS network (i.e., an approximately equal amount of CSF flows through both a distal penetrating PVSs and a proximal one). Also, the uncertainties in the cavity fraction and gap width of the astrocyte endfeet lead to a huge range in possible values of R AE ( Figure 3Y). In the R max and Intermediate 1 scenarios, the astrocyte endfeet constitute a significant barrier to flow entering the parenchyma

Glymphatic flow during wakefulness versus sleep
We carried out additional calculations with our model aimed at investigating the increase in tracer influx during sleep/anesthesia reported by several studies (Xie et al., 2013;Hablitz et al., 2019Hablitz et al., , 2020. The CSF simulations presented up to this point (Figures 2 and 3) correspond to sleep conditions (or, comparably, conditions under ketamine-xylazine anesthesia). To model the change in flow during wakefulness, we used the Kozeny-Carman equation (see STAR Methods) to estimate that k par decreases by a factor of about 5.5 in wakefulness, compared to sleep. We repeated the simulations of the eight scenarios presented in Figure 3 using a parenchymal permeability that was 5.5 times smaller, but all other parameters (including the imposed pressure drop) were left unchanged for each scenario. We then compared these results to the results from each corresponding simulation under sleep conditions. The total volume flow rates through the entire model network for wake and sleep in each of the eight scenarios are plotted in Figures 4A-4H. The combined flow for awake conditions (open gray diamonds) varies substantially across different scenarios, whereas for sleep (filled gray diamonds) all correspond to Q total of either 0.063 or 0.089 mL/min, consistent with Figure 2G.
We quantified the sleep/wake change in flow by plotting the ratio of volume flow rates Q sleep total =Q awake total , shown in Figures 4I-4P. For sleep compared to wakefulness, in every scenario the volume flow rate decreases for precapillary PVSs and increases for the parenchyma, leading to an overall increase in the combined volume flow rate. This is expected because the increased parenchymal permeability during sleep leads to an overall reduction in the hydraulic resistance of the network, and locally this change will reroute some precapillary PVS flow through the parenchyma. We find that the wake/sleep increase in combined volume flow rate is largest for small precapillary PVSs (Figures 4I, 4K, 4M, and 4O); this combined increase is small for Intermediate scenario 2 (0.8%), R max (13%-20%), and R min (21%), but up to 220% for Intermediate scenario 1. In this latter scenario, however, there is substantial sensitivity to the size of the precapillary PVSs ( Figures 4E-4F and 4MÀ4N), which arises because parenchymal flow dominates the combined transport for the G precap = 0:07 case (and is therefore sensitive to wake/sleep changes in k par ), whereas precapillary PVS flow dominates the combined transport for the G precap = 0:36 case (and is therefore insensitive to wake/ sleep changes in k par ); this observation is consistent with Figure 2M. Of all of our wake/sleep simulations, none exhibits an increase in combined flow greater than 220%.

DISCUSSION
In this study, we have developed a numerical model of a substantial portion of the glymphatic system in the murine brain. This model is based on an idealized vascular geometry inspired by detailed measurements reported by Blinder et al. (2010Blinder et al. ( , 2013, and we characterized the effects of idealizing the vascular geometry by first performing simulations of blood flow (Figure 1). In modeling CSF flow through the glymphatic pathway, we matched median pial CSF velocity to experiments (Mestre et al., 2018b), we realistically modeled pial PVSs as open (non-porous) (Min Rivas et al., 2020) and oblate (Tithof et al., 2019), and we used experimentally measured mean vessel diameters and lengths. To overcome the multiple uncertainties in other parameters, we set reasonable bounds (Table 1) and performed several simulations corresponding to different combinations of the extreme values of the uncertain parameters (Figure 2). It should be noted that these bounds are not strict extrema, but rather correspond to maximum/minimum values of each quantity as reported in various experimental studies. This ''bracketing'' approach included upper and lower bounds on the hydraulic resistance for penetrating PVSs, precapillary PVSs, and the parenchyma (based on a lumped model of astrocyte endfeet and the parenchymal ECS). Our model assumes CSF passes from penetrating PVSs to either precapillary PVSs or through the parenchymal ECS via a paracellular route through gaps between astrocyte endfeet . Ultimately, our goal was to investigate different scenarios to test which parameter regimes are feasible and explain as much experimental data as possible. We focused primarily on quantifying the required pressure drops, flow fraction and speed, cortical CSF perfusion, and sleep/wake changes in volumetric flow rate.
The pressure drops and total volumetric flow rates we computed ( Figures 2F-2G) provide novel insights. The two scenarios with high penetrating PVS resistance R pen (R max and Intermediate scenario 2) require infeasibly large pressure drops between 30 and 43 mmHg. This renders both scenarios very unlikely because such a large pressure drop is even greater than the typical systolic-diastolic variation in blood pressure of iScience Article about 20 mmHg (Mattson, 2001), which is thought to provide an absolute upper bound for the pressure drop driving glymphatic flow (Faghih and Sharp, 2018). The R min scenario, however, requires the lowest pressure drop by definition, which is only 0.21 mmHg. Such a pressure drop is feasible and in line with estimates for the transmantle pressure difference (Penn and Linninger, 2009) (i.e., that between the SAS and lateral ventricles); note however that Penn and Linninger (2009)   iScience Article from 0.36 to 0.07, respectively. Such a pressure drop is perhaps of marginal feasibility, possibly requiring driving mechanism(s) beyond simply a transmantle pressure difference (additional mechanisms are discussed further below).
Because we matched the median pial CSF velocity to experimental measurements (Mestre et al., 2018b), we find Q total = 0:064mL/min for every scenario, except the R min scenario in whichQ total = 0:089 ( Figure 2G). The reason Q total is moderately larger in the R min scenario is because the minimal resistances of penetrating PVSs and parenchyma ( Figures 3S-3T) allow more fluid to exit the network along the parenchymal nodes most proximal to the inlet, which is perhaps discernible in Figure S5B. Our model represents approximately one-fifth of the cortical glymphatic network (e.g., in the vicinity of one MCA), so the total CSF volume flow rate through cortical PVSs would be approximately 0.32 mL/min, much larger than the CSF production rate of the choroid plexus, which has recently been measured to be about 0.1 mL/min for young, healthy, anesthetized mice . Although this measurement involves invasive techniques, Karimy et al. (2015) (who developed the technique used by Liu et al. (2020) in rats) reported that results were consistent with a prior method; still, this measurement may be an underestimate, as the approach excludes CSF production at the 4th ventricle.
Multiple potential explanations exist for the discrepancy between estimates of CSF production and the larger volume flow rate from our model, some of which depend on the details of pial PVSs. The pial PVSs that we have modeled are extensions of the SAS, and prior studies have suggested that not all fluid in pial PVSs continues to penetrating PVSs but rather a portion of the flow continues directly from PVSs of pial arteries to those of veins (Bedussi et al., 2017;Pizzo et al., 2018;Ma et al., 2019). Furthermore, not all CSF from the SAS enters pial PVSs; Lee et al. (2018) delivered a tracer to the cisterna magna in rats and determined that approximately 20% reached the parenchyma, with the rest following CSF efflux routes, including the arachnoid villi, cribriform plate, and cranial and spinal nerves. Hence, it is likely that only a fraction of the total CSF enters pial PVSs, and perhaps not all CSF in pial PVSs continues through penetrating PVSs and into the parenchyma. We note that our model does not include direct flow from the SAS into penetrating PVSs . Our prediction of a volume flow rate larger than CSF production thus suggests that either (1) published in vivo measurements of fluid velocities (Mestre et al., 2018b;Bedussi et al., 2017) are inaccurately large, (2) the pial geometry in our model is too idealized and greatly overestimates the volume flow rate, (3) published measurements of CSF production rates are inaccurately small, and/or (4) the fraction of CSF in pial PVSs which does not enter penetrating PVSs is able to flow back into the SAS and reenter pial PVSs of arteries, forming a kind of recirculation along the surface of the brain. Future studies could test these possibilities. Option (1) seems unlikely because of the good agreement between two independent studies (17 mm/s versus 18.7 mm/s reported by Bedussi et al. (2017) and Mestre et al. (2018b), respectively). Option (2) perhaps plays a role, and future numerical studies with improved modeling of the pial PVS geometry should investigate this possibility. It is likely that option (3) might contribute to the discrepancy, but such experiments are challenging and always have confounding factors. Option (4) may contribute as well; future particle tracking experiments should investigate possibilities (1) and (4).
The values of hydraulic resistance computed with our model can be directly compared to those of prior work. Faghih and Sharp (2018) developed a network model of flow through periarterial spaces and computed a total network resistance of 1.14 mmHg,min/mL. This value is about 2000 times lower than the lowest hydraulic resistance we compute, R = 2300 mmHg,min/mL for the R min scenario. This discrepancy is probably primarily because Faghih and Sharp modeled glymphatic flow in a human, with far more parallel channels than we have considered. Vinje et al. (2020) developed a compartmental model to estimate how elevated intracranial pressure may affect CSF outflow pathways. Although their study modeled human anatomy, they used parameters similar to Intermediate scenario 1 in this study and reported that the hydraulic resistance of the parenchyma was comparable to that of the PVSs, which is in good agreement with our observations.
We find that a substantial fraction of the CSF flowing through penetrating PVSs continues through the parenchyma in every scenario, with values ranging from 32%  (Figures 2J, 2L, 2N, and 2P). Our upper bound is in agreement with the lower bound of flow speeds, 0.083 mm/s, reported by . In addition, our lower bound is in agreement with results from Holter et al. (2017), in which parenchymal flow speed near the outer wall of the PVS is about 0.035 mm/s (see Figure 3 in Holter et al. (2017)). For cases in which the precapillary PVS flow fraction is non-negligible (> 0:5%; R max and Intermediate 1 scenarios), the speeds are also fairly robust, ranging from 2.7 to 20 mm/s (Figures 2J and 2N). This moderate insensitivity to precapillary PVS size (G precap ) -especially for the R max scenario -can be understood as follows: as the crosssectional area A PVS increases, the hydraulic resistance R precap decreases causing the volume flow rate Q to increase, rendering the flow speed ( = Q=A PVS ) approximately constant. To the best of our knowledge, this is the first time precapillary PVS flow speed has been predicted.
We assessed whether each scenario exhibits uniformity in cortical CSF perfusion, which we expect based on reports of tracer penetration below the brain's surface (Gaberel et al., 2014;Ringstad et al., 2018;Eide et al., 2021;Taoka and Naganawa, 2020) and evidence that flow is important for metabolic waste removal (Iliff et al., 2012;Xie et al., 2013;Roberts et al., 2014;Shokri-Kojori et al., 2018;Koundal et al., 2020;Gu et al., 2020;Eide et al., 2021). Our simulations revealed near-perfect cortical CSF perfusion for Intermediate scenario 1, moderately uniform CSF perfusion for the R min scenario, fairly poor CSF perfusion for the R max scenario, and negligible CSF perfusion below the brain surface for Intermediate scenario 2 ( Figure 3). As discussed above, good uniformity in CSF perfusion can be understood as a consequence of scale separation in the hydraulic resistance of the three sequential CSF routes: pial PVSs, penetrating PVSs, and parenchyma/precapillary PVSs ( Figures 3Q-3X). Poor CSF perfusion occurs if these resistances are comparable ( Figures 3Q-3R) or do not increase in the aforementioned order ( Figures 3W-3X). This observation provides an argument in favor of large parenchymal resistance, which could arise because of tight astrocyte endfeet gaps, a low-permeability parenchymal ECS, or a combination of the two. Furthermore, the separation in scale between pial and penetrating PVS resistance ensures that CSF is uniformly perfused from the pial PVSs to the penetrating PVSs. This need for separation in scale may explain why pial PVSs have an oblate shape that minimizes their hydraulic resistance (Tithof et al., 2019).
We performed simulations aimed at capturing the increase in CSF flow during sleep compared to wakefulness (Xie et al., 2013;. Multiple studies demonstrate that glymphatic transport is enhanced under ketamine/xylazine (K/X) anesthesia, resembling natural sleep, and inhibited under isoflurane, resembling wakefulness (Xie et al., 2013;Hablitz et al., 2019;Stanton et al., 2021); indeed, both the prevalence of slow (delta) waves and the ECS porosity under K/X are comparable to natural sleep (Xie et al., 2013). These studies comparing K/X and isoflurane highlight the heterogeneity of tracer transport in different regions of the brain, often with two-to four-fold greater tracer influx under K/X, compared to isoflurane. We found that R min , R max , and Intermediate scenario 2 all exhibit less than a 22% increase in combined volume flow rate during sleep compared to wakefulness ( Figures 4I-4L and 4O-4P); however, we found a 3.2-fold increase in combined volume flow rate for Intermediate scenario 1 with small precapillary PVSs ( Figures 4M). Increased tracer transport can be estimated from increased CSF flow based on the theory of Taylor dispersion (Taylor, 1953;Troyetsky et al., 2021), which describes the effective diffusion coefficient D eff characterizing the rate at which a tracer spreads in a shear flow because of the combined effect of advection and diffusion. For measured pial PVS size and flow speed (Mestre et al., 2018b) and a diffusion coefficient of D = 1310 À 11 m 2 /s (Asgari et al., 2016), D eff =D = 3:8 ( Figure S6A), suggesting pial CSF flow enhances transport 3.8-fold greater than diffusion alone. When the awake-to-sleep volume flow rate is increased less than 22% (R min , R max , and Intermediate scenario 2), the enhanced tracer transport is less than D sleep eff =D awake eff = 32%, whereas a 3.2-fold increase in awake-to-sleep volume flow rate (Intermediate scenario 1 with G precap = 0:07) leads to D sleep eff =D awake eff = 300% ( Figure S6B). Note that prior studies computed enhancement factors based on oscillatory (zero mean) flow (Sharp et al., 2019;Asgari et al., 2016), whereas our calculations are based on steady (nonzero mean) flow, which we have previously argued is more effective for dispersive transport (Thomas, 2019;Troyetsky et al., 2021). Although Taylor dispersion in pial PVSs is unlikely to account for the entirety of tracer transport observed in experiments, these estimates generally suggest that Intermediate scenario 1 with small precapillary PVSs (G precap = 0:07) is the only scenario with sleep/awake variations in volume flow rate large enough to explain tracer transport reported in several experiments (Xie et al., 2013;Hablitz et al., 2019;Stanton et al., 2021 Overall, we find that parameters in the general range of Intermediate scenario 1 will satisfy the majority of experimental observations described in this article. We have found that a network with low PVS resistance (high PVS permeability) and high parenchymal resistance (whether from tight gaps between astrocyte endfeet, low parenchymal permeability, or both) requires a reasonably low pressure drop ( Figure 2F), exhibits nearly perfect cortical CSF perfusion ( Figures 3J and 3L), and -for small precapillary PVSs -most closely captures the observed increase in CSF influx during sleep compared to wakefulness ( Figure 4M). In addition, Intermediate scenario 1 with G precap z0:27 is the only case which exhibits an equal 50/50 flow through precapillary PVSs and parenchyma. It is enticing to speculate that such a parameter regime may enable dynamic regulation of CSF transport; in this scenario, if parenchymal resistance were dominated by astrocyte endfeet, small changes in the endfoot gap could substantially shift CSF perfusion between slower parenchymal flow and faster precapillary PVSs flow. We caution that the parameter space is large, so Intermediate scenario 1 does not provide the only possible case that satisfies the aforementioned criteria, but rather points to a general parametric regime.

Limitations of the study
There are numerous limitations in this study that are noteworthy. Perhaps the most consequential limitation is the uncertainty in several parameters that affect CSF transport through the glymphatic pathway, which we attempted to address by considering different limiting parametric scenarios. We restricted ourselves to a moderate number of cases for the sake of clarity, and we did so by lumping some parameters together, such as the astrocyte endfoot geometry and parenchymal permeability (Table 1). It is important to note that the parametric bounds we employ correspond to extreme values reported in or inferred from the literature (and therefore do not represent strict limits on feasible parameter ranges). Future experimental studies aimed at refining uncertain parameters will be of tremendous value for constructing predictive models. In particular, the hydraulic resistance of gaps between astrocyte endfeet is especially uncertain, with our estimates here ranging over almost seven orders of magnitude ( Figure 3Y). The low end of this range suggests the astrocyte endfeet play no role in limiting CSF transport from the penetrating PVS to the parenchymal ECS ( Figures 3S-3T and 3W-3X), whereas the upper limit has hydraulic resistance comparable to that of the parenchymal ECS ( Figures 3Q-3R, 3U-3V), suggesting the astrocyte endfeet play a critical role. The bounds for endfeet gap sizes that we consider (20 nm-5.1 mm) come from studies based on chemical (Mathiisen et al., 2010) and cryo (Korogod et al., 2015) fixation. The former exhibits overlapping between the endfeet which is likely a consequence of shrinkage of the PVS during the fixation process (Mestre et al., 2018b), casting some doubt on this lower bound. However, a recent study (Wang et al., 2021) reported heterogeneity in the size of astrocyte endfeet, with larger endfeet (i.e., smaller and/or fewer gaps) surrounding larger vessels. This observation provides a mechanism that improves the uniformity of cortical CSF perfusion if the endfoot gaps are small enough to provide resistance comparable to that of the parenchymal ECS. Future studies should characterize astrocyte endfeet gap dimensions, ideally in vivo.
In our model, CSF flow is driven by the simplest possible mechanism -an externally applied pressure drop across the entire network. However, other potential driving mechanisms (e.g., pressure gradients generated by arterial pulsations (Mestre et al., 2018b), functional hyperemia (Kedarasetti et al., 2020), or osmotic effects Halnes et al., 2019;Rasmussen et al., 2021)) could be tested with this network model approach by implementing pressure sources (i.e., ''batteries'') throughout the network. In particular, incorporation of osmotic effects could be leveraged to investigate the mechanisms by which aquaporin-4 facilitates glymphatic flow (Iliff et al., 2012;Asgari et al., 2015;Mestre et al., 2018a;Hablitz et al., 2020), although there is some debate about this point (Smith et al., 2017;Mestre et al., 2018a). In this study, we have chosen an applied pressure drop in each scenario such that the median pial PVS flow speed matches the average value measured in experiments (Mestre et al., 2018b;Bedussi et al., 2017); future experiments that obtain flow speeds at multiple PVS locations will prove useful for validating the accuracy of our model. Yet another important limitation to our approach, already touched on in the fourth paragraph of the Discussion involves the connectivity of pial PVSs at the surface of the brain. By introducing ''short-circuit'' connections between PVSs of pial arteries and pial veins, our model could be adapted to estimate the fraction of CSF that continues along the surface of the brain (and/or through stomata) versus the fraction that continues through deeper PVSs and the parenchyma. Such a model would greatly benefit from experimental estimates of how many such connections typically exist. Finally, we highlight that our model can be generalized to predict transport of dye, metabolic waste, drugs, or any other molecules because of advection-diffusion. Such future studies will contribute to the substantial ongoing debate regarding the nature of transport in penetrating PVSs (Asgari et al., 2016;Sharp et al., 2019;Kedarasetti et al., 2020 In future work, we intend to implement numerous refinements to our simulation, but many will likely offer improvements that are of secondary importance compared to obtaining better estimates of critical parameters (as discussed above). The idealized geometry we have adopted has a regular, repeating structure composed of four different types of homogeneous channels and consequently lacks the high spatial variability characteristic of the true network. Future models could use randomly sampled statistical distributions to assign geometric parameters (Blinder et al., 2010(Blinder et al., , 2013 or directly implement the geometry of a synthetic (Linninger et al., 2019) or real (Kirst et al., 2020;Mestre et al., 2020) vascular network. We restricted our model to the arterial side of the network whereas we relied on assumptions about PVSs at the capillary and venous level to enable lumped modeling, but future studies could include substantially greater detail.
In this study, we predicted CSF transport throughout a mouse brain, but our network could be expanded to model a human brain by adding more vascular generations. Such an approach would be more challenging because of the fewer measurements available for constraining the parameter space in humans, compared to mice. However, many parameters may be conserved across species (e.g., porosity, PVS area ratios, endfoot gap size). Development of such a model has tremendous clinical value, as it could offer insight into a myriad of neurological disorders. Conditions such as Alzheimer's disease, traumatic brain injury, and subarachnoid hemorrhage are all known to coincide with disrupted glymphatic transport (Rasmussen et al., 2018).

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: d KEY RESOURCES where ε is the porosity, t is the tortuosity, and S is the specific surface area for a porous medium (Hommel et al., 2018). Xie et al. (2013) reported an increase of ε from 0.14 during wakefulness to 0.23 during sleep, with no change in tortuosity. Assuming S remains approximately constant, this suggests k sleep par =k wake par = 5:5.