Functional hyperemia drives fluid exchange in the paravascular space

The brain lacks a conventional lymphatic system to remove metabolic waste. It has been proposed that directional fluid movement through the arteriolar paravascular space (PVS) promotes metabolite clearance. We performed simulations to examine if arteriolar pulsations and dilations can drive directional CSF flow in the PVS and found that arteriolar wall movements do not drive directional CSF flow. We propose an alternative method of metabolite clearance from the PVS, namely fluid exchange between the PVS and the subarachnoid space (SAS). In simulations with compliant brain tissue, arteriolar pulsations did not drive appreciable fluid exchange between the PVS and the SAS. However, when the arteriole dilated, as seen during functional hyperemia, there was a marked exchange of fluid. Simulations suggest that functional hyperemia may serve to increase metabolite clearance from the PVS. We measured blood vessels and brain tissue displacement simultaneously in awake, head-fixed mice using two-photon microscopy. These measurements showed that brain deforms in response to pressure changes in PVS, consistent with our simulations. Our results show that the deformability of the brain tissue needs to be accounted for when studying fluid flow and metabolite transport.


A Kinematics in ALE
To describe the motion of a continuum(solid, fluid), we start with two configurations, the deformed configuration B t , the reference configuration B p . As shown in Fig 1, the deformed configuration B t is the image of B p under the map χ p . The points in B p , B t are denoted by X p , x respectively. The displacement of material particles from the reference to the deformed configuration is given by.
Alternatively, the Eulerian description of the displacement is given by Where X p (x, t) is the image of x under the map χ −1 p at time t. Denoting the material time derivative by D t , the material particle velocity v is related to the Eulerian description of the displacement field as follows v(x, t) : where ∇ x is the gradient with respect to x over B t . The above equation for the material time derivative holds true for the Eulerian description of any vector quantity. It can be viewed as a property of the map χ p Additionally, for the ALE computation, we have the computational configuration B c . The points in B c are denoted by X c . As shown in Fig. A.1, the deformed configuration B t is the image of B c under the map χ c . In a finite element computation, the map χ c can be interpreted as the motion of the mesh. The displacement of the mesh gives us the relation between position of a particle in the computational domain and its position in the deformed configuration.
The deformation gradient and the Jacobian determinant of the mesh motion are given by where, ∇ Xc is the gradient with respect to X c over B c . Now, we take the material time derivative of all the quantities in eq. A.4. The material time derivative of x(X c , t) is the ALE representation of material particle velocity.
We want to find the solution to an initial boundary value problem(or a boundary value problem) on the deforming domain. However, since our finite element computation is performed on the computational domain B c instead of the deforming domain B t , the primary unknown fields,for example, the stress tensor(σ) fluid Figure A.1: Current configuration B t is the common image of particle motion χ p and mesh motion χ c velocity (v) and pressure (p) will be calculated as functions of the computational coordinates X c instead of the physical coordinates x, and time t.
We will need to express the governing equations in continuum mechanics, that contain spatial gradients and material time derivatives, in ALE coordinates. The following transformations use the deformation gradient defined in eq. A.5 and the Jacobian determinant defined in eq. A.6 to rewrite the spatial derivatives(that commonly occur in mass and momentum balance laws) in the computational 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) estimated in the physical and computational coordinates is given as follows.
The derivation of eq. A.16 is not immediately obvious. A detailed derivation is given by Donea et al. [3]. Further, any integral over the deformed domain can be transformed to an integral over the computational domain using the Jacobian determinant defined in eq. A.6.
where ν andν are infinitesimal volumes is B t and B c respectively. An integral over the boundary of the deformed domain can be transformed into an integral over the boundary of the computational domain as follows where a andâ are infinitesimal areas in ∂B t and ∂B c respectively andn is a unit normal in the computational domain.

B Particle tracking in ALE
To track fluid particles in the computational domain(in ALE), we want to find the fluid particle coordinates (X c ) as a function of time. The material time derivative of X c is the particle velocity observed from the computational domain B c is defined as: The material time derivative of the mesh displacement is obtained by a method similar to that in eq. A.3. In this case we use the map χ −1 c oχ p instead of the map χ p in eq.A.3.
From the finite element computation in ALE, we know the fieldsv(X c , t) and u m (X c , t). We want to findẊ c as a function of these known quantities. Using eq.s A.7, B.1 and B.2 in the material time derivative of eq. A.4 and rearranging the terms, we get We can use the definition of F m from eq. A.5 to obtain the particle velocity as seen from the computational domainẊ Here, we present the initial value problem to calculate fluid particle trajectories in ALE framework. This calculation is proposed as a post-processing step after the actual finite element computation. Therefore, it is assumed that the material particle velocityv(X c , t) and the mesh displacement u m (X c , t) (and subsequently F m and ∂um ∂t ) are known quantities in the computational domain.
The above initial value problem can be solved using a forward Euler integration scheme to obtain fluid particle position in the computational and deformed domains. The time integration scheme is implemented as follows to calculate the position in the computational domain at time t+dt when the position at time t is known.
The particle position in the deformed domain is simply a consequence of the first equation of the initial value problem(eq. B.6). We can track fluid particles till the time they exit the computational domain.
There are two main reasons for choosing a forward Euler integration scheme over a backward Euler integration scheme. Firstly, the ODE described in eq. B.6 does not necessarily have a fixed point. Secondly, the values of the quantities on the right side of the ODE only exist within the computational domain, and therefore the equation is not solvable using a backward Euler scheme at the timestep when a particle leaves the computational domain.

C Initial conditions for problems
In all the intial-boundary value problems described below, the initial values and initial time-derivatives of all variables are always set to zero. This is chosen for convenience since specifying a non-zero initial condition for one variable would require the knowledge the corresponding initial values for all the other variables with coupled physics. The parameter values used in Dirichlet boundary conditions, are ramped up from zero to the specified values using step functions in Comsol multiphysics. These functions have continuous first and second derivatives with respect to time. In the subsequent sections, we will use step1(t) in the equations to indicate wherever these functions are used. The results for these simulations are shown after the effect of the initial conditions are insignificant. We show the results after 20 cycles of the peristaltic wave, where the variations in the velocity and pressure fields from cycle to cycle are less than 10 −6 % of peak value.

D Darcy-Brinkman Flow in ALE coordinates
We want to solve for the fluid velocity v f and pressure p f in a deforming domain B t , representing the PVS. The displacement field that defines the time-dependent deformation of the domain is denoted by u m (same as eq. A.4). The finite element calculations are done in the computational domain B c , where we solve for the velocity and pressure fields as a function of the computational coordinates X c . These fieldsv f and p f are defined according to equations A.9 and A.8 respectively.
The computations are done in an axisymmetric framework. The inner and outer radius of the computational domain are R 1 & R 2 respectively, where R 2 = R 1 + wd. The domain has a length L a in the z direction. See Table 1

for the list of parameters
We use a harmonic model for the mesh motion [1,2].The governing equation for the mesh displacement, u m on B c is given by From here on in the document, since all the equations are presented in the computational coordinates, the X c subscript on the spatial derivatives (gradient, divergence and laplacian) is dropped.

The boundary of B c (∂B c =: Γ c )is divided into four non-overlapping regions such that Γ
where, ∅ is the empty set. The division of the boundary is shown in Fig D.1. The mesh displacement on the boundaries should reflect the motion of the physical domain. Here, the outer wall of the artery(Γ D1 c ) moves according to the sinusoidal wave described by eqn D.2.
on Γ D1 c : on Γ D2 c : By converting all the derivatives to the computational coordinates(eqs A.11-A.18), we get the ALE formulation.
Whereâ f is the acceleration of fluid particles defined according to eq. A.16 and I is the identity tensor. The body force, b is zero for our problem.
We use a no-slip boundary condition at the arterial wall(Γ D1 c ) and the brain tissue(Γ D2 c ) boundaries of the PVS. At the axial ends of the PVS, we assume that there is no applied traction.
The weak form of equations D.1, D.8 and D.10 are solved with the boundary conditions in equations D.2-D.4, D.11 and D.12. We discretize the 2D(axisymmetric) geometry using Lagrange polynomials of second order for u m andv f , and Lagrange polynomials of first order forp f . We use a fully coupled time dependent solver and a backward difference time integration scheme with a time step of 0.001s.

E Fluid Structure Interaction
We have two domains B f and B s , representing the PVS and the Brain tissue respectively. We want to solve for the velocity(v f ) and pressure(p f ) fields in B f . We continue to use the hat notation for these fields because they are calculated in the computational domain that represents the un-deformed PVS. We use the Darcy-Brinkman flow model for the fluid dynamics. The displacement of the computational domain is given the field u m . In B s , we want to calculate the solid displacement u s and solid velocity v s . Saint-Venant-Kirchoff model is used for the solid elasticity.
For this problem, we usen f andn s to represent the outward normals to the boundaries of B f and B f respectively.
The function b(t) for vasodilation starts with zero initial value and is shown in main Fig 4. on Here, Q 1 is the flow rate out of the PVS through the face Γ N 1 f . R s is the resistance of the subarachnoid space, which is taken to be 1% the flow resistance of the PVS. Equation E.6 is a Robin boundary condition, where the traction at a surface is proportional to the flow rate out of the surface. It serves as a lumped model for the subarachnoid space.
For the domain B s , the mesh motion is given by the solid displacement u s . This is the Lagrangian framework, used commonly in solid mechanics. The deformation gradient F s and the Jacobian determinant J s are defined similar to those in eqs A.5 and A.6. The strain energy for the Saint-venant-Kirchoff elastic model and the first Piola-Kirchoff stress are given by here, T r is the trace operator in tensors.
Since the mesh in the domain B s follows the material particles, the material time derivatives in B s are equal to the partial derivatives. The governing equations for u s , v s are given by, For our problem, the body forces b s are zero in B s .
The boundary conditions on Γ s − Γ f s are as follows.
on Γ D s : At the solid fluid interface, we enforce continuity of velocity and traction.
The weak form of equations D.1, D.8, D.10 on B f and E.11 and E.12 on B s are solved with the boundary conditions in equations E.1-E.7 and E.13-E.17. We discretize the 2D(axisymmetric) geometry using Lagrange polynomials of second order for u m ,v f , u s and v s and Lagrange polynomials of first order forp f , p s , q s and j s . We use a fully coupled time dependent solver and a backward difference time integration scheme with a time step of 0.0001s for arterial pulsations and 0.0001s for functional hyperemia. Since only the fluid leaves the PVS, the change in the fluid volume of the PVS is equal to the change in the total volume of the PVS. The calculation of the initial and deformed volume of the PVS is straight forward from A.17. The change in the PVS volume, δV t at any time t can be calculated over the undeformed configuration of the fluid domain B f .
From these two quantities, we can calculate the fraction of PVS fluid volume exchanged with the SAS at any time t. The volume exchange fraction Q f is the maximum value of this quantity over the entire time period of the simulation (0 < t < T ).