Optimal Mass Transport with Lagrangian Workflow Reveals Advective and Diffusion Driven Solute Transport in the Glymphatic System

Abstract

The glymphatic system (GS) hypothesis states that advective driven cerebrospinal fluid (CSF) influx from the perivascular spaces into the interstitial fluid space rapidly transport solutes and clear waste from brain. However, the presence of advection in neuropil is contested and solutes are claimed to be transported by diffusion only. To address this controversy, we implemented a regularized version of the optimal mass transport (rOMT) problem, wherein the advection/diffusion equation is the only a priori assumption required. rOMT analysis with a Lagrangian perspective of GS transport revealed that solute speed was faster in CSF compared to grey and white matter. Further, rOMT analysis also demonstrated 2-fold differences in regional solute speed within the brain. Collectively, these results imply that advective transport dominates in CSF while diffusion and advection both contribute to GS transport in parenchyma. In a rat model of cerebral small vessel disease (cSVD), solute transport in the perivascular spaces (PVS) and PVS-to-tissue transfer was slower compared to normal rats. Thus, the analytical framework of rOMT provides novel insights in the local dynamics of GS transport that may have implications for neurodegenerative diseases. Future studies should apply the rOMT analysis approach to confirm GS transport reductions in humans with cSVD.

Introduction

The glymphatic system is described as a perivascular transit passageway for cerebrospinal fluid (CSF) for exchange with interstitial fluid (ISF), thereby facilitating waste drainage from the brain^1,2. Investigations of glymphatic system (GS) function have escalated given its important role in Aβ¹ and tau³ clearance from brain and the inferred implication for neurodegeneration, including Alzheimer’s disease^2,4,5,6,7. The GS is made up by the perivascular spaces (PVS), which connect with ISF via the aquaporin 4 (AQP4) water channels on astrocytic end-feet and through the small gaps between the overlapping astrocytic end-foot processes¹. The GS hypothesis states that advective CSF influx from the PVS rapidly drives interstitial solutes and waste products out via peri-venous channels^1,8. Although solute transport in the PVS along pial arteries on the surface of the brain is advective (bulk flow) and driven by cardiac pulsatility^1,9,10,11, the presence of advective streams in parenchyma^12,13,14 is contested with the argument that advection does not occur in the neuropil and solutes are transported by diffusion only^{15,16,17,18,19,20,21}. No single study has unequivocally determined which mode of solute transport prevails inside the brain proper. This gap in knowledge has significantly impeded exploration of the GS and has raised questions regarding its importance for brain health, disease, and its potential as a therapeutic target against neurodegeneration.

Recently, various computational analyses and simulations have been performed to better understand the driving mechanisms of solute transport across PVS boundaries and through ISF space^22,23,24. Conservation principles are typically used to derive these simulations, a popular one being the porous media model^25,26,27. Ray et al.²⁸ proposed adding an advection term to the (diffusive) porous media model and their simulations involving molecules such as Aβ indicated the presence of advective transport in the ISF. Simulations are very effective for distinguishing between advection and diffusion, but when applied in vivo, inferable insight becomes limited. The key reason for this problem lies in the fact that modeling in live brain requires accurate delineation of anatomical boundaries and knowledge of various kinematic parameters, many of which are unknown.

To provide local dynamic analysis and inform on GS transport modes across compartments in the live brain we now introduce novel computational methods for studies based on dynamic contrast enhanced magnetic resonance imaging (DCE-MRI)^8,29,30. Specifically, we implemented techniques from the theory of optimal mass transport (OMT)^31,32,33. The OMT theory is centuries old and has impacted many topics of research in the physical sciences, economics, computer science and image processing^31,34,35,36. Given initial and final mass distributions, OMT provides a powerful dynamic computational fluid formulation that gives an optimal interpolation path of minimal energy among all possible interpolations that preserve mass³⁷. In this way, OMT can be formulated as a variational problem based on the principle of least action. A particularly attractive feature of OMT is that one obtains unique solutions to the model, not by imposing smoothness, but rather by seeking flows that minimize the expended transport energy. In the classical formulation of OMT, the continuity equation involves advection only. Considering that diffusion is assumed to be always occurring in the brain, we have added a diffusion term in the work presented here. We refer to our modified OMT formulation with the advection/diffusion equation as the regularized OMT problem (rOMT). We utilized an associated Lagrangian formulation of rOMT for the construction of ‘pathlines’ to effectively extract and visualize GS transport flows over a predetermined set of time frames in one comprehensive figure. Time-varying particle (a.k.a. solute) attributes associated with the pathlines, such as ‘speed’ were also computed.

Here we apply our new rOMT formulation to address three overarching and unresolved questions: 1) Can we confirm the existence of the two types of transport (diffusion and advection) in the different GS compartments? 2) Can we obtain voxel-level information on solute transport across all tissue compartments, circumventing the problem of varying GS transport kinetics and avoiding the need to specify all of the material properties for simulations? and 3) Given that GS dysfunction is reported in neurodegeneration particularly where vascular dysfunction may be involved^7,38,39, can our rOMT framework dissect different modes of solute transport in an animal model of cerebral small vessel disease?

Results

Introducing rOMT for tracking advection and diffusion modes of GS transport

GS transport was measured in live rats using DCE-MRI and administration of gadoteric acid (Gd-DOTA) into CSF via the cisterna magna^29,30,40. The DCE-MRI rat brain data were input into the rOMT framework for dissecting GS transport modes of the Gd-DOTA solute. The rOMT formulation with the inclusion of both advection and diffusion terms in the constraint (continuity) equation is described in detail in Methods. Here we highlight that the Lagrangian formulation was used to construct dynamic ‘pathlines’ for visualizing GS transport flows in one comprehensive figure, derived from the rOMT returned velocity field and interpolated images (Methods and Fig. 1). As such, a Lagrangian pathline traces the trajectory of a specific particle over a pre-defined time interval. A particle may refer to a parcel of mass or an individual substance. For our purposes, we use particle interchangeably with solute. Observations of particle attributes, such as speed, are made from the vantage point of the particular particle as it traverses its trajectory. In this way, ‘pathlines’ and ‘speed-lines’ encapsulate the moving particle dynamics in the GS network in a single entity. Additional information, such as flow volume (size of the pathline network), reflecting the total flux of moving particles, is extracted and used in the following for analyzing GS transport across tissue compartments and brain regions. For ease of narration, we refer to the extraction and visualization of the Lagrangian representation of Glymphatic Dynamics as our GLaD framework (Fig. 1). rOMT as well as GLaD analysis was performed on DCE-MRI images taken over the 120 min interval starting at the time of peak signal (Fig. 1a). This time interval afforded the best representation of GS transport because the signal-to-noise ratio (peak) and redistribution of Gd-DOTA tracer into brain parenchyma is maximized.

Figure 1

Processing steps for rOMT Lagrangian Glymphatic System transport flow derivatives. Illustration of the regularized optimal mass transport (rOMT) and Lagrangian representation of Glymphatics Dynamics (GLaD) pipeline for visualizing transport flows (a) DCE-MRIs over a 120 min experimental window (50–160 min) were fed as inputs into the rOMT model (b) The model returns the interpolated density images and velocity fields which are (c) subsequently processed by our GLaD visualization framework. (d) Advective and diffusive glymphatic system transport is expressed by the pathlines captured in the Lagrangian analysis with corresponding dynamic speedlines. Binary pathlines were used to extract flow volume and speedlines were used to assess spatial distribution of speed trajectories in whole brain as well as to derive mean speed values. (e) Thereafter, compartmental analysis is carried out based on tissue masks derived from the morphometric analysis.

Note also that the rOMT algorithm captures the Fick’s law aspect of the flow model, namely the diffusion mode of particle transport controlled by the gradient of the concentration. This was accomplished using a custom built ‘diffusion’ phantom incorporating the microdialysis technique⁴¹ (Supplementary Fig. 1 and Movie 1). The rOMT model was applied to the DCE-MRI data acquired in the diffusion phantom. Knowing that the solute transport in the phantom was uniquely diffusive, we examined how our model performed if we forced fictitious advective transport. This was done by using very low diffusivity values that were intentionally insufficient to mimic the diffusing intensity patterns observed in DCE-MRI acquired in the phantom. This numerical intervention thereby forced advection to compensate for the lack of diffusion in the modeling of interpolated images. Although the advective and diffusive fluxes influence each other (see Methods), the advective flux is primarily influenced by time-varying changes in voxel intensity while the diffusive flux is primarily influenced by spatial variations. Attempting to account for both temporal and spatial changes in intensity with advection alone resulted in an irregular flow field in the interpolated images. Therefore, it was not surprising that efforts to extract meaningful (smooth) pathlines with our aforementioned GLaD framework failed on the phantom data under these conditions.

Next, we repeated the rOMT model on the DCE-MRI diffusion phantom data, but this time with an appropriate diffusivity value (experimentally determined). The resulting flow field was governed by the diffusive flux, as expected. As shown (Supplementary Movie 1), diffusive flux vectors point from voxels with higher intensity toward voxels with lower intensity and there is near-perfect consistency between the acquired DCE-MRI images (Fig. 2a) and the returned interpolated model images (Fig. 2b). This numerical experiment proved informative for understanding the advantages of adding the diffusion term to the classical dynamic OMT formulation, and the roles of each transport mode (see Methods).

Figure 2

Diffusion is required for rOMT model data to capture glymphatic transport patterns. (a) Agar distribution after ~ 2 hrs of diffusion, as captured with contrast-enhanced MRI, is shown in grayscale. The diffusive direction, proportional to the negative of the image gradient, is shown by the green arrows. (b) The corresponding rOMT-returned image is shown in grayscale with diffusive flux vectors characterizing the flow behavior and in alignment with Fick’s law are overlaid in green for comparison to ‘ground truth’ (a). (c) Whole brain % signal from baseline at 100 min after contrast injection is shown for context. Corresponding color-coded flow speeds along streamlines computed from the rOMT-derived vector field associated with 100 min after contrast injection for multiple diffusion levels are compared, illustrating how an appropriate diffusion value was selected. (d) With no diffusion, the flow shows exaggerated peripheral advective behavior and central parenchymal transport is unaccounted for. (e) Using moderate but sufficient diffusion yields feasible peripheral behavior and accounts for transport in deeper parenchymal regions. (f) With high diffusion, more activity is accounted for but exaggerated peripheral behavior is also seen.

Diffusion and advection are both required for rOMT modeled images to match in vivo images

Adding diffusion in the rOMT model was required for matching GS transport patterns observed in the live rat brain by DCE-MRI (Fig. 2c) with the modelled rOMT images (Fig. 2d–fs). To choose an appropriate diffusivity value, we tested multiple values and examined how the flow fields changed. Specifically, we computed streamlines from the rOMT derived velocity fields at each time step along with the corresponding speeds. Streamlines are curves that are tangent to the velocity field at a fixed time, informing on the collective instantaneous behavior of the flow. We chose this approach over the GLaD-pathline analysis to investigate the ‘smoothness’ of the flow field at individual time steps, which should be directly affected by diffusion. A representative example of speed (color-coded) from a normal Wistar Kyoto (WKY) rat associated with streamlines is shown in Fig. 2d–f, illustrating the effect of increasing the strength of the diffusion term ‘σ²’ in the rOMT algorithm. With either no diffusion (Fig. 2d) or too much diffusion (Fig. 2f), the captured behavior in the CSF compartment is erroneously overactive while information about the flow in the deeper regions, where parenchymal transport of interest occurs, is altogether unrealized. We determined that an intermediate amount of diffusion (σ² = 0.002) was required (Fig. 2e) to best replicate observed GS transport in both peripheral and central regions (see the Methods section for a more detailed discussion). This led us to conclude that parenchymal GS transport is governed by both advection and diffusion. However, it is important to emphasize that parenchymal GS transport of the extracellular tracer, Gd-DOTA captured within a given DCE-MRI defined voxel (e.g., range of 0.003–0.027 mm³) represents transport within an ISF compartment of 18–22%^42,43 in addition to the PVS spaces. The volume fraction of (abnormally) enlarged PVS by MRI is reported in the range of ~2–3% in the basal ganglia of human brain⁴⁴ and we therefore estimate the PVS volume fraction in normal rodent brain to be of the order of 1% or less.

We should note that the speed pathlines do not address the issue of ‘concentration’ and over the time window, there will be declining contrast in both compartments. We are not assuming that the concentration is the same in the two compartments when we apply the algorithm to the data. The data in its current form does not allow us to distinguish concentrations. Further, when we implement the mathematical algorithm, the contrast of the images is normalized so that the total contrast (sum of the intensities) is the same, even though locally there are of course changes. We are now working on employing an unbalanced version of the algorithm³⁶, for which the normalization (nor any assumption of mass preservation) will be necessary.

We also note that future work will entail improvements in the diffusion model. At the current stage of method development, we have chosen the σ² value empirically based on ‘best’ match of simulated rOMT data to raw DCE-MRI data. We plan to train our method on a number of GS transport samples and then use standard machine learning techniques to design a more automatic method for tuning σ² as well as consider possible spatially dependent diffusion covariance matrices.

Solute speed varies across tissue compartments in normal and hypertensive rats

To derive GS solute transport across tissue compartments at the voxel-level, we applied the proposed rOMT + GLaD analytical pipeline to segmented whole brain DCE-MRIs acquired in normotensive WKY rats and in spontaneously hypertensive stroke prone (SHRSP) rats (Methods). The SHRSP rat is a naturally bred animal which mimics human hypertension and associated cerebral small vessel disease (cSVD)⁴⁵. Examples of pathlines (color-coded for speed) derived in CSF, WM and grey matter (GM) compartments from a WKY and a SHRSP rat are shown in Fig. 3a–c. Anatomically, pathlines in GM were more prominent in the cerebellum, hippocampus and hypothalamus of WKY rats (Fig. 3b). In WM of WKY rats, pathlines dominated in the brainstem and midbrain (Fig. 3b). In SHRSP rats, subarachnoid CSF pathlines showed overall more heterogeneous distribution and were most prominent ventrally (Fig. 3c, top). In WKY rats, total pathline volume in the CSF compartment (representing the total flux of ‘moving’ solutes in CSF over 120 min) was 20% higher when compared to SHRSP rats (7124 ± 605 voxels vs 5822 ± 636 voxels, p = 0.001, Mann-Whitney, Fig. 3d), suggesting restricted CSF transport in SHRSP compared to WKY rats. In SHRSP rats, the lower than normal CSF pathline volume was associated with reduced pathline volume in both GM and WM compartments when compared to WKY rats (Fig. 3d), signifying overall decreased GS transport in SHRSP rats. In the corpus callosum, pathline volume was variable with a trend towards decreased solute transport in the SHRSP compared to WKY rats (582 ± 236 vs 330 ± 120, p = 0.051, Mann Whitney).

Figure 3

rOMT whole brain analysis reveals varying solute speed across tissue compartments. (a) Cerebrospinal fluid (CSF), grey matter (GM) and white matter (WM) masks overlaid on the brain anatomical template to illustrate the spatial distribution of the brain tissue compartments used for the rOMT analysis. (b) Pathlines - color-coded to denote speed - overlaid on corresponding anatomical brains are shown from a WKY rat in CSF, WM and GM compartments. Red and blue colored pathlines denote fast and slow pathline speeds, respectively. In the CSF compartment, pathlines are fast moving at the ventral surface. Pathlines with slower speed lines can be seen inside the cerebral ventricles (cVtr). Speed pathlines in GM and WM compartments are variable: faster speed in the periphery and slower speed profiles as they penetrate into deeper parts of parenchyma. CTX = cortex; CB = cerebellum; Hip = hippocampus; HT = hypothalmus; MB = midbrain; BS = brain stem. (c) Corresponding pathlines in the three tissue compartments from a SHRSP rat are shown. It is clear that pathlines in CSF are sparser and more heterogeneously distributed in the SHRSP rat when compared to the WKY rat. For example, there are no pathlines above cortex and very few superficial to the cerebellum. For SHRSP, there are also less pathlines in the GM and WM, as compared to the WKY rat. (d) Analysis of the total flux of solutes in the pathline network (number of voxels) revealed decreases in all three compartments in SHRSP when compared to WKY rats. (e) In both WKY and SHRSP strains, mean pathline speed in CSF was significantly faster than mean pathline speed in the WM compartment. Data are mean ± SD. Mann-Whitney’s test was used for cross-strain analysis. Friedman’s ANOVA was used for within-strain regional analysis. *significant at level alpha = 0.05.

We next analyzed the time-varying particle attributes associated with the pathlines, that is, ‘speed’ given by the magnitude of the transport vectors derived by the rOMT + GLaD procedure (Fig. 1). Speed trajectories of the pathlines are displayed as color-coded maps for each tissue compartment (Fig. 3b,c). In WKY rats, fast speed (red color) was evident in the CSF subarachnoid compartment on the ventral surface of the brain, superficial to the cerebellum and cerebral cortex (Fig. 3b, top). Slower (blue color) speed trajectories were observed inside the cerebral ventricles (Fig. 3b top). In WKY and SHRSP rats, mean pathline speed was ~40% faster in CSF when compared to the GM and WM compartments (Fig. 3e). Mean pathline speed in CSF of SHRSP rats was within a similar range to those of WKY rats (0.10 ± 0.026 versus 0.09 ± 0.02, p = 0.558) and no differences in mean speed were observed in GM and WM compartments between the two strains. We also compared directional tracer movement from CSF into the tissue compartment between WKY and SHRSP rats (Supplementary Fig. 2). Using the rOMT-derived velocity field, we evaluated the total (advective + diffusive) flux across this surface with attention to the net mass moving from the CSF compartment into the tissue compartment. This analysis showed that the average net mass transported into the tissue compartment was significantly higher for WKY when compared to SHRSP rats (Supplementary Fig. 2).

Solute speed varies across brain regions

We next considered regional GS transport using GLaD analysis applied to DCE-MRI data extracted from pre-selected brain regions distributed in the GS network including the hippocampus, cerebellum, superior colliculus and basal forebrain regions (Fig. 4a). Using the same visualization scheme as for the whole brain data, Fig. 4 shows pathlines color-coded for speed captured inside cerebellum, superior colliculus, hippocampus and basal forebrain from a representative WKY (Fig. 4b) and SHRSP rat (Fig. 4c). In the SHRSP rats, the total volume of pathlines was significantly reduced in the cerebellum, hippocampus and superior colliculus when compared to WKY rats (Fig. 4d). Furthermore, GS transport varied across brain regions as evidenced by the varying distribution of pathlines across the brain in both strains. In the SHRSP rats, pathlines were sparser in the hippocampus compared to the basal forebrain (Fig. 4d). For both WKY and SHRSP rats, it was evident that pathline speed was not uniform across regions (Fig. 4e). Within normotensive WKY rats, mean pathline speed in the basal forebrain region was significantly increased when compared to the hippocampus (0.102 ± 0.025 vs 0.055 ± 0.016, p = 0.006, Friedman test) and superior colliculus (0.102 ± 0.025 vs 0.053 ± 0.012, p = 0.021, Friedman test). Within the SHRSP rats, mean pathline speed also varied significantly across regions (Fig. 4f). Across strains, the analysis revealed that mean speed in the superior colliculus was significantly reduced in the SHRSP when compared to WKY rats (0.053 ± 0.012 vs 0.040 ± 0.005, p = 0.022, Mann Whitney).

Figure 4

Solute speed is not uniform across brain regions. (a) Volume rendered anatomical delineation of brain regions including the cerebellum (Cb), superior colliculus (Sc), hippocampus (Hip) and basal forebrain (Bf) that were used for the analysis. Top and lateral views of the regions are shown. (b) Color-coded pathlines depicting speed captured over 120 min are shown from a normal WKY rats inside the regions of interest. In the normal WKY rat, speed trajectories in the cerebellum penetrate the entire region. Pathlines at the level of the dorsal hippocampus are sparser when compared to the ventral hippocampus. (c) Color-coded pathlines depicting speed captured over 120 min are shown from a SHRSP rat for the same regions. Note that in this particular SHRSP rat, pathlines barely reach the superior colliculus and are not observed in the center of the cerebellum. (d) The total pathline volume for each region was calculated and compared between WKY and SHRSP rats. Except for the Bf, all regions showed lower pathline volume in SHRSP when compared to WKY rats. (e) Mean pathline speed varied across brain regions in both strains. The fastest mean speed was observed in the basal forebrain for both strains. Data are mean ± SD. Mann-Whitney’s test was used for cross-strain analysis. Friedman’s ANOVA was used for within-strain regional analysis. *Significant at level alpha = 0.05.

Solute transfer from PVS-to-tissue is impaired in hypertensive SHRSP rats

To further explore GS solute transport dynamics between the two strains, we acquired DCE-MRIs at higher spatial resolution using a smaller 1-cm RF-surface coil positioned above the left hemisphere of the rat’s head to capture transport in the PVS of pial arteries at the circle of Willis and along the entire middle cerebral artery (MCA) of the left hemisphere (Movie 2). To boost signal-to-noise ratio in the PVS, we administered 20 µL of 1:5 Gd-DOTA diluted with sterile water at a rate of 1.5 µL/min into the CSF. The higher Gd-DOTA levels in CSF enabled us to track solute transport along the entire left hemispheric MCA at a voxel resolution of 0.15 × 0.15 × 0.15 mm³. We first extracted time signal curves (TSC) from the PVS along the MCA on the ventral surface (a.k.a. ‘MCA-PVS root’, Fig. 5a) and observed strikingly different signal dynamics between SHRSP and WKY rats (Fig. 5b). While the MCA-PVS root TSCs from SHRSP and WKY rats were characterized by similar time-to-peak and peak magnitude, the signal ‘relaxation’ phase was faster in WKY rats, suggesting faster transit of solutes into parenchyma compared to SHRSP rats (Fig. 5b). This interpretation was supported when comparing DCE-MRIs at the base of the brain after ~3 hrs of CSF Gd-DOTA circulation (Fig. 5c,d): Significantly higher Gd-DOTA signal was evident along the root PVS-MCA of SHRSP rats (Fig. 5d) but not along the MCA distal to this point when compared to WKY rats (Fig. 5e), suggestive of physical obstruction (or ‘slow-down’) of MCA-PVS flow in this area.

Figure 5

Solute transfer from the perivascular space (PVS) along the middle cerebral artery. (a) Contrast enhanced MRIs from a WKY rat at the level of the circle of Willis near the take-off of the middle cerebral artery (MCA). Contrast in the perivascular space along the circle of Willis and MCA is evidenced as a high intensity signal along the vasculature, which is black. Regions of interest along the MCA (white lines) were used to capture time signal curves (TSC) from the dynamic contrast enhanced MRIs acquired in WKY and SHRSP rats. (b) Mean TSC extracted from WKY (N = 4) and SHRSP (N = 4) rats along the first segment of the MCA (‘root’) are displayed. Data are mean ± SD. It is evident that while ‘time-to-peak’ and peak magnitude are within same range for the two strains, the signal decline (representing redistribution of contrast from PVS to tissue and tissue clearance) is strikingly faster in WKY when compared to SHRSP rats. (c) Volume-rendered, color-coded signal map after ~180 min of circulation of CSF contrast in a WKY rat overlaid onto the corresponding anatomical template. High and low contrast signal amplitude is given by red and blue colors respectively. Contrast is observed along the entire MCA during its trajectory across the left hemisphere. (d) Corresponding color-coded signal map from a SHRSP rat showing much higher levels (red color) of contrast in the PVS along the circle of Willis and along the MCA root. Heterogeneously distributed contrast along the MCA is more evident in the SHRSP when compared to WKY rat. (e,f) The color-coded speed maps along the MCA are overlaid on the anatomical volume rendered MRIs of the left hemisphere of a WKY (e) and SHRSP (f). While pathline speed appeared to be homogeneous along the MCA in WKY rats, pathline speed in SHRSP rats was high in areas localized to the ‘root’ area and lower distal to this point.

The rOMT + GLaD pipeline was applied to extract pathline speed information along the MCA. Pathline start points were selected from binary masks capturing tissue along the MCA in order to analyze transport behavior stemming from this specific region (see Methods). The color-coded speed maps along the MCA were overlaid on the anatomical MRIs of the left hemisphere (Fig. 5e,f). While pathline speed appeared to be homogeneous along the MCA in WKY rats (Fig. 5e), pathline speed varied along the MCA in SHRSP rats (Fig. 5f). Specifically, in SHRSP rats, pathline speed was increased by 3-fold along the MCA ‘root’ when compared to speed along the MCA distal to this area (0.12 ± 0.04 vs 0.04 ± 0.01, p = 0.029, Student’s t test). Corresponding measurements in WKY rats did not reveal significant differences in speed along the MCA (p = 0.25). We also noted that in WKY pathlines had migrated >4 mm away from the MCA-PVS into the adjacent tissue (Fig. 5e). In contrast, in SHRSP rats, pathlines were observed in closer proximity to the PVS-MCA indicating slower transfer of solutes into the ISF (Fig. 5f).

We also performed directional analysis of the advective flux direction along the MCA in WKY rats (Fig. 6). Advective flux vectors were calculated across the boundaries of the binary mask created along the MCA from the rOMT output. Figure 6 shows advective flux vectors (normalized to have unit length) along the MCA overlaid on corresponding anatomical MRI templates from WKY (Fig. 6a,b) and SHRSP (Fig. 6c,d) rats. Fig. 6a–d highlight directions of the advective flux vectors from the ventral surface of the brain at the level of the circle of Willis from where the MCA takes off (white boxes in Fig. 6a,c). In the WKY rat the advective flux vectors’ direction is following the path along the MCA in a congruent pattern (Fig. 6a,b). However, in the SHRSP rat, the advective flux vectors are following the direction of the MCA only part of the way. Specifically, at the inflection point where the MCA change direction the advective flux vectors deviate from the PVS along the MCA (Fig. 6c,d). These deflections in advective flux vector directions along the MCA in SHRSP compared to WKY rats is also clearly appreciated in caudal projection views (Fig. 6e,f).

Figure 6

Advective flux vectors along the MCA show aberrant flow in SHRSP. (a,b) Advective flux vectors across the boundary along the MCA, determined from the binary mask, were calculated from the rOMT output. The advective flux direction, determined by normalizing the advective flux vectors to have unit length, are shown overlaid on an anatomical mask of a WKY rat. The ventral surface at the level of the circle of Willis and MCA origin is highlighted. In the WKY rat, the directions of the advective flux vectors follow the path along the MCA in a continuous pattern. (c,d) In a corresponding view of a SHRSP rat, the advective flux vectors are pointing in the same direction along the first segment of the MCA along the ventral surface. However, distal to this point, the vectors start diverging away from the MCA. (e,f) The deflections in the direction of the advective flux vectors at the transition from ventral to lateral segments of the MCA in the SHRSP rat compared to the WKY rat is also clearly visible in ‘caudal’ projection views.

Mechanisms underlying the decreased GS transport in SHRSP compared to WKY rats

We next explored major drivers of GS transport to uncover possible mechanisms underlying the differences observed between SHRSP and WKY rats. Increases⁹ and decreases^1,10 in vascular pulsatility have been shown to accelerate and decelerate GS influx, respectively. Only subtle differential effects of rat strain on the mean heart rate recorded continuously during the glymphatics MRI experiments were noted (Supplementary Table 1). In fact, mean heart rate was slightly increased in SHRSP compared to WKY rats, leading us to conclude that changes in heart rate were unlikely to contribute to GS transport impairment observed in the SHRSP rats. In a series of randomly selected awake, non-anesthetized rats from age-matched cohorts of WKY and SHRSP rats, non-invasive blood pressure was measured using the tail-cuff method. Systolic and pulse pressures were significantly higher in SHRSP compared to WKY rats (Supplementary Table 1). Anesthesia generally lowers blood pressure when compared to the awake state and the anesthetized SHRSP rats used for the glymphatics experiments might therefore not have been hypertensive. Thus, in a separate series of rats with femoral artery catheterization, we documented that mean arterial blood pressure indeed remained higher in SHRSP rats, even under anesthesia, when compared to WKY rats (Supplementary Table 1). Based on these data, we concluded that higher than normal pulse pressures in SHRSP rats might have contributed to the impaired GS transport observed in our study. However, the anesthetic regimen did lower the blood pressure in the SHRSP rats when compared to non-anesthetized SHRSP rats, and it is possible that this change in blood pressure might have blunted the pathophysiological consequence of chronic hypertension on GS transport.

State of arousal⁴⁶ and type of anesthesia⁴⁷ have also been shown to affect GS transport. We did not measure electroencephalogram in the experiments, however anesthetic depth was evaluated indirectly by quantifying anesthetics administered and was found to be similar between strains (Supplementary Table 1). Therefore, we concluded that differences in depth of anesthesia were unlikely to have contributed to the GS transport alterations observed.

Polarized perivascular expression of AQP4 water channels has been associated with rapid transport of small molecular weight solutes from PVS to brain ISF^1,39,48. We performed immunohistochemistry on a smaller series of WKY (N = 4) and SHRSP (N = 4) rats to examine the AQP4 expression patterns associated with small vessels and capillaries in select brain regions. In WKY rats, the distribution of AQP4 immunofluorescence appeared homogeneously distributed in the occipital cortex (Fig. 7a) and ventral hippocampus (Fig. 7c). In SHRSP rats, increased immunofluorescence was noted in patchy areas associated with penetrating cortical and hippocampal vessels (Fig. 7b,d). AQP4 immunofluorescence intensity was measured along linear ROIs crossing the lumen of small vessels or cerebral capillaries (separately) and converted into a polarization index (Supplementary Fig. 3 and Methods). The AQP4 polarization index for capillaries was of the same order of magnitude as that of small vessels. This analysis showed that AQP4 immunofluorescence surrounding the small vessels and capillaries was more polarized in WKY (higher polarization index) when compared to SHRSP rats (Fig. 7g,h). Further, the AQP4 polarization index for capillaries in ventral hippocampus of SHRSP rats was significantly lower when compared to WKY rats (WKY AQP4 polarization index: 2.55 ± 0.08 versus SHRSP AQP4 polarization index: 2.28 ± 0.058, p-value = 0.0079). We also executed qPCR to measure the mRNA expression levels of GFAP and AQP4 in a separate series of WKY (N = 8) and SHRSP (N = 8) rats in the following regions: cerebellum, cortex, hippocampus, striatum and thalamus. This analysis revealed no statistical differences in AQP4 or GFAP mRNA expression between WKY and SHRSP rats (Supplementary Fig. 4). Collectively, based on these data, we conclude that impaired perivascular AQP4 polarization in SHRSP rats (but not increased total tissue AQP4 expression or astrogliosis) might have contributed to the reduction in PVS-to-brain solute transfer observed in the SHRSP rats when compared to WKY rats.

Figure 7

Perivascular AQP4 water channel expression is altered in SHRSP rats. (a) Confocal representative micrographs show AQP4 immunoreactivity extending from the occipital cortical surface towards the corpus callosum (CC) and shows that AQP4 was distributed uniformly. (b) Equivalent confocal montage from the ventral hippocampus of a WKY rat. AQP4 is distributed homogeneously in the Stratum Oriens (Or) and Stratum Radiatum (SR). (c) AQP4 distribution in occipital cortex from a SHRSP rat showing higher intensity expression along penetrating vessels (white arrows). (d) In SHRSP rats, intensely expressed AQP4 can be noted in patchy areas surrounding vessels in the SR (arrows) and occasionally in the hilus (H) region. (e,f) AQP4 intensity was measured along linear ROIs crossing the lumen of small vessels or cerebral capillaries in WKY and SHRSP rats and converted into a polarization index (Supplementary Fig. 4 and Methods for more information). (g) Perivascular AQP4 polarization index of small vessels in the ventral hippocampus is significantly reduced in SHRSP when compared to WKY rats. (h) In the occipital cortex, the perivascular AQP4 polarization index of small vessels was also significantly reduced in SHRSP when compared to WKY rats. Scale bar = 250 µm.

Finally, we explored vascular and structural anatomy at the level of the circle of Willis and area of the MCA root where GS transport in the PVS along the MCA appeared to be ‘obstructed’. Supplementary Fig. 6 shows representative brain sections from WKY and SHRSP rats that were immunolabeled with an antibody to collagen IV to identify cerebral vessels and micro-vessels. In WKY rats, smaller networks of micro-vessels were observed in distinct areas at the base of the brain surrounding penetrating arteries (Supplementary Fig. 6d–f). However, in SHRSP rats, corresponding microvascular networks were larger in size and the microvasculature were denser (Supplementary Fig. 6g–i). We speculate that the microvascular network remodeling might have affected the PVS transport from the circle of Willis into downstream PVS conduits including MCA. However, more investigations are needed to fully understand mechanisms underlying the ‘obstructed’ PVS flow in SHRSP rats compared to WKY rats.

Discussion

Our rOMT Lagrangian workflow method applied to DCE-MRI images of GS transport in whole rat brain provides evidence of both diffusion and advective solute transport modes in the brain parenchyma. First, in normal brains, rOMT analysis revealed that both advection and diffusion terms were required for pathlines with trajectory speed in the modeled images to resemble GS transport patterns in live brain. Second, particle speed varied strikingly across tissue compartments with faster solute passage in CSF compared to GM and WM compartments. Third, large differences in regional solute speed within the brain were also documented. Specifically, in regions close to the ventral surface near large pulsatile arteries (e.g. basal forebrain) solute speed was 2-fold higher when compared to midbrain and hippocampal regions. Collectively, these novel in vivo results demonstrate that advective transport is dominating in the non-cellular CSF spaces and that slower diffusion transport contributes to GS transport as the solute is transferred into parenchyma. Differences in solute speed between CSF and parenchyma can be explained by more restricted transport in the smaller ISF space⁴⁹ when compared to non-cellular CSF spaces, given that the solute used in the DCE-MRI studies (Gd-DOTA) is an extracellular tracer. We note however, that parenchymal GS transport of the extracellular tracer represents transport in ISF (volume fraction of ~18–22%) as well as the PVS spaces (volume fraction of ~1% or less). Alternative explanations would be less overall advection as the solute moves away from larger pulsatile pial vessels and/or increased diffusion mode of solute transport in deeper portions of parenchyma.

The OMT problem was first posed by Gaspard Monge in 1781⁵⁰ and may be formulated as one of transporting one distribution of mass to another in a manner that minimizes a given cost functional. Mass is used in a general sense and originally referred to soil, as Monge was interested in the civil engineering problem of leveling the ground. In our work on the GS, mass is represented by signal intensity induced by a paramagnetic contrast solute (reflecting Gd-DOTA concentrations²⁹) administered into CSF and captured by DCE-MRI. The model is related to the Navier-Stokes model but instead of considering momentum, one considers the flow of density^51,52. The advection/diffusion rOMT model is simpler to analyze, but nevertheless provides physically meaningful information. Interestingly, there are strong connections between rOMT and the porous media and kinetic models previously applied to characterize GS transport^30,40,53. However, in contrast to these models, the rOMT problem does not require a priori identification or delineation of boundaries, time-varying input or tissue function, or specification of kinematic parameter values. Other than the assumption of the underlying advection/diffusion equation, no prior assumptions are imposed on the GS flow. Finally, rOMT may be considered to be an optical flow technique⁵⁴, but instead of imposing smoothness on the solutions of the continuity equation, one employs smoothness as a constraint in minimizing a quadratic energy function.

GS dysfunction is reported in neurodegeneration particularly where vascular dysfunction may be involved^7,38,39. We therefore applied the rOMT framework to dissect modes of solute transport changes in an animal model of cSVD. cSVD is the major pathology in vascular dementia, and GS dysfunction has recently been highlighted as a novel pathogenetic factor in cSVD⁵. The hypothesized link between GS dysfunction and cSVD is based on the observation of abnormally enlarged PVS in clinical cases of cSVD⁵⁵. We utilized rOMT analysis to characterize modes of GS solute transport across tissue compartments in WKY and SHRSP rats. The SHRSP rat is a naturally bred animal with a parent control strain, which mimics human hypertension and cSVD⁴⁵. Another reason that we used the SHRSP rat is that they do not have hydrocephalus – a drawback inherent to other rat strains with chronic hypertension (e.g. SHR rat⁵⁶) which may affect CSF flow dynamics⁷. The rOMT analysis of DCE-MRI data from the brain showed that the total flux of moving solutes in all tissue compartments (including CSF) was reduced in SHRSP when compared to normotensive WKY rats.

Interestingly, mean solute speed in the CSF and tissue compartments was not significantly reduced in SHRSP when compared to WKY rats. However, a separate analysis of CSF transport in the PVS along the MCA obtained with DCE-MRI acquired with increased spatial resolution and higher contrast-to-noise ratio (increased Gd-DOTA in CSF) revealed altered CSF flow patterns in SHRSP when compared to WKY rats. In WKY rats, solute speed along the MCA was homogeneous whereas in SHRSP rats solute flow became overtly heterogeneous along the MCA. The rOMT analysis of the SHRSP MCA data revealed that solute speed was 3-fold higher at the MCA-PVS ‘root’ area when compared to speed distal to this area, suggesting a physical or physiological obstruction at this site.

We investigated possible mechanisms underlying the reduced CSF solute flux in specific regions such as the hippocampus of SHRSP compared to WKY rats. We documented altered perivascular AQP4 expression in SHRSP compared to WKY rats in hippocampus and occipital cortex which might have contributed to the reduced GS transport patterns including reduced PVS-to-tissue transfer of solutes. Reduced perivascular expression of AQP4 water channels has been documented in postmortem human brain with Alzheimer’s disease⁶ and in aging rodents³⁹.

In summary, the new rOMT Lagrangian framework analysis introduced here established the presence of diffusion and advective solute transport in brain parenchyma in part supporting the GS hypothesis and its relevance for waste clearance. In addition, altered solute speed profiles across tissue beds consistent with slower transfer of solutes from PVS to ISF was revealed in chronic hypertension. Future studies should be conducted to further optimize the rOMT Lagrangian analysis developed here, and to confirm the observed reductions in solute transport in human cSVD and other neurodegenerative states. In addition, future rOMT work will consist of including a source term in our set-up to include physical parameters (e.g., pressure) as well as combining OMT with other metrics, such as Fisher-Rao, in the cost functional to alleviate the explicit constraint of mass preservation in the model⁵⁷.

Methods

Animals

The local institutional animal care and use committees at Yale University, New Haven, USA approved all animal work. 8-month old WKY female rats and aged-matched female SHRSP rats (Charles River, Wilmington, MA, USA) were included in the study. Rats were given standard rat chow and water ad libitum. All methods were performed in accordance with the relevant guidelines and regulations.

Blood pressure measurements in SHRSP and WKY rats

Non-invasive blood pressure measurements were performed in a selection of awake WKY and SHRSP rats using tail cuff method (CODA monitor, Kent Scientific Corporation, USA). Three measurements were taken at 5-min intervals from each rat and the average of systolic and diastolic pressures were used for the analysis. In a separate series of anesthetized WKY (N = 6) and SHRSP rats (N = 5), mean arterial blood pressure was measured using a catheter inserted into the femoral artery. Anesthesia for these rats was performed with isoflurane and dexmedetomidine as described below.

Anesthesia and cisterna magna catheter

For all glymphatic experiments, rats were anesthetized with 3% isoflurane delivered into an induction chamber in 100% O₂ and after loss of the righting reflex received dexmedetomidine (0.015 mg/kg i.p.) and glycopyrrolate (0.2 mg/kg i.p.). Surgical anesthesia was maintained with 2–2.5% isoflurane and the rats were breathing spontaneously throughout the experiments. A small 5-mm copper tube (0.32 mm o.d., Nippon Tockushukan, MFG. CO., LTD, Tokyo, Japan) attached to a PE-10 microcatheter was implanted into the CSF via the cisterna magna and secured in place using cyanoacrylate glue. Following surgery, the rats were transferred to the MRI instrument. Anesthesia was maintained throughout the experiment by 0.5–1.0% isoflurane in 1:1 Air:O₂ mixture delivered via a nose cone, and a subcutaneous infusion of dexmedetomidine (0.015–0.020 mg/kg/hr)⁵⁸. During imaging, respiratory rate, oxygen saturation, body temperature and heart rate were continuously monitored using an MRI compatible monitoring system (SA Instruments, Stony Brook, NY, USA). Body temperature was kept within a range of 36.5 °C–37.5 °C using a heated waterbed.

MRI imaging

All MRI acquisitions were performed on a Bruker 9.4 T/16 magnet (Bruker BioSpin, Billerica MA) interfaced to an Avance III console controlled by Paravision 6.0 software. All rats were imaged in the supine position. For the whole-brain glymphatic experiments, we used a 3-cm planar receive surface radio-frequency (RF) coil (Bruker) which was positioned under the head of the rat with a separate preamplifier used as a receiver. For rats where glymphatic transport was explored along the middle cerebral artery, a 1-cm planar RF coil (Bruker) was placed above the left side of the rat’s head. A custom-made volume diameter volume coil (50-mm i.d.) was used as a transmit RF coil.

DCE-MRI for GS transport

Following an anatomical localizer scan, acquisition of the DCE-MRI started using a VFA-SPGR sequence. For whole brain glymphatics studies: TR (repetition time) = 25 msec, TE (echo time) = 4 msec, 0.30 × 0.30 × 0.30 mm FA (flip angle) = 15°, number of averages = 2, Time/scan = 5 min). For the MCA glymphatics studies (left hemisphere): TR = 25 msec, TE = 4 msec, 0.15 × 0.15 × 0.15 mm FA = 20°, number of averages = 1, scan time = 4 min 10 s). A reference phantom filled with 0.1 mM Gd-DOTA was placed in the field of view to allow both intensity normalization between each scan and receiver gain adjustment. For measurement of GS transport, DCE-MRI were acquired before (3 baseline scans), during, and after the infusion of gadoteric acid (Gd-DOTA, DOTAREM, Guerbert LLC, Carol Stream IL) into the CSF via the CM catheter³⁰. Specifically, for whole brain studies we infused Gd-DOTA (0.5 mmol/ml) in a 1:37 dilution: 1-part Gd-DOTA (0.5 mmol/ml):18-parts of sterile water:18-parts of 0.9% saline. This Gd-DOTA solution has a density of approximate 1.005 g/ml, which is isobaric with respect to CSF (CSF specific gravity = 1.004–1.007). The solution was infused at a constant rate of 1.5 µL/min and a total volume of 20 µl was infused. Intracranial pressure in the CM was measured in separate rats positioned in a stereotaxic frame under the same anesthesia regimen and with physiological parameters similar to those of rats undergoing glymphatic transport inside the MRI. The ICP was 4–5 mmHg prior to infusion and increased by ~1 mmHg during infusion. For the higher resolution MCA scans, we infused 20 µL, of 1:5 Gd-DOTA diluted with sterile water at a rate of 1.5 µL/min.

Morphometry for segmentation of tissue compartments and ROI

To study tissue compartment specific glymphatic transport, WKY (N = 7, 8-months old) and SHRSP (N = 7, 8-months old) rats were scanned using the 3D VFA-SPGR sequence (TR/TE/FA/number of averages = 15 ms/4 ms/5~30°/2 scanning time = 30 min 0.30 × 0.30 × 0.30 mm³) as described previously^4,56. Proton density weighted images calculated from VFA-SPGR were segmented using a custom made SHRSP-WKY rat brain template in deriving spatial registration parameters⁴. Spatially normalized segmented images were population averaged and thresholded (0.5) to create a population averaged binary mask for GM, WM and CSF, and inverse warped to match the individual rat brain. The GM, WM and CSF compartments were carefully segmented based on the pre-contrast images. CSF compartment contains PVS associated with pial vessels on the surface of the brain. However, the outer boundaries of the GM mask excluded pial vessels on the brain surface and therefore also PVS associated with pial vasculature. We also carefully avoided including large veins such as the strait sinus, transverse and sagittal sinuses. We note that the PVS associated with smaller vasculature inside GM and WM masks are unavoidably included, but as mentioned represents only a minor volume fraction (< 1–2% in normal brain) when compared to the ISF fluid space (18–22%). To study ROI specific glymphatic transport, the publicly available Waxholm Wistar rat brain atlas⁵⁹ was downloaded, and ROIs in the atlas were warped onto the individual SHRSP-WKY rat brains in two steps. First, the high resolution T2 weighted Waxholm Wistar rat brain atlas was segmented using the SHRSP-WKY template to derive the registration parameters between the atlas and the SHRSP-WKY template. The registration parameters were then applied onto the ROIs to construct the ROIs in the template space. ROIs in the template space were then inverse warped, which was derived from the tissue compartment analyses, to create the ROIs in the individual brain. Nineteen regions were included in the atlas: anterior commissure, basal forebrain region, brainstem, hippocampus, corpus callosum, descending corticofugal pathway, fimbria of hippocampus, hypothalamus, superior colliculus, cerebellum, cortex, olfactory bulb, penial, periaqueductal grey, septal region, striatum, subiculum, substantia nigra and thalamus. Only four of these (basal forebrain, cerebellum hippocampus and superior colliculus) were used for GS analysis purposes.

rOMT analysis of DCE-MRI data

In the next few sections, we provide the necessary details pertaining to our regularized OMT framework before describing how it was utilized for analyzing the DCE-MRI data.

Regularized OMT (eulerian framework)

The regularized OMT (rOMT) problem⁶⁰ considers the minimization of the energy functional

$$ {\mathcal R} [\mu ,v]=\,{\int }_{0}^{1}\mathop{\int }\limits_{{\mathfrak{X}}}\frac{1}{2}\mu (t,x)||v(t,x)|{|}^{2}dxdt$$

(1)

overall time-varying densities $\mu =\mu (t,x)\ge 0$ and sufficiently smooth velocity fields $v=v(t,x)\in {{\mathbb{R}}}^{n}\,$ subject to the advection/diffusion (constraint) equation

$${\partial }_{t}\mu +\nabla \cdot (\mu v)=\nabla \cdot (D\nabla \mu ),$$

(2)

for all $x\in {\mathfrak{X}}$, a connected Euclidean subdomain of ${{\mathbb{R}}}^{n}$, over the normalized time interval $t\in [0,1]$ with initial and final conditions

$$\mu (0,x)=\,{\mu }_{init}(x),\,\mu (1,x)={\mu }_{final}(x),\,\forall x\in {\mathfrak{X}}.$$

(3)

Here, ${\mu }_{init}$ and ${\mu }_{final}$ are given observed states of the density at respective times $t=0$ and $t=1$, and Ddenotes a symmetric positive definite diffusivity matrix, which may be spatially dependent or independent. When $D=0$, the constraint (2) is simply the advection equation and we recover the classical dynamic-OMT formulation³⁷.

From the Calculus of Variations^37,60, the Euler-Lagrange equations associated with the regularized OMT problem (Eqs. (1–3)) shows that the optimal velocity, ${v}_{opt}$, has the form

$${v}_{opt}=-\nabla \eta ,$$

(4)

where $\eta$ is the Lagrange multiplier of constraints (2 and 3) and satisfies the diffusion Hamilton-Jacobi equation,

$${\partial }_{t}\eta -\frac{1}{2}||\nabla \eta |{|}^{2}=\nabla \cdot D\nabla \eta .$$

(5)

The diffusion term plays a dual role. First of all, it regularizes the flow. In fact, it makes the problem equivalent to the so-called Schrödinger bridge problem, and therefore to a certain type of regularization (related to Fisher information) of the cost function^37,60. This leads to smoother flow fields. Secondly, it captures a key feature in modelling the underlying dynamics that may best be understood via compartmental models that will be described in the next section.

Regularized OMT and compartmental models

Compartmental models⁶¹ are comparable to spatially averaging the regularized OMT continuity Eq. (2). More precisely, one may average the continuity equation by integrating it over a control volume (i.e., compartment) C:

$$\mathop{\int }\limits_{C}{\partial }_{t}\mu \,dC+\mathop{\int }\limits_{C}\nabla \cdot (\mu v)\,dC=\mathop{\int }\limits_{C}\nabla \cdot (D\nabla \mu )\,dC.$$

(6)

Considering the advective term, we see that

$$\mathop{\int }\limits_{C}\nabla \cdot (\mu v)\,dC=\mathop{\int }\limits_{B}\mu v\cdot \overrightarrow{n}\,dB=-\,\mathop{\int }\limits_{B}\mu \nabla \eta \cdot \overrightarrow{n}\,dB$$

(7)

by (4) and the divergence theorem, where $\overrightarrow{n}$ is the unit outward normal of B, the compartment boundary. In this form, the OMT-advection may be interpreted as an expression of Darcy’s law, which relates flow velocity to the pressure gradient, as used in the porous media model. We note that the OMT-velocity does not explicitly account for pressure differences while the porous-media-velocity does. One can supplement a pressure dependent source term in the continuity equation, an extension we plan to explore in future work. Considering the unresolved true physical parameter values, this exemplifies the OMT framework’s utility for capturing, representing and informing on meaningful dynamic behavior implicitly encoded in the given observed density distributions, along with interesting insight into how the method works.

Applying the divergence theorem to the diffusive term we find:

$$\mathop{\int }\limits_{C}\nabla \cdot (D\nabla \mu )\,dC=\,\mathop{\int }\limits_{B}D\nabla \mu \cdot \overrightarrow{n}\,dB,$$

(8)

and thus, we see that diffusion is influenced by the density gradient, in accordance with Fick’s law. Noting that the total mass in C at time t is given by

$${\varrho }(t):=\mathop{\int }\limits_{C}\mu \,dC,$$

(9)

the averaged form of the diffusion-OMT constraint over the given compartment C relates the change in its density over time to the flux across its boundary, and (6) may be written as

$$\frac{d}{dt}{\varrho }(t)=\mathop{\int }\limits_{B}(D\nabla \mu -\mu v)\cdot \overrightarrow{n}.$$

(10)

Here, ${j}_{AD}=-(D\nabla \mu -\mu v)$ is the total flux, comprised of both advective and diffusive components. We note that the two are not independent of each other, as is seen by (5).

To illustrate the compartmental regularized OMT formulation, we compare tracer movement from CSF into the tissue compartment between normotensive (WKY) and hypertensive (SHRSP) cases. As such, we take the tissue compartment to be our control volume C and consider the surface boundary B between the CSF and tissue compartment (Supplementary Fig. 3a). For every boundary point $x\in B$, we can define the unit normal vector $\overrightarrow{n}(x)$ that is orthogonal to the surface at x and we orient the surface so that $\overrightarrow{n}$ points from the CSF compartment into the tissue compartment (Supplementary Fig. 3b). Using the OMT-derived velocity field and smooth images, we can then look at the total (advective + diffusive) flux ${j}_{AD}$ across this surface at each time step by considering the component of the flux in the normal direction. At a given time, (10) describes the net mass flowing through the surface. However, we are specifically interested in the net mass that moves from the CSF compartment into the tissue compartment. We are therefore interested in the normal component of the flux such that ${j}_{AD}\cdot \overrightarrow{n} > 0$. Due to the directional information provided by the rOMT, this can be readily be determined. Fluctuations in voxel-level information given by rOMT are clearly lost when averaged, indicative of the well-known challenge for the kinetic analysis approach for local insight.

Regularized OMT-lagrangian framework

The Eulerian framework records measurements at fixed locations for multiple points in time, informing on the collective instantaneous behavior of the flow. We now describe our Lagrangian framework for expressing the trajectories of fixed parcels or particles. In this work we have assumed that the diffusion is constant, $D={\sigma }^{2}I$. While this assumption is somewhat restrictive, it nevertheless gives physically meaningful results. All the theory goes through in the non-constant spatially dependent setting. In future work, we plan to examine this extension in more detail. We note that

$$\nabla \cdot (\mu v)-\nabla \cdot ({\sigma }^{2}\nabla \mu )=\nabla \cdot (\mu v)-{\sigma }^{2}\nabla \cdot (\mu \nabla \,\log \,\mu )$$

(11)

$$=\nabla \cdot [\mu (v-{\sigma }^{2}\nabla \,\log \,\mu )].$$

(12)

Defining the augmented velocity

$${v}^{aug}\,:\,=v-{\sigma }^{2}\nabla \,\log \,\mu ,$$

(13)

we derive the conservation form of the regularized constraint (2),

$${\partial }_{t}\mu +\nabla \cdot (\mu {v}^{aug})=0.$$

(14)

The Lagrangian coordinates $L=L(t,x)$ of the optimal trajectory for the regularized OMT system are given by

$${\partial }_{t}L={v}_{opt}^{aug}(t,L(t,x),\,L(0,x)=x,$$

(15)

where

$${v}_{opt}^{aug}={v}_{opt}-{\sigma }^{2}\nabla \,\log \,{\mu }_{opt},$$

(16)

and ${v}_{opt}$ and ${\mu }_{opt}$ denote the optimal arguments of the regularized OMT problem (1–3).

Numerical regularized OMT model for glymphatic system

In this section, we discuss the considerations taken for applying the regularized OMT model to the glymphatic system and its numerical implementation. For our purposes, MR signal enhancement induced by Gd-DOTA is treated as density (the relationship between mass and signal enhancement has been previously validated⁶²). The images then provide discrete observations ${\mu }_{init}$ and ${\mu }_{final}$ of the tracer distribution at some initial time $t=0$ and a final time $t=1$, respectively. Space is discretized as a cell-centered grid of size ${m}_{1}\times {m}_{2}\times {m}_{3}$ with a total number of $M$ cells, each with uniform length $dx$, as delineated by the voxels. Time is split into $n$ intervals of length $dt$ and the $N+1$ time steps are denoted by the subscript $n$ where $n=0$ corresponds to the time $t=0$, $n=N$ corresponds to the time $t=1$ and $dt\ast N=1$. Bold font is used to denote linearized variables and we use ${\boldsymbol{\mu }}={[{{\boldsymbol{\mu }}}_{0}^{T},\ldots ,{{\boldsymbol{\mu }}}_{N}^{T}]}^{T},\,{\boldsymbol{v}}={[{{\boldsymbol{v}}}_{1}^{T},\ldots ,{{\boldsymbol{v}}}_{N}^{T}]}^{T}$ to represent the temporal concatenation of the linearized densities and velocities, respectively. Note that v_n is the velocity field that describes the evolution between ${\mu }_{n-1}$ and ${\mu }_{n}$ subject to the advection/diffusion Eq. (2), for $n=1,\ldots ,N$.

Following the work of Steklova and Haber⁶³, the numerical method used here was chosen for its nice properties for differentiating and optimizing the resulting discrete regularized OMT problem, which will become clear in the discussion to follow.

We start by discretizing the advection/diffusion Eq. (2) as

$$({{\boldsymbol{I}}}_{M}-dt{\boldsymbol{D}}){{\boldsymbol{\mu }}}_{n+1}={\boldsymbol{A}}({{\boldsymbol{v}}}_{n}){{\boldsymbol{\mu }}}_{n},\,n=0,\ldots ,N-1,$$

(17)

where ${{\boldsymbol{I}}}_{s}$ is the s × s identity matrix, ${\boldsymbol{D}}$ is the discrete diffusion operator $\nabla \cdot D\nabla$ and ${\boldsymbol{A}}$ is an averaging matrix that assigns the density of the advected particle in space to the cell centers. This is derived via operator splitting where the advective step is performed using a particle-in-cell method and the backward Euler method is used for the diffusive step, see Steklova and Haber⁶³ for more details. Notice, the density for any time step now depends only on the velocity ${\boldsymbol{v}}$ and initial density ${{\boldsymbol{\mu }}}_{0}$. Considering that the initial density is fixed by the given image ${{\boldsymbol{\mu }}}_{0}={{\boldsymbol{\mu }}}_{init}$, we have reduced the parameter space of our objective cost functional (1), so that we only need find the optimal velocity.

A straightforward discretization of the OMT energy functional ${\mathcal R} [v]$ (1) (as just noted we ignore the density term) is given as

$${\boldsymbol{ {\mathcal R} }}[{\boldsymbol{v}}]\,:=dt{(dx)}^{3}{{\boldsymbol{\mu }}}^{T}({{\boldsymbol{I}}}_{N}\otimes [{{\boldsymbol{I}}}_{M}|{{\boldsymbol{I}}}_{M}|{{\boldsymbol{I}}}_{M}])({\boldsymbol{v}}\odot {\boldsymbol{v}}),$$

(18)

where $\otimes$ denotes the Kronecker product, $[\cdot |\cdot ]$ denotes block matrices and $\odot$ denotes the Hadamard product. We seek the velocity ${\boldsymbol{v}}$ that characterizes the density ${\boldsymbol{\mu }}$ which interpolates between the given initial and final densities, ${{\boldsymbol{\mu }}}_{init}$ and ${{\boldsymbol{\mu }}}_{final}$, seen to evolve according to the advection/diffusion Eq. (17) with minimal energy, as described by the discrete objection function ${\boldsymbol{ {\mathcal R} }}[{\boldsymbol{v}}]$ (18). Lastly, in anticipation of noise in the signal intensity, we do not attempt to match the final density exactly (3). Instead, we explicitly regard the final image as a noisy state of the tracer distribution

$${\mu }_{N}+\psi ={\mu }_{final}$$

(19)

where $\psi$ is independent, identically distributed Gaussian noise with covariance $\Psi$, treated as a free temporal endpoint of the partial differential equation. As such, the discrete regularized OMT problem is given as

$${\rm{\inf }}\,{\boldsymbol{ {\mathcal R} }}[{\boldsymbol{v}}]+\xi ||{{\boldsymbol{\mu }}}_{N}-{{\boldsymbol{\mu }}}_{final}|{|}_{\Psi }^{2}$$

(20)

such that

$$\{({{\boldsymbol{I}}}_{M}-dt{\boldsymbol{D}}){{\boldsymbol{\mu }}}_{n+1}={\boldsymbol{A}}({{\boldsymbol{v}}}_{n}){{\boldsymbol{\mu }}}_{n},\,n=0,\ldots ,N-1$$

(21)

$$\{{{\boldsymbol{\mu }}}_{0}={{\boldsymbol{\mu }}}_{init}.$$

(22)

$$\{{\boldsymbol{\mu }}\ge 0$$

(23)

Here, the free-endpoint condition added to the objective function (20) inhibits the erroneous influence of noise on the velocity, inversely weighted by the covariance and the parameter $\xi$ indicates the relative significance of fitting the data to minimizing the transport energy.

Using no diffusion or too small of a diffusion value yields implausible results (e.g. non-smooth velocity fields due to overcompensation by advection which therefore cannot reflect proper trajectories), but too large a value will have an accumulated effect that may end up over-smoothing the images over time (e.g. intensity gradients will be blurred yielding fictitious advective behavior attempting to reverse the spreading and subtler information will be lost). Multiple levels of diffusion were tested to determine values that afforded an appropriate balance. We note that for small changes in diffusion values (i.e. changes of the same order of magnitude), we found the model to be robust to the diffusivity. The aforementioned implausible behavior was only observed when the magnitude of the diffusivity was changed by an order of 2 or more, as expected. One would expect that increasing diffusion in a system of diffuse particles would allow for more random or unexpected behavior to occur. In this manner, diffusion allows for a natural description of stochastic particle trajectories thereby regularizing the flow. Accordingly, adding the diffusion term to the classical OMT constraint enables our model to account for more behavior than would be possible with advection alone. We therefore refer to our system of Eqs. (20–23) as the free regularized OMT problem (frOMT), used interchangeably with rOMT), which turns out to be a quadratic optimization problem with linear constraints with respect to ${\boldsymbol{v}}$. The optimal velocity and interpolated densities are then found using a Gauss-Newton method. See Table 1 for a list of symbols and their interpretations.

Table 1 Symbol reference list.

rOMT data analysis

Image processing of the series of DCE-MRI data of each rat consisted of head motion correction, intensity normalization, smoothing, and voxel-by-voxel conversion to percentage of baseline signal using SPM12^64,65 (https://www.fil.ion.ucl.ac.uk/spm/). The processed images were converted into NIfTI format for further analysis. DCE-MRI images were fed as inputs into our rOMT model in order to capture advective and diffusive transport in different compartments including CSF spaces, grey matter (GM), white matter (WM), brain regions (and PVS) along the MCA. Unless otherwise specified, all computation was carried out in MATLAB (Version 2016b). The rOMT procedure was run in parallel on all 14 whole brain datasets (8 WKY, 6 SHRSP) on the SeaWulf (Stony Brook University) CPU cluster (with a total of 4592 nodes and 128 gigabytes RAM per node) to model a 120 min time period of tracer transport captured across 21 images, which took approximately 52 hours to execute. We note that fewer interpolated time steps could have been used to achieve a shorter runtime, but we felt a thorough investigation warranted the extended time. rOMT was performed on the 8 higher resolution MCA datasets (4 WKY and 4 SHRSP) in the same manner over a time period of 144 min captured across 25 images with a total runtime of approximately 66 hrs. More time was required due to the longer time period being considered and decreased signal-to-noise ratio. Again, fewer intermediate steps could be used with minimal differences in output for future applications.

We used the Lagrangian OMT formulation (vide supra) to construct the pathlines for visualizing GS transport flows over the 120 min time period in one comprehensive figure, derived from the rOMT returned velocity field and interpolated images. In order to extract and visualize trajectories from 4D data in a meaningful way, the several steps were performed on each dataset. First, pathline start points were selected uniformly among voxels that showed at least 12% increase in signal as compared to baseline throughout the imaging acquisition interval. This was done in order to reduce the complexity of the visualization, as the trajectory of a parcel from a location with such low fluctuation may be attributed to noise or insignificant behavior. Pathlines were then parametrized by the Lagrangian coordinates (21) associated with each start point and were computed by integrating the augmented time-varying velocity field (19), as described above. Particle attributes such as speed, informing on the state of the particle along its trajectory, were simultaneously recorded and referred to as speed-lines. The pathlines were subsequently clustered by spatial proximity using the QuickBundles algorithm⁶⁶ in Python (Version 3.5.4), and pathlines of small clusters deemed ‘insignificant’, as determined by a predefined threshold, were removed from consideration. ‘Significant’ pathlines were then converted back from spatial coordinates to the cell centered grid using a tri-linear averaging matrix and saved in nifti format for visualization in Amira. Runtime for the Lagrangian procedure is ~3 min per dataset on a standard laptop with a 2.2 GHz Intel Core i7 processor.

Time-varying particle attributes associated with the pathlines, such as speed and density, provide more detailed information regarding the particle dynamics. We selectively focused on the speed for the novel insight it provides. Specifically, speeds that make up the aforementioned speedlines were computed as the magnitudes of the rOMT-derived transport vectors ascribed to the time and place of each point along the particle’s path and then visualized to analyze speed trajectories in different brain regions and across normotensive and hypertensive states.

For each rat, we computed one Lagrangian rOMT image reflecting pathlines and one reflecting pathline speed (i.e. speedlines), at a voxel resolution of 0.3 × 0.3 × 0.3 for whole brain and 0.15 × 0.15 × 0.15 mm³ for hemispheric data. For the whole brain data, total pathline volume and mean pathline speed were extracted from CSF, GM, WM and four predetermined brain regions (basal forebrain, superior colliculus, hippocampus and cerebellum). For MCA data, pathline volume and mean pathline speed were extracted in regions along the MCA segment on the base of the brain (a.k.a. MCA ‘root’) and along the left MCA.

Immunohistochemistry

Animals were deeply anesthetized with ketamine/dexmedetomidine (150 mg kg⁻¹/2 mg kg⁻¹ i.p) and were transcardially perfused with phosphate buffered saline (PBS) followed by paraformaldehyde (PFA; 4% in PBS). The brains were extracted, post-fixed overnight in PFA at 4 °C and transferred to PBS. All incubations were carried out at room temperature unless specified. Brains were sectioned coronally (50 µm) with a vibratome (Leica VT1000 S) and underwent heat-induced epitope retrieval in citrate buffer (10 mM citric acid, 0.05% Tween 20, pH 6.0) at 85 °C for 30 min. The sections were then blocked and permeabilized for 1.5 h (1% BSA, 1% normal goat serum, and 0.3% Triton) and incubated for 2 h at 25 °C with primary antibodies: rabbit anti-AQP4 (Novus, NBP1–87679, 1:500) and chicken anti-GFAP (Abcam, ab4674, 1:1000) in PBS, 0.05% Tween 20, and 1% BSA. Following this, sections were incubated for 2 h at 25 °C with Alexa Fluor-Plus goat anti-rabbit 555 and Alexa Fluor goat anti-chicken 647 secondary antibodies (Invitrogen, 1:1000) in PBS with 0.05% Tween 20. The sections were then mounted on slides and cover slipped.

The MCA brains were also processed for GFAP, AQP4 and collagen IV. The PFA fixed brains were cryoprotected in 30% (w/v) sucrose in PBS, embedded in O.C.T. compound and snap-frozen at −80 °C. Sections were cut coronally (20 µm) using a cryostat, mounted on Colorfrost/Plus slides and stored at −80 °C. The sections were rehydrated in PBS, incubated with proteinase K for 15 min (0.2 mg/ml) and incubated in blocking buffer for 30 min (0.3% Triton in SuperBlock buffer). The sections were then incubated overnight at 4 °C with primary antibodies: goat anti-collagen IV (Novus, NBP-126549,1:200), rabbit anti-GFAP (Dako, Z0334, 1:250), rabbit anti-AQP4 (Novus, NBP-1-87679) and goat anti-Iba-1 (Novus, NBP-100-1028, 1:250). The next day, sections were incubated with Alexa Fluor 594 or 488 secondary antibodies (Invitrogen, 1:1000). The sections were then cover slipped and sealed. Images were captured using a Keyence fluorescence microscope BZ-X700 system.

Quantification of AQP4 polarization

Images of AQP4 immunofluorescence were captured on a Keyence BZ-X700 automated fluorescence microscope using 4–20 × air objectives as z-stack tile-scans with 3 µm z-steps in maximum contrast projections. Quantification of AQP4 polarization was done in the ventral hippocampus and occipital cortex using ImageJ software (Image J 1.52i). The AQP4 digitized images were imported into ImageJ, scaled and converted to optical density images. In order to unbias the analysis, a grid (150 numbered boxes) was overlaid over the ventral hippocampus (each box area = 228 × 228 pixels) and occipital cortex (each box area = 282 × 282 pixels). 15 boxes from hippocampal or cortical sections were randomly selected for a total of ~30 small vessels and ~30 capillaries for WKY rats and the same number of vessels in SHRSP rats. Representative line segments (300 µm for small vessels and 100 µm for capillaries) were used to extract immuno-intensity across the micro-vessels selected in each grid-box. The polarization index for each vessel was calculated (peak-to-baseline level).

Quantification of GFAP and AQP4 mRNA by real-time quantitative RT-PCR (qPCR)

Astrogliosis and AQP4 water channel expression were evaluated in the brains of WKY and SHRSP rats suing real-time quantitative RT-PCR (qPCR). Specifically, for this purpose, mRNA levels of the two targets (GFAP and AQP4) by qPCR was performed with the QuantStudio 3 Real-Time PCR system (Applied Biosystems). Brain tissues of interest (WKY, N = 8; SHRSP, N = 8) were quickly dissected, snap frozen in liquid nitrogen and stored at −80 °C. Total RNA extraction was carried out using the PARIS kit (Invitrogen) following brief homogenization using the TissueRuptor II (Qiagen). Total RNA concentrations were measured using the SpectraMax i3x Multi-Mode microplate reader (Molecular Devices). Samples were then treated with ezDNase enzyme (Invitrogen) to degrade any genomic DNA and subsequently reverse transcribed to cDNA using the Superscipt IV VILO kit (Invitrogen). Quantification of target mRNA was performed using TaqMan gene expression assays (ThermoFisher Scientific). The PCR cycle parameters were set according to the recommended guidelines for the TaqMan Fast Advanced Master Mix (Applied Biosystems), which was used to run the qPCR reactions (95 °C for 2 min [polymerase activation] followed by 40 cycles of 95 °C for 1 s [denture] and 60 °C for 20 s [anneal/extend]). TaqMan primer/probe sets for AQP4 (Rn00563196_m1) and GFAP (Rn00566603_m1) were used. mRNA target signals were normalized to the housekeeping gene RPLP0 (Rn00821065_g1). Relative fold change values were calculated using the 2^−ΔΔCt method⁶⁷. The mRNA target signals were first normalized to the reference housekeeping gene RPLP0. Following this, relative fold change values were calculated using the 2^−ΔΔCt method by normalizing the data to the cerebellum, which was used a reference region.

Statistical information

Group sizes for the whole brain glymphatics studies were estimated based on previous studies^40,58. For peak magnitude point comparisons, group sample sizes of 6 (WKY) and 6 (SHRSP) achieve 80% power to detect a Cohen’s effect size of 2.0, and with a significance level (α) of 0.025 using a two-sided two-sample unequal-variance t-test. All data are presented as mean ± SD. If data showed normal distribution, a paired or unpaired Student’s t-test was used to compare groups. When data lacked normal distribution, groups were compared using non-parametric testing (Mann-Whitney test for unpaired comparison, Wilcoxon signed rank test for paired comparisons). For within-group regional differences Friedman’s ANOVA was used. All testing was two-tailed and exact P values were calculated at a 0.05 level of significance and stated in the figure legends. Blinding for data analysis was done where possible. All statistical analyses were performed using XLSTAT statistical and data analysis solution (Addinsoft (2019). XLSTAT statistical and data analysis solution. Boston, USA. https://www.xlstat.com) or GraphPad Prism 7 (GraphPad Software).