Fluid flow in the ventricular system and subarachnoid spaces

Introduction

In many medical studies, tracers are used to understand how solutes move through the human brain [], [], []. They are preferred as one can use medical imaging techniques, such as PET and MRI, while CT uses X-ray images to construct a 3D model of an organ.

Mathematical model

We start by describing the model of cerebrospinal fluid (CSF), based on []. This network consist of two connected regions, namely the ventricular system \(\Omega_V\) and the subarachnoid space \(\Omega_{SAS}\). The subarachnoid space is the space between the membrane of the brain (pia matter) and the membrane towards the skull (dura matter). The ventrical system is a cavity within the brain where cells produce CSF. The CSF enters the brain at the top of the spinal coord. We will call this boundary \(\Gamma_{SSAS}\) (Spinal subarachnoid space). We consider the brain paraenchyma (the brain tissue) as being surrounded by CSF. The interface between the parenchyma and the CSF can be split into two regions, the surface of the lateral ventricles \(\Gamma_{LV}\) and the pial surface \(\Gamma_{pia}\). We consider no fluid flow through the pial surface. However, the ventricles produce fluid, and therefore we consider an influx through \(\Gamma_{LV}\). The CSF is surrounded by dura matter, which is denoted as \(\Gamma_{AM}\). This region is subdivided into an upper (\(\Gamma_{AM-U}\)) and lower region (\(\Gamma_{AM-L}\)). Fluid can penetrate the the upper region (and exit the system), while no flow occurs through the lower region. Summarizing, we can write the fluid equations as

\[\begin{split} \begin{align} \mu \Delta \mathbf{u} + \nabla p &= 0 &&\text{in } \Omega_{CSF}\\ \nabla \cdot \mathbf{u} &=0 && \text{in } \Omega_{CSF}\\ \mathbf{u} &= 0 &&\text{in } \Omega_{PAR}\\ -\mu \left((\nabla \mathbf{u})\cdot\mathbf{n} + p \mathbf{n}\right)\cdot \mathbf{n}&= R_0(\mathbf{u}\cdot \mathbf{n}) && \text{on } \Gamma_{AM-U}\\ \mathbf{u} \cdot \mathbf{t} &=0 &&\text{on }\Gamma_{AM-U}\\ \mathbf{u}&= 0 &&\text{on } \Gamma_{AM-L}\cup\Gamma_{pia}\cup\Gamma_{SSAS}\\ \mathbf{u}\cdot \mathbf{n} &=\frac{1}{\vert\Gamma_{LV}\vert}u_{in} && \text{on } \Gamma_{LV}\\ \mathbf{u}\cdot\mathbf{t} &= 0 && \text{on } \Gamma_{LV} \end{align} \end{split}\]

We solve the set of linear equation with the FEM and the following variational form:

Find \(u\in H_g^1(\mathrm{div}, \Omega_{\mathrm{CSF}})\), \(p\in L^2(\Omega_{\mathrm{CSF}})\) such that

\[ a(u, p, v, q) = L(v, q) \quad \forall v\in H_0^1(\mathrm{div}, \Omega_{\mathrm{CSF}}), q\in L^2(\Omega_{\mathrm{CSF}}) \]

where

\[\begin{split} \begin{align*} a(u, p, v, q)&= \int_{\Omega_{CSF}} \mu \nabla \mathbf{u} \cdot \nabla \mathbf{v} - \nabla \cdot \mathbf{v} p ~\mathrm{d}x + \int_{\Gamma_{AM-u}} R_0 \mathbf{u}\cdot \mathbf{n} \mathbf{n} \cdot\mathbf{v} ~\mathrm{d}s+\\ &+\int_{\mathcal{F}} -\mu (\{\nabla \mathbf{u}\}\cdot\mathbf{n}_F) \cdot \left[\mathbf{v}_t\right] -\mu (\{\nabla \mathbf{v}\}\cdot \mathbf{n}_F) \cdot \left[\mathbf{u}_t\right] + \frac{\sigma \mu}{h_F}\left[\mathbf{u}_t\right]~\mathrm{d}s\\ &+\int_{\Omega_{CSF}} \nabla\cdot \mathbf{u}q~\mathrm{d} x = 0. \end{align*} \end{split}\]

