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 [RGB+85], [ZRKW92], [NGB+13]. 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 [CKMR25]. 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
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
where
and
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.mesh.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 [HKXZ16]:
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:
from packaging.version import Version
import numpy.typing as npt
def strong_bc_bdm_function(
Q: dolfinx.fem.FunctionSpace,
expr: ufl.core.expr.Expr,
facets: npt.NDArray[np.int32],
) -> dolfinx.fem.Function:
"""
Create a function $u_h\in Q$ such that $u_h=\text{expr}$ for all dofs belonging to a subset of ``facets``.
All other dofs are set to zero.
Args:
Q: The function space to create the function $u_h$ in.
expr: The expression to evaluate.
facets: The facets on which to evaluate the expression.
"""
domain = Q.mesh
Q_el = Q.element
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
# Compute integration entities (cell, local_facet index) for all facets
if Version(dolfinx.__version__) > Version("0.9.0"):
boundary_entities = dolfinx.fem.compute_integration_domains(
dolfinx.fem.IntegralType.exterior_facet, domain.topology, facets
)
else:
boundary_entities = dolfinx.fem.compute_integration_domains(
dolfinx.fem.IntegralType.exterior_facet,
domain.topology,
facets,
domain.topology.dim - 1,
)
interpolation_points = Q_el.basix_element.x
fdim = domain.topology.dim - 1
c_el = domain.ufl_domain().ufl_coordinate_element()
ref_top = c_el.reference_topology
ref_geom = c_el.reference_geometry
facet_types = set(
basix.cell.subentity_types(domain.basix_cell())[domain.topology.dim - 1]
)
assert len(facet_types) == 1, "All facets must have the same topology"
# Pull back interpolation points from reference coordinate element to facet reference element
facet_cmap = basix.ufl.element(
"Lagrange",
facet_types.pop(),
c_el.degree,
shape=(domain.geometry.dim,),
dtype=np.float64,
)
facet_cel = dolfinx.cpp.fem.CoordinateElement_float64(facet_cmap.basix_element._e)
reference_facet_points = None
for i, points in enumerate(interpolation_points[fdim]):
geom = ref_geom[ref_top[fdim][i]]
ref_points = facet_cel.pull_back(points, geom)
# Assert that interpolation points are all equal on all facets
if reference_facet_points is None:
reference_facet_points = ref_points
else:
assert np.allclose(reference_facet_points, ref_points)
# Create expression for BC
normal_expr = dolfinx.fem.Expression(expr, reference_facet_points)
points_per_entity = [sum(ip.shape[0] for ip in ips) for ips in interpolation_points]
offsets = np.zeros(domain.topology.dim + 2, dtype=np.int32)
offsets[1:] = np.cumsum(points_per_entity[: domain.topology.dim + 1])
values_per_entity = np.zeros(
(offsets[-1], domain.geometry.dim), dtype=dolfinx.default_scalar_type
)
entities = boundary_entities.reshape(-1, 2)
values = np.zeros(entities.shape[0] * offsets[-1] * domain.geometry.dim)
for i, entity in enumerate(entities):
insert_pos = offsets[fdim] + reference_facet_points.shape[0] * entity[1]
# Backwards compatibility
try:
normal_on_facet = normal_expr.eval(domain, entity.reshape(1, 2))
except (AttributeError, AssertionError):
normal_on_facet = normal_expr.eval(domain, entity)
# NOTE: evaluate within loop to avoid large memory requirements
values_per_entity[insert_pos : insert_pos + reference_facet_points.shape[0]] = (
normal_on_facet.reshape(-1, domain.geometry.dim)
)
values[
i * offsets[-1] * domain.geometry.dim : (i + 1)
* offsets[-1]
* domain.geometry.dim
] = values_per_entity.reshape(-1)
qh = dolfinx.fem.Function(Q)
qh._cpp_object.interpolate(
values.reshape(-1, domain.geometry.dim).T.copy(), boundary_entities[::2].copy()
)
qh.x.scatter_forward()
return qh
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 = strong_bc_bdm_function(
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 [ADLExcellentK01] 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
)
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 = 20 * 10**9
else:
min_remaining_ram_required = 80 * 10**9
memory_stats = psutil.virtual_memory()
if memory_stats.available > min_remaining_ram_required:
wh = 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")
child_cells = np.arange(len(cell_map), dtype=np.int32)
u_full.interpolate(u_in_fluid, cells0=child_cells, cells1=cell_map)
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#
Patrick R. Amestoy, Iain S. Duff, Jean-Yves L'Excellent, and Jacko Koster. Mumps: a general purpose distributed memory sparse solver. In Tor Sørevik, Fredrik Manne, Assefaw Hadish Gebremedhin, and Randi Moe, editors, Applied Parallel Computing. New Paradigms for HPC in Industry and Academia, 121–130. Berlin, Heidelberg, 2001. Springer Berlin Heidelberg. doi:10.1007/3-540-70734-4_16.
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.
Qingguo Hong, Johannes Kraus, Jinchao Xu, and Ludmil Zikatanov. A robust multigrid method for discontinuous galerkin discretizations of stokes and linear elasticity equations. Numerische Mathematik, 132(1):23–49, Jan 2016. doi:10.1007/s00211-015-0712-y.
Ruiqing Ni, Per-Göran Gillberg, Assar Bergfors, Amelia Marutle, and Agneta Nordberg. Amyloid tracers detect multiple binding sites in Alzheimer's disease brain tissue. Brain, 136(7):2217–2227, 06 2013. doi:10.1093/brain/awt142.
Marshall L. Rennels, Thomas F. Gregory, Otis R. Blaumanis, Katsukuni Fujimoto, and Patricia A. Grady. Evidence for a "paravascular" fluid circulation in the mammalian central nervous system, provided by the rapid distribution of tracer protein throughout the brain from the subarachnoid space. Brain Research, 326(1):47–63, 1985. doi:https://doi.org/10.1016/0006-8993(85)91383-6.
E. T. Zhang, H. K. Richards, S. Kida, and R. O. Weller. Directional and compartmentalised drainage of interstitial fluid and cerebrospinal fluid from the rat brain. Acta Neuropathologica, 83(3):233–239, Feb 1992. doi:10.1007/BF00296784.