Tracer modelling within the brain parenchyma

In the previous sections, we explored how to generate a computational model of the brain parenchyma (\(\Omega_{PAR}\)) and the surrounding spaces fluid with CSF (\(\Omega_{CSF}\)). In this section, we will consider how tracers, such as Gabudutrol is moving through the brain parenchyma and CSF, given intrathecal injection.

The mathematical model for the concentration of a species \(c\) is given as

\[\begin{split} \begin{align} \frac{\partial (\phi c)}{\partial t} - \nabla \cdot (D \nabla (\phi c)) + \nabla \cdot (\mathbf{u}\phi c)&=0\quad \text{in } \Omega_{CSF}\cup\Omega_{PAR}&&\\ (-D \nabla (\phi c) + \mathbf{u} \phi c)\vert_{\Omega_{PAR}}\cdot \mathbf{n}\vert_{\Omega_{PAR}} &=\\ (-D \nabla (\phi c) + \mathbf{u} \phi c)\vert_{\Omega_{CSF}}\cdot \mathbf{n}\vert_{\Omega_{PAR}}&=\beta_{pia} (c_{PAR} - c_{CSF}) &&\text{on } \Gamma_{pia}\cup \Gamma_{LV}\\ (-D\nabla (\phi c) + \mathbf{u} \phi c) \cdot \mathbf{n} &= \begin{cases}-g_{influx} &&\text{on } \Gamma_{SSAS}\\ \beta_{exit} c &&\text{on } \Gamma_{AM-U}\\ 0 &&\text{on } \Gamma_{AM-L}\cup \Gamma_{SC}\\ \end{cases}&& \end{align} \end{split}\]

where \(\mathbf{n}\) is the outward pointing normal, \(\mathbf{n}_{\Omega_{PAR}}\) the outwards pointing normal from the brain parenchyma, and \(\phi\) is the fluid volume fraction, defined as

\[\begin{split} \begin{align} \phi = \begin{cases} 1&\quad\text{in } \Omega_{CSF}\\ \phi_{PAR}&\quad \text{in } \Omega_{PAR} \end{cases} \end{align} \end{split}\]

where \(\phi_{PAR}\ll 1\) and is a known quantity. \(\Omega_{PAR}\) refers to all the regions with brain tissue, while \(\Omega_{CSF}\) refers to the fluid cavities (ventricles) and networks within and surrounding (subarachnoic space) the brain. The choroid plexus is a region within the ventricles where cells produce and secrete cerebrospinal fluid (CSF). The dura matter, is the outer set of membranes between the SAS space and the human skull.

We will set \(\phi_{PAR}=0.2\).

The boundary/interface conditions above can be interpreted as: The brain is covered by Pia mater, which interfaces with the subarachnoid space and the ventrical system, which contains the cerebral spinal fluid. We let this membrane be semi-permeable, i.e the movement of Gadubutrol is related to the difference in velocity in the CSF and Parenchym

We additionally have that the effective diffusion coefficient \(D\) is defined as

\[\begin{split} \begin{align} D = \begin{cases} D_{CSF} &\quad\text{in } \Omega_{CSF}\\ D_{PAR}&\quad \text{in } \Omega_{PAR} \end{cases} \end{align} \end{split}\]

The finite element method formulation is based on A.1.3 in [CKMR25], using a weighted interior penalty discontinuous Galerkin method to model the diffusion, an upwinding scheme to stabilize the advective part of the equation, and a backward Euler temporal discretization. This results in the following variational formulation:

\[ a(c^n, c^{n-1}, v; \mathbf{u}) = l_{BC}(v) \]

where

\[ a(c^n, c^{n-1}, v; \mathbf{u}^n) = \sum_{K\in\mathcal{K}} \frac{1}{\delta t} \int_{K} \phi(c^n-c^{n-1})\cdot v ~\mathrm{d}x + a_{DG}(c^n, v) + a_{UP}(c^n, v \mathbf{u}^n) + a_{BC}(c^n, v). \]

Following is a breakdown of how to implement \(a\) and \(L\) in FEniCS using the unified form language UFL.

Interior penalty discontiuous Galerkin formulation

We start by defining the weighted \(\{\cdot\}_w\) average over a facet shared between two cells \(K_1\), \(K_2\) as