and

\[ L(v, q) = 0. \]

Reading in the mesh and creating submeshes

We start by reading in the mesh generated in Generating brain meshes. The refined mesh requires 80 GB RAM for running this demo, while the un-refined mesh requires around 25 GB RAM.

from mpi4py import MPI
import dolfinx
from pathlib import Path
import os

use_refined_mesh = True
if use_refined_mesh:
    mesh_file = "refined_marked_brain.xdmf"
else:
    mesh_file = "marked_brain.xdmf"
folder = Path(os.environ["WILDFENICS_DATA_PATH"])
assert folder.exists(), "Could not find surface files"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, folder / mesh_file, "r") as xdmf:
    brain_with_csf = xdmf.read_mesh(dolfinx.cpp.mesh.GhostMode.none)
    region_markers = xdmf.read_meshtags(brain_with_csf, name="mesh_tags")
    brain_with_csf.topology.create_connectivity(
        brain_with_csf.topology.dim, brain_with_csf.topology.dim - 1
    )
    interface_markers = xdmf.read_meshtags(
        brain_with_csf, name="interfaces_and_boundaries"
    )

We define the same maps for easy look-up of the various marked regions.

subdomain_map = {
    "PAR": (2,),
    "SAS": (1,),
    "LV": (3,),
    "V34": (4,),
}
interface_map = {
    "LV_PAR": 1,
    "V34_PAR": 2,
    "PAR_SAS": 5,
    "AM_U": 3,
    "AM_L": 4,
}

As we would only like to simulate the fluid flow in the ventricles and SAS, and not through the brain parenchyma, we extract these cells into a new mesh. We also transfer the interface markers.

import scifem

fluid_domains = subdomain_map["LV"] + subdomain_map["SAS"] + subdomain_map["V34"]

csf_mesh, cell_map, vertex_map, node_map, csf_markers = scifem.extract_submesh(
    brain_with_csf,
    region_markers,
    fluid_domains,
)
fluid_interface_marker, _ = scifem.transfer_meshtags_to_submesh(
    interface_markers, csf_mesh, vertex_map, cell_map
)
del region_markers, interface_markers

Defining the variational formulation

Next, we create the mixed function space, using “BDM”-1 and “DG”-0 as the stable finite element pair as illustrated in []:

import basix.ufl
import ufl
import numpy as np

degree = 1
element_u = basix.ufl.element(basix.ElementFamily.BDM, csf_mesh.basix_cell(), degree)
element_p = basix.ufl.element(
    basix.ElementFamily.P,
    csf_mesh.basix_cell(),
    degree=degree - 1,
    discontinuous=True,
)
me = basix.ufl.mixed_element([element_u, element_p])
W = dolfinx.fem.functionspace(csf_mesh, me)

Next, we define the integration measures we will use in the variational formulation

dx = ufl.dx(domain=csf_mesh, subdomain_data=csf_markers)  # All fluid cells
ds = ufl.ds(
    domain=csf_mesh, subdomain_data=fluid_interface_marker
)  # All exterior facets
dAM_U = ds(interface_map["AM_U"])  # Upper skull
dWall = ds(
    (
        interface_map["V34_PAR"],
        interface_map["PAR_SAS"],
        interface_map["AM_L"],
        interface_map["LV_PAR"],
    )  # Remaining exterior faces
)
dS = ufl.dS(domain=csf_mesh)  # All interior facets

Next, we define some of the problem specific parameters

mu = dolfinx.fem.Constant(csf_mesh, dolfinx.default_scalar_type(7e-3))
R = dolfinx.fem.Constant(csf_mesh, dolfinx.default_scalar_type(1e4))
sigma = dolfinx.fem.Constant(csf_mesh, dolfinx.default_scalar_type(100.0))
U_in = dolfinx.fem.Constant(csf_mesh, dolfinx.default_scalar_type(4.63e-7))

The tangential projection operator is defined as

def tangent_projection(u, n):
    return u - ufl.dot(u, n) * n

giving us the intermediate variables

n = ufl.FacetNormal(csf_mesh)
u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)
u_t = tangent_projection(u, n)
v_t = tangent_projection(v, n)

which yields the following variational form

a = mu * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx - ufl.div(v) * p * dx
a += -ufl.div(u) * q * dx
a += R * ufl.dot(u, n) * ufl.dot(v, n) * dAM_U
hF = ufl.FacetArea(csf_mesh)
hA = ufl.avg(2.0 * ufl.Circumradius(csf_mesh))
# Wall condition (slip condition)
a += (
    -ufl.inner(ufl.dot(mu * ufl.grad(v), n), u_t) * dWall
    - ufl.inner(ufl.dot(mu * ufl.grad(u), n), v_t) * dWall
    + 2 * mu * sigma / hF * ufl.inner(u_t, v_t) * dWall
)
# Weak enforcement of tangential continuity
a += (
    -ufl.inner(
        ufl.dot(ufl.avg(mu * ufl.grad(u)), n("+")),
        ufl.jump(v_t),
    )
    - ufl.inner(
        ufl.dot(ufl.avg(mu * ufl.grad(v)), n("+")),
        ufl.jump(u_t),
    )
    + 2 * mu * (sigma / hA) * ufl.inner(ufl.jump(u_t), ufl.jump(v_t))
) * dS

# Zero right-hand side
f = dolfinx.fem.Constant(
    csf_mesh, dolfinx.default_scalar_type(np.zeros(csf_mesh.geometry.dim))
)
L = ufl.inner(f, v) * dx

Defining strong boundary conditions

In the formulation, we have that \(\mathbf{u}\cdot\mathbf{n} = \frac{1}{\vert\Gamma_{LV}\vert}u_{in}\). As we are using BDM-spaces to enforce this, we require some special handling:

With this function in hand, we can define the Dirichlet condition. We first compute the area \(\Gamma_{LV}\):

surface_area_form = 1 * ds(interface_map["LV_PAR"])
surface_area_value = scifem.assemble_scalar(surface_area_form)
gamma_LV = dolfinx.fem.Constant(csf_mesh, surface_area_value)

Next, we create the symbolic expression to enforce

inlet_expr = -U_in / gamma_LV * n

We find the facets on the interface of the left ventricle

LV_facets = fluid_interface_marker.indices[
    np.isin(fluid_interface_marker.values, interface_map["LV_PAR"])
]

and finally create the Dirichlet boundary condition

V, _ = W.sub(0).collapse()
u_bc = scifem.interpolate_function_onto_facet_dofs(
    V,
    inlet_expr,
    LV_facets,
)
LV_dofs = dolfinx.fem.locate_dofs_topological(
    (W.sub(0), V), fluid_interface_marker.dim, LV_facets
)
bc_LV = dolfinx.fem.dirichletbc(u_bc, LV_dofs, W.sub(0))

We additionally create the slip conditions at the remaining walls, where \(\mathbf{u}\cdot \mathbf{n} = 0\)

u_wall = dolfinx.fem.Function(V)
u_wall.x.array[:] = 0.0
wall_facets = fluid_interface_marker.indices[
    np.isin(
        fluid_interface_marker.values,
        [interface_map["V34_PAR"], interface_map["PAR_SAS"], interface_map["AM_L"]],
    )
]
wall_dofs = dolfinx.fem.locate_dofs_topological(
    (W.sub(0), V), fluid_interface_marker.dim, wall_facets
)
bc_wall = dolfinx.fem.dirichletbc(u_wall, wall_dofs, W.sub(0))
bcs = [bc_LV, bc_wall]

Solving the linear variational problem

We use PETSc and MUMPS [] to solve the arising linear system. We turn on several optimizations to speed up the assembly process.

