Use case 1: Molecular tracer diffusion in a human brain

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

\[\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. \]

We start by reading in the mesh generated in Generating brain meshes.

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

folder = Path(os.environ["WILDFENICS_DATA_PATH"])
assert folder.exists(), "Could not find surface files"

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, folder / "brain.xdmf", "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")

We define the following types for looking up information regarding subdomains or surfaces/interfaces where there are boundary conditions

import typing

subdomains = typing.Literal["SAS", "LV", "V34"]
interfaces = typing.Literal["LV_PAR", "V34_PAR", "AM_U", "AM_L", "EXTERNAL"]

The following maps the various regions to an integer marker

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,
}

The fluid regions can be gathered in a sequence of integers as

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

We will refine all cells (tetrahedra) whose facet is considered external, i.e. it does not connect to another fluid cell. We create the following convenience function to locate these cells.

import scifem
import numpy as np


def compute_subdomain_exterior_cells(
    mesh: dolfinx.mesh.Mesh, ct: dolfinx.mesh.MeshTags, markers: tuple[int, ...]
):
    """Compute the exterior boundary of a set of subdomains.

    Args:
        mesh: Mesh to extract subdomains from
        ct: MeshTags object marking subdomains
        markers: The tags making up the "new" mesh
    Returns:
        Cells which has a facet on the exterior boundary of the subdomains.
    """
    # Find facets that are considered exterior
    subdomain_exterior_facets = scifem.mesh.compute_subdomain_exterior_facets(
        mesh, ct, markers
    )
    tdim = mesh.topology.dim
    assert ct.dim == tdim
    sub_cells = dolfinx.mesh.compute_incident_entities(
        mesh.topology,
        subdomain_exterior_facets,
        tdim - 1,
        tdim,
    )
    full_subdomain = ct.indices[
        np.isin(ct.values, np.asarray(markers, dtype=ct.values.dtype))
    ]
    cell_map = mesh.topology.index_map(tdim)
    return scifem.mesh.reverse_mark_entities(
        cell_map, np.intersect1d(full_subdomain, sub_cells)
    )

With this function at hand, we can locally refine our mesh.

# Refine parent mesh within ventricles
num_refinements = 2
if num_refinements == 0:
    refined_mesh = brain_with_csf
    refined_regions = region_markers
    num_cells = brain_with_csf.topology.index_map(3).size_global
    print(f"Mesh, #Cells={num_cells}")
for i in range(num_refinements):
    # Refine parent mesh within ventricles
    refine_cells = region_markers.indices[
        np.isin(
            region_markers.values,
            np.asarray(subdomain_map["V34"] + subdomain_map["LV"]),
        )
    ]

    # Find all cells associated with outer boundary (dura) and refine the cells they correspond to
    brain_with_csf.topology.create_connectivity(
        brain_with_csf.topology.dim - 1, brain_with_csf.topology.dim
    )
    fmap = brain_with_csf.topology.index_map(brain_with_csf.topology.dim - 1)
    exterior_facet_indices = scifem.mesh.reverse_mark_entities(
        fmap, dolfinx.mesh.exterior_facet_indices(brain_with_csf.topology)
    )
    boundary_cells = dolfinx.mesh.compute_incident_entities(
        brain_with_csf.topology,
        exterior_facet_indices,
        brain_with_csf.topology.dim - 1,
        brain_with_csf.topology.dim,
    )

    fluid_boundary_cells = compute_subdomain_exterior_cells(
        brain_with_csf, region_markers, fluid_domains
    )

    # For any further refinement, only refine the boundary of the fluid domains, not the interior
    if i < 1:
        cells_to_refine = np.unique(
            np.hstack([boundary_cells, fluid_boundary_cells, refine_cells])
        ).astype(np.int32)

    else:
        cells_to_refine = refine_cells

    edges_to_refine = dolfinx.mesh.compute_incident_entities(
        brain_with_csf.topology, cells_to_refine, brain_with_csf.topology.dim, 1
    )
    edge_map = brain_with_csf.topology.index_map(1)
    edges_to_refine = scifem.mesh.reverse_mark_entities(edge_map, edges_to_refine)
    refined_mesh, parent_cell, _ = dolfinx.mesh.refine(
        brain_with_csf,
        edges_to_refine,
        partitioner=None,
        option=dolfinx.mesh.RefinementOption.parent_cell,
    )
    refined_regions = dolfinx.mesh.transfer_meshtag(
        region_markers, refined_mesh, parent_cell
    )
    brain_with_csf = refined_mesh
    region_markers = refined_regions
    num_cells = brain_with_csf.topology.index_map(3).size_global
    print(f"Mesh refinement {i + 1}, #Cells={num_cells}")

Mesh refinement 1, #Cells=1405474

Mesh refinement 2, #Cells=1724676

We need to ensure that the mesh is sufficiently fine to run fluid simulations on. We ensure this by using local refinement of the ventricular system and the regions of the CSF space that is either close to the brain parenchyma or the dura matter.

References

[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.

[NGB+13]

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.

[RGB+85]

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.

[ZRKW92]

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.