\[ \{v\}_{w} = \frac{w\vert_{K_2}}{w\vert_{K_1}+w\vert_{K_2}}v\vert_{K_1} + \frac{w\vert_{K_1}}{w\vert_{K_1}+w\vert_{K_2}}v\vert_{K_2}, \]

where an unweighed average is recovered in the case \(w=\frac{1}{2}\). In FEniCS we use ("+") and ("-") to restrict to \(K_1\) and \(K_2\) respectively. We use the shorthand \(\{\cdot\}\) to denote \(\{\cdot\}_{\frac{1}{2}}\).

import ufl


def avg(v: ufl.core.expr.Expr, w: ufl.core.expr.Expr = 0.5) -> ufl.core.expr.Expr:
    """Weighted average of an operator $v$.

    Note:
        If `w` is not supplied an unweighted average is returned
    """
    weight = ufl.as_ufl(w)  # Overload scalar valued kappa
    wp = weight("+") / (weight("+") + weight("-"))
    wm = weight("-") / (weight("+") + weight("-"))
    return wm * v("+") + wp * v("-")

Following [CKMR25], we use a symmetric interior penalty discontinuous Galerkin method with weighted residuals to discretize the diffusion term

\[ a_{DG}(c^n, v) = \sum_{K\in\mathcal{K}} \int_{K} \kappa \nabla c^n \cdot \nabla v - \sum_{f\in\mathcal{F}}\left(\{\kappa \nabla c^n \}_\kappa \cdot \mathbf{n}_F [v] + \{\kappa \nabla v \}_w \cdot \mathbf{n}_F [c^n] + \eta\frac{\gamma_{D,F}}{h_F}[c^n][v]\right)~\mathrm{d}s, \]

where \(\kappa=D\phi\). The set \(\mathcal{K}\) is the set of all tetrahedrons used to discretize \(\Omega_{CSF}\cup\Omega_{PAR}\), while \(\mathcal{F}\) is the set of all facets shared between two cells, excluding those on the interface \(\Gamma_{pia}\cup\Gamma_{LV}\).

When we want to restrict integrals to sub-sets of facets in the domain, we supply subdomain_data and subdomain_id to the integration measure. We create a convenience function that performs this task:

def define_volume_and_surface_measures(
    domain: ufl.Mesh,
    cell_marker=None,
    facet_marker=None,
):
    # Check for consistency of mesh markers
    tdim = domain.ufl_cell().topological_dimension()
    if cell_marker is not None and cell_marker.dim != tdim:
        raise ValueError(
            f"Volume marker must be of dimension {tdim}, but got {cell_marker.dim}"
        )
    if facet_marker is not None and facet_marker.dim != tdim - 1:
        raise ValueError(
            f"Surface marker must be of dimension {tdim - 1}, but got {facet_marker.dim}"
        )

    dx = ufl.dx(domain=domain, subdomain_data=cell_marker)
    ds = ufl.ds(domain=domain, subdomain_data=facet_marker)
    dS = ufl.dS(domain=domain, subdomain_data=facet_marker)
    return dx, ds, dS

def a_DG(
    c: ufl.Coefficient | ufl.Argument,
    v: ufl.Argument,
    kappa: ufl.core.expr.Expr,
    eta: ufl.core.expr.Expr,
    facet_markers,
    facet_values: tuple[int, ...],
) -> ufl.Form:
    """
    Args:
        c: The concentration field
        v: Test-function for concentration function space
        kappa: Diffusion coefficient weighted by porosity
        eta: Stability parameter for interior penalty/DG method
        facet_marker: Indicator for facets of the mesh
        facet_values: A list of values that indicates which facets in `facet_marker` the
            DG method should be applied to.
    """
    # Extract integration domain and create integration measures
    mesh = c.ufl_domain()
    dx, _, dS = define_volume_and_surface_measures(mesh, facet_marker=facet_markers)

    # Geometric quantities
    h_f = ufl.avg(2 * ufl.Circumradius(mesh))
    n_f = ufl.FacetNormal(mesh)("+")

    # Porosity weighted-diffusion coefficient
    kappa = ufl.as_ufl(kappa)  # Overload scalar valued kappa
    gamma = kappa("+") * kappa("-") / (kappa("+") + kappa("-"))

    # Variational form
    a = ufl.inner(kappa * ufl.grad(c), ufl.grad(v)) * dx
    a += gamma * eta / h_f * ufl.inner(ufl.jump(c), ufl.jump(v)) * dS(facet_values)
    a -= ufl.inner(ufl.dot(avg(kappa * ufl.grad(c), kappa), n_f), ufl.jump(v)) * dS(
        facet_values
    )
    a -= ufl.inner(ufl.dot(avg(kappa * ufl.grad(v), kappa), n_f), ufl.jump(c)) * dS(
        facet_values
    )
    return a

The advection term and upwinding scheme

Using the upwinding technique over the same facets \(\mathcal{F}_i\) as above, we get

\[ a_{UP}(c^n, v, \mathbf{u}^n) = - \sum_{k\in\mathcal{K}} \int_{k} c^n \mathbf{u}^n \cdot \nabla v ~\mathrm{d}x + \sum_{F\in \mathcal{F}}\int_{F} \{\mathbf{u}c\}\cdot \mathbf{n}_F + \frac{1}{2}\vert \mathbf{u}^n\cdot \mathbf{n}_F\vert[c][v]~\mathrm{d}s \]

where we note that for the last term to be well-defined \(\mathbf{u}\cdot\mathbf{n}\) has to be continuous across \(\mathcal{F}_i\).

The mathematical formulation above corresponds to the following code:

def a_up(
    c: ufl.Coefficient | ufl.Argument,
    v: ufl.Argument,
    u: ufl.core.expr.Expr,
    facet_markers,
    facet_values: tuple[int, ...],
) -> ufl.Form:
    """Return the upwinding residual form the advection-diffusion problem

    Note:
        It is assumed that u is in an H(div) function space, i.e. dot(u, n) is continuous over a facet.

    Args:
        c: The concentration field
        v: Test-function for concentration function space
        u: The velocity field
        facet_marker: Indicator for facets of the mesh
        facet_values: A list of indices that indicates which facets in `facet_marker` the
            DG method should be applied to.
    """
    mesh = c.ufl_domain()
    dx, _, dS = define_volume_and_surface_measures(mesh, facet_marker=facet_markers)
    n_f = ufl.FacetNormal(mesh)("+")

    F = -ufl.inner(u * c, ufl.grad(v)) * dx
    F += ufl.inner(avg(u * c), n_f) * ufl.jump(v) * dS(facet_values)

    u_n = ufl.dot(u("+"), n_f)
    F += ufl.inner(abs(u_n) / 2 * ufl.jump(c), ufl.jump(v)) * dS(facet_values)
    return F

Boundary and interface conditions

The boundary and interface conditions are resolved by integration by parts, and substituted to replace the arising terms over \(\Gamma_{pia}\cup\Gamma{LV}\), \(\Gamma_{SSAS}\) and \(\Gamma_{AM-U}\):

\[\begin{split} \begin{align*} a_{BC}(c^n, v) &= \sum_{F\in\Gamma_{pia}\cup\Gamma_{LV}}\int_{F} \beta_{pia}[c^n][v]~\mathrm{d}s + \sum_{F\in \Gamma_{AM-U}} \int_{F} \beta_{exit} c^n v~\mathrm{d}s\\ l(v)&= \sum_{F\in \Gamma_{SSAS}} g_{influx}v~\mathrm{d}s \end{align*} \end{split}\]

def a_bc(
    c: ufl.Coefficient | ufl.Argument,
    v: ufl.Argument,
    beta_pia: ufl.core.expr.Expr,
    beta_exit: ufl.core.expr.Expr,
    facet_marker,
    am_u: tuple[int, ...] | int,
    pia: tuple[int, ...] | int,
    lv: tuple[int, ...] | int,
) -> ufl.form.Form:
    """Returns the residual form of the interface and boundary conditions

    Args:
        c: The concentration field
        v: Test-function for concentration function space
        beta_pia: Permeability coefficient at the pial surface
        beta_exit: Permeability coefficient at the upper arachnoid and dura surface
        facet_marker: Indicator for facets of the mesh
        am_u: Tag(s) for the upper dura matter
        pia: Tag(s) for the pia matter
        lv: Tag(s) for the lateral ventricle
    """
    mesh = c.ufl_domain()
    _, ds, dS = define_volume_and_surface_measures(mesh, facet_marker=facet_marker)
    F = beta_pia * ufl.inner(ufl.jump(c), ufl.jump(v)) * (dS(pia) + dS(lv))
    F += beta_exit * ufl.inner(c, v) * ds(am_u)
    return F


def l_bc(
    g_influx: ufl.core.expr.Expr,
    v: ufl.Argument,
    facet_markers,
    ssas: tuple[int, ...],
):
    """Influx term over a set of marked facets

    Args:
        g_influx: Influx at the spinal subarachnoid space
        v: Test-function for concentration function space
        facet_markers: Indicator for facets of the mesh
        facet_values: A list of indices that indicates which facets in `facet_marker` the influx should be applied to.
    """
    mesh = v.ufl_domain()
    _, ds, _ = define_volume_and_surface_measures(mesh, facet_marker=facet_markers)
    return g_influx * v * ds(ssas)

With all these functions in place, we can load our previously generated mesh

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.shared_facet)
    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"
    )
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,
    "SSAS": 6,
}

As we would like to integrate over \(\mathcal{F}\), all those interior facets that are not on the interface between the brain and CSF flow, we extend our markers and interface map with additional tags.

import numpy as np

fdim = brain_with_csf.topology.dim - 1
fmap = brain_with_csf.topology.index_map(fdim)
num_facets_local = fmap.size_local + fmap.num_ghosts
interior_tag = 12
values = np.full(num_facets_local, interior_tag, dtype=np.int32)
for value in interface_map.values():
    values[interface_markers.find(value)] = value
extended_marker = dolfinx.mesh.meshtags(
    brain_with_csf, fdim, np.arange(num_facets_local, dtype=np.int32), values
)
del interface_markers
interface_map["DG"] = interior_tag

We define the problem specific parameters

Dh = dolfinx.fem.Constant(brain_with_csf, 1.2 * 1e-4)
phih = dolfinx.fem.Constant(brain_with_csf, 0.2)
beta_piah = dolfinx.fem.Constant(brain_with_csf, 2.6 * 1e-8)
beta_exith = dolfinx.fem.Constant(brain_with_csf, 1e-4)
eta = dolfinx.fem.Constant(brain_with_csf, 1.0)
dt = dolfinx.fem.Constant(brain_with_csf, 10.0)
gh = dolfinx.fem.Constant(brain_with_csf, 0.1)

Define the discrete function spaces for the concentration \(c^n\) and \(\mathbf{u}\). For simplicity, we choose \(\mathbf{u}=\mathbf{0}\) below:

V = dolfinx.fem.functionspace(brain_with_csf, ("DG", 1))
c_n = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
c_nm1 = dolfinx.fem.Function(V)
u = dolfinx.fem.Constant(
    brain_with_csf,
    np.zeros(brain_with_csf.geometry.dim, dtype=dolfinx.default_scalar_type),
)

Next, we define the porosity parameter \(\phi\)

import numpy as np

Q = dolfinx.fem.functionspace(brain_with_csf, ("DG", 0))
phi = dolfinx.fem.Function(Q)
par_cells = region_markers.indices[np.isin(region_markers.values, subdomain_map["PAR"])]
phi.x.array[par_cells] = 0.2
csf_cells = region_markers.indices[
    np.isin(
        region_markers.values,
        subdomain_map["LV"] + subdomain_map["SAS"] + subdomain_map["V34"],
    )
]
phi.x.array[csf_cells] = 1.0

With that definition at hand, we are ready to define \(a\) and \(L\)

dx = ufl.dx(domain=brain_with_csf)
kappa = Dh * phi
a = ufl.inner(c_n - c_nm1, v) * dx
a += dt * a_DG(c_n, v, kappa, eta, extended_marker, interface_map["DG"])
a += a_up(c_n, v, u, extended_marker, interface_map["DG"])
a += a_bc(
    c_n,
    v,
    beta_piah,
    beta_exith,
    extended_marker,
    interface_map["AM_U"],
    interface_map["V34_PAR"] + interface_map["PAR_SAS"],
    interface_map["LV_PAR"],
)
L = l_bc(gh, v, extended_marker, interface_map["SSAS"])