jit_options = {
    "cffi_extra_compile_args": ["-Ofast", "-march=native"],
    "cffi_libraries": ["m"],
}
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "ksp_error_if_not_converged": True,
    "pc_factor_mat_solver_type": "mumps",
    "mat_mumps_icntl_14": 250,
    "mat_mumps_icntl_24": 1,
    "mat_mumps_icntl_4": 2,
}
wh = dolfinx.fem.Function(W)
problem = dolfinx.fem.petsc.LinearProblem(
    a,
    L,
    u=wh,
    bcs=bcs,
    petsc_options=petsc_options,
    jit_options=jit_options,
    petsc_options_prefix="StokesSolver_",
)

We only solve the problem if we have sufficient amount of RAM on our system.

import psutil

if use_refined_mesh:
    min_remaining_ram_required = 80 * 10**9
else:
    min_remaining_ram_required = 25 * 10**9
memory_stats = psutil.virtual_memory()
if memory_stats.available > min_remaining_ram_required:
    problem.solve()
else:
    # Implement a checkpoint here
    pass

As we would like to use the solution on the parent mesh, that includes the brain parenchyma, we transfer the solution to this domain, and store it as a “checkpoint” that can be read in alongside the mesh from the input file.

import adios4dolfinx

u_in_fluid = wh.sub(0).collapse()
V_full = dolfinx.fem.functionspace(brain_with_csf, element_u)
function_filename = Path(os.environ["WILDFENICS_DATA_PATH"]) / "velocity_checkpoint.bp"

# To transfer the solution to the parent, we use interpolation between the parent and child mesh
u_full = dolfinx.fem.Function(V_full, name="u")
tdim = csf_mesh.topology.dim
num_fluid_cells = (
    csf_mesh.topology.index_map(tdim).size_local
    + csf_mesh.topology.index_map(tdim).num_ghosts
)
fluid_cells = np.arange(num_fluid_cells, dtype=np.int32)
submesh_to_mesh = cell_map.sub_topology_to_topology(fluid_cells, inverse=False)
child_cells = np.arange(len(submesh_to_mesh), dtype=np.int32)
u_full.interpolate(u_in_fluid, cells0=fluid_cells, cells1=submesh_to_mesh)
adios4dolfinx.write_function_on_input_mesh(
    function_filename.with_suffix(".bp"),
    u_full,
    mode=adios4dolfinx.adios2_helpers.adios2.Mode.Write,
    time=0.0,
)

To visualize the solution, we transfer the velocity to a compatible discontinuous Lagrange space

# For visualization of fluid flow within fluid cavities
V_out = dolfinx.fem.functionspace(
    brain_with_csf, ("DG", degree, (csf_mesh.geometry.dim,))
)
u_out = dolfinx.fem.Function(V_out, name="Velocity")
u_out.interpolate(u_full)
with dolfinx.io.VTXWriter(csf_mesh.comm, "uh.bp", [u_out], engine="BP4") as bp:
    bp.write(0.0)

References

[BFRSC24]

Pietro Benedusi, Paola Ferrari, Marie E. Rognes, and Stefano Serra-Capizzano. Modeling excitable cells with the emi equations: spectral analysis and iterative solution strategy. Journal of Scientific Computing, 98(3):58, Feb 2024. doi:10.1007/s10915-023-02449-2.

[CKMR25]

Marius Causemann, Miroslav Kuchta, Rami Masri, and Marie E. Rognes. In-silico molecular enrichment and clearance of the human intracranial space. bioRxiv, 2025. doi:10.1101/2025.01.30.635680.

[KMR21]

Miroslav Kuchta, Kent-André Mardal, and Marie E. Rognes. Solving the EMI Equations using Finite Element Methods, pages 56–69. Springer International Publishing, Cham, 2021. doi:10.1007/978-3-030-61157-6_5.

[TMR21]

Aslak Tveito, Kent-Andre Mardal, and Marie E. Rognes, editors. Modeling Excitable Tissue: The EMI Framework. Springer International Publishing, Cham, 2021. ISBN 978-3-030-61157-6. doi:10.1007/978-3-030-61157-6.