Arterial pulsations drive oscillatory flow of CSF but not directional pumping

The brain lacks a traditional lymphatic system for metabolite clearance. The existence of a “glymphatic system” where metabolites are removed from the brain’s extracellular space by convective exchange between interstitial fluid (ISF) and cerebrospinal fluid (CSF) along the paravascular spaces (PVS) around cerebral blood vessels has been controversial. While recent work has shown clear evidence of directional flow of CSF in the PVS in anesthetized mice, the driving force for the observed fluid flow remains elusive. The heartbeat-driven peristaltic pulsation of arteries has been proposed as a probable driver of directed CSF flow. In this study, we use rigorous fluid dynamic simulations to provide a physical interpretation for peristaltic pumping of fluids. Our simulations match the experimental results and show that arterial pulsations only drive oscillatory motion of CSF in the PVS. The observed directional CSF flow can be explained by naturally occurring and/or experimenter-generated pressure differences.

The velocity v is the time derivative of the displacement. The following relation holds for the Eulerian description of displacement v(x, t) : where ∇ x is the gradient with respect to x For calculations in ALE, we define the mesh configuration, where the points are denoted by X m . The deformed configuration B t is related to the mesh configuration B m by the map χ m . The displacement of the mesh domain (mesh motion for finite element calculations) is given by where, ∇ Xm is the gradient with respect to X m over B m . The primary unknown fields,for example, the stress tensor(σ) fluid velocity (v) and pressure (p) in the mesh coordinates X m are given by:p The material time derivative of x(X m , t) is the material particle velocity in mesh Coordinates, which should be the same as the material particle velocity defined in eq. A.3 The deformation gradient defined in eq. A.5 and the Jacobian determinant defined in eq. A.6 are used for transforming the spatial derivatives to the mesh domain.
For a vector v: where, Tr is the trace operator on tensors.
For a scalar p: For a tensor σ The material time derivative of velocity(the acceleration) is calculated aŝ A detailed derivation of eq. A.15 was done by Donea et al. [5].
Integrals over the deformed domain can be transformed to the mesh domain as follows where ν and ν m are infinitesimal volumes is B t and B m respectively. Integrals over the boundaries are transformed as follows where a and a m are infinitesimal areas in ∂B t and ∂B m respectively and n m is a unit normal in the mesh domain.

B Anisotropic non-dimensionalization
We transform the mesh domain into the non-dimensionalized computational domain with coordinatesX c . The transformation is defined by the characteristic length L o and scaling tensor G. For our problem, the characteristic length and the scaling tensor are constants.
Any field in X m can be transformed into a field in X c using the relation in B.1 The spatial derivatives are transformed as follows.
For a vector v: For a scalar p: For a tensor σ: The time derivatives are unaltered because there is no relative motion between X m and X c .
The volume and surface integrals are transformed as follows.
Integrals over the boundaries are transformed as:

C Initial conditions for problems
In all the intial-boundary value problems used in the study, the initial values and initial time-derivatives of all variables are set to zero because the initial conditions of the system are unknown. The parameter values used in the boundary conditions, are ramped up from zero to the specified values using step functions available in Comsol. These functions have continuous first and second derivatives with respect to time. In this document we will use step(t) in the equations to indicate wherever these functions are used. The results for these simulations are shown after 4 cycles of the periodic wall motion, where the variations in the velocity and pressure values from cycle to cycle are less than 10 −4 % of peak value.

D Navier-Stokes in ALE coordinates
We solve for the fluid velocity v f and pressure p f in a deforming domain B t , representing the PVS. The time-dependent deformation of the domain is represented by u m (same as eq. A.4).We present the strong form equations here in the mesh coordinates X m .The equations can be verified with previously published literature [3,2].The appropriate transformations for change of coordinates are applied (eq. A.11 -eq. A.17) and the X m subscript on the spatial derivatives (gradient, divergence and laplacian) is dropped.
We use a linear elliptical model for the mesh motion [1,4].The governing equation for the mesh displacement, u m on B m is given by The equations in the mesh coordinates are written for the velocity and pressure in the mesh coordinates(v f andp f ). The governing equation for the fluid in the PVS(Navier-Stokes) are: The boundary of the computational domain is divided into three subsets, Γ D , Γ P and Γ N , representing the surfaces with Dirichlet, Periodic and Neumann Boundary conditions respectively. The periodic boundary, Γ P , is further divided into Γ P 1 and Γ P 2 , representing the source and destination for the periodic boundary conditions. We define the following functional spaces in the computational domain (B c ) for velocity, pressure and mesh displacement respectively.
Additionally, we define the spaces In equations E.1 to E.6, d is the number of dimensions of the problem. The abstract weak formulation of the problem is as follows: For cases where there are surfaces with a non-zero traction boundary conditions(Γ D1 ), an additional term T BC f (eq. E.11) is added to left side of equation E.7.

F Particle tracking in ALE
To find the fluid particle trajectories in the deforming domain, B t , we first solve for the fluid particle coordinates in the mesh domain (X m ) as a function of time. The material time derivative of X m is the particle velocity observed from the mesh domain B m is defined as: We need to calculate the material time derivative of the mesh displacement. We use a method similar to that in eq. A.3, where we use the map χ −1 m oχ p instead of the map χ p in eq.A.3.
The finite element computations in the computational coordinates give us the fieldsv(X c , t) and u m (X c , t). These fields can be transformed into mesh coordinates using the inverse of the transformation in eq. B.2. We calculateẊ m as a function of these known quantities using equations A.10, F.1 and F.2 in the material time derivative of eq. A. 4.v Using the definition of F m from eq. A.5 we obtain the particle velocity in the mesh domaiṅ We use the following initial value problem to calculate fluid particle trajectories the deforming domain. Here, it is assumed that the material particle velocityv(X m , t) and the mesh displacement u m (X m , t) (and subsequently F m and ∂um ∂t ) are known fields in the mesh domain. We first find the particle trajectory in the mesh coordinates and add the mesh displacement to find the particle trajectory in the deformed coordinates.
f ind x(X m , t) such that x(X m , t) = X m + u m (X m , t) With the initial condition: x(X m , 0) = X m + u m (X m , 0) We solve this problem using a forward Euler integration scheme and calculate the fluid particle trajectories. When the position in the mesh domain at time t is known, the time integration scheme allows us to find the position at time t+dt.
The fluid particles can only be tracked till the time they exit the computational domain since the velocity field outside is unknown.