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}\]

With the fluid model at hand, we can describe the spread of tracer throughout the CSF and brain parenchyma:

\[\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} &=\\ -(-D \nabla (\phi c) + \mathbf{u} \phi c)\vert_{\Omega_{CSF}}\cdot\mathbf{n}&=\beta_{pia} (c_{PAR} - c_{CSF})\qquad \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 \(c\) is the tracer (Gadubutrol) concentration field, \(\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}\]

For modelling the equations above, we will employ the finite element method (FEM). To use this method, we transform the equations into their variational form.

\[\begin{split} \int_{\Omega_{CSF}} -\mu \Delta \mathbf{u}\cdot \mathbf{v}~\mathrm{d}x + \int_{\Omega_{CSF}} \nabla p \cdot \mathbf{v}~\mathrm{d}x = 0\\ \int_{\Omega_{CSF}} \nabla\cdot \mathbf{u}q~\mathrm{d} x = 0 \end{split}\]

Integrating the first equation by parts, and choosing \(\mathbf{u}\in H_{div}(\Omega)\), we get

\[ \int_{\Omega_{CSF}} \mu \nabla \mathbf{u} \cdot \nabla \mathbf{v} - \nabla \cdot \mathbf{v} p ~\mathrm{d}x + \int_{\partial\Omega_{CSF}} \left(-\mu \nabla \mathbf{u} \cdot \mathbf{n} + p\mathbf{n}\right) \cdot \mathbf{v} ~\mathrm{d}s=0. \]

We observe that we can split the boundary integral into its normal and tangential component \(\mathbf{A}_t = \mathbf{A} - (\mathbf{A}\cdot \mathbf{n}) \mathbf{n}\)

\[ \int_{\partial\Omega_{CSF}} \left(-\mu \nabla \mathbf{u} \cdot \mathbf{n} + p\mathbf{n}\right) \cdot \mathbf{v} ~\mathrm{d}s= \int_{\partial\Omega_{CSF}} \left(-\mu \nabla \mathbf{u} \cdot \mathbf{n} + p\mathbf{n}\right)\cdot \mathbf{n} \mathbf{n} \cdot\mathbf{v} + \left(-\mu \nabla \mathbf{u} \cdot \mathbf{n} + p\mathbf{n}\right)_t \cdot \mathbf{v} ~\mathrm{d}x. \]

As we enforce \(\mathbf{u}\cdot n=\mathbf{0}\) on \(\Gamma_{AM-L}\cup\Gamma_{pia}\cup\Gamma_{SSAS}\) and \(\mathbf{u}\cdot n =\frac{1}{\vert \Gamma \vert}u_{in}\) through our choice of function space, these integrals dissapear on these boundaries.

On \(\Gamma_{AM-u}\) we use the boundary condition relation to rewrite the equation above as

\[ \int_{\partial\Omega_{CSF}} \left(-\mu \nabla \mathbf{u} \cdot \mathbf{n} + p\mathbf{n}\right) \cdot \mathbf{v} ~\mathrm{d}s= \int_{\Gamma_{AM-u}} R_0 \mathbf{u}\cdot \mathbf{n} \mathbf{n} \cdot\mathbf{v}~\mathrm{d}s + \int_{\partial \Omega_{CSF}} \left(-\mu \nabla \mathbf{u} \cdot \mathbf{n} + p\mathbf{n}\right)_t \cdot \mathbf{v} ~\mathrm{d}x. \]

It is clear that \((p\mathbf{n})_t=0\). As \(\mathbf{u}\cdot \mathbf{t} =0\) across the whole boundary, we have that \((\nabla \mathbf{u})_t=0\), and we are left with only the first integral.

This leaves us with the variational form

\[\begin{split} \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=0\\ \int_{\Omega_{CSF}} \nabla\cdot \mathbf{u}q~\mathrm{d} x = 0 \end{split}\]

Discretized problem

When the domain is discretized into tetrahedra, and a discrete \(H_{div}\) space is chosen, one cannot guarantee the tangential continuity of \(\mathbf{u}\) across element surfaces \(\mathcal{F}\). To enforce this, we add a standard interior penalty method of the form

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

try:
    import ufl
except ImportError:
    import ufl_legacy as ufl

from typing import TypedDict


def weighted_avg(a, kappa):
    wp = kappa("+") / (kappa("+") + kappa("-"))
    wm = kappa("-") / (kappa("+") + kappa("-"))
    return wp * a("+") + wm * a("-")


class MarkerLookup(TypedDict):
    pia: tuple[int]  # Pial surface
    lv: tuple[int]  # Lateral ventricle
    am_u: tuple[int]  # Arachnoid and dura upper
    am_l: tuple[int]  # Arachnoid and dura lower
    ssas: tuple[int]  # Spinal subarachnoid space


def diffusion_residual(
    c: ufl.Coefficient | ufl.Argument,
    v: ufl.Argument,
    D: ufl.core.expr.Expr,
    phi: ufl.core.expr.Expr,
    eta: ufl.core.expr.Expr,
) -> ufl.Form:
    """Returns the residual formulation of the "standard" DG method for a diffusion problem
    with a weighted average over the flux conditions.

    Args:
        c: Unknown concentration field
        v: Test-function for concentration function space
        D: Diffusion coefficient
        phi: Fluid-solid fraction
        eta: Stability parameter for interior penalty/DG method
    """
    mesh = c.ufl_domain()

    # Integration measures
    dx = ufl.Measure("dx")
    dS = ufl.Measure("dS")

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

    # Porosity weighted-diffusion coefficient
    kappa = D * phi
    gamma = kappa("+") * kappa("-") / (kappa("+") + kappa("-"))
    # Variational form
    F = ufl.inner(D * ufl.grad(c), ufl.grad(v)) * dx
    F += gamma * eta / h_f * ufl.inner(ufl.jump(c), ufl.jump(v)) * ufl.dS
    F -= (
        ufl.inner(
            ufl.dot(weighted_avg(D * phi * ufl.grad(c), D * phi), n_f), ufl.jump(v)
        )
        * dS
    )
    F -= (
        ufl.inner(
            ufl.dot(weighted_avg(D * phi * ufl.grad(v), D * phi), n_f), ufl.jump(c)
        )
        * dS
    )
    return F


def boundary_and_interface_conditions(
    c: ufl.Coefficient | ufl.Argument,
    v: ufl.Argument,
    beta_pia: ufl.core.expr.Expr,
    beta_exit: ufl.core.expr.Expr,
    g_influx: ufl.core.expr.Expr,
    markers: MarkerLookup,
) -> ufl.form.Form:
    """Returns the residual form of the interface and boundary conditions

    Args:
        c: The unknown 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
        g_influx: Influx at the spinal subarachnoid space
        markers: Lookup table for boundary and interface markers        _
    """
    dS = ufl.Measure("dS")
    ds = ufl.Measure("ds")
    F = (
        beta_pia
        * ufl.inner(ufl.jump(c), ufl.jump(v))
        * dS(markers["pia"] + markers["lv"])
    )
    F += beta_exit * ufl.inner(c, v) * ds(markers["am_u"])
    F -= g_influx * v * ds(markers["ssas"])
    return F


def upwinding(
    c: ufl.Coefficient | ufl.Argument, v: ufl.Argument, u: ufl.core.expr.Expr
) -> 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 unknown concentration field
        v: Test-function for concentration function space
        u: The velocity field
    """
    dx = ufl.Measure("dx")
    dS = ufl.Measure("dS")
    mesh = c.ufl_domain()
    n_f = ufl.FacetNormal(mesh)("+")

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

    # NOTE: Continuity assumption on u
    u_n = ufl.dot(u("+"), n_f)
    F += ufl.inner(abs(u_n) / 2 * ufl.jump(c), ufl.jump(v)) * dS
    return F


def advection_diffusion_residual(
    c: ufl.Coefficient,
    c_n: ufl.Coefficient,
    v: ufl.Argument,
    D: ufl.core.expr.Expr,
    phi: ufl.core.expr.Expr,
    eta: ufl.core.expr.Expr,
    beta_pia: ufl.core.expr.Expr,
    beta_exit: ufl.core.expr.Expr,
    tau: ufl.core.expr.Expr,
    g_influx: ufl.core.expr.Expr,
    u: ufl.Coefficient,
    markers: MarkerLookup,
) -> ufl.Form:
    """Return the residual form of the advection-diffusion problem

    Args:
        c: The unknown concentration field
        c_n: The concentration field at the previous time-step
        v: Test-function for concentration function space
        D: Diffusion coefficient
        phi: Fluid-solid fraction
        eta: Stability parameter for interior penalty/DG method
        beta_pia: Permeability coefficient at the pial surface
        beta_exit: Permeability coefficient at the upper arachnoid and dura surface
        tau: Time-step
        g_influx: Influx at the spinal subarachnoid space
        markers: Lookup table for boundary and interface markers
        u: The velocity field
    """
    F = (
        1 / tau * phi * ufl.inner(c - c_n, v) * ufl.dx
        + diffusion_residual(c, v, D, phi, eta)
        + boundary_and_interface_conditions(
            c, v, beta_pia, beta_exit, g_influx, markers
        )
        + upwinding(c, v, u)
    )
    return F


def advection_diffusion_system(
    c: ufl.Argument,
    c_n: ufl.Coefficient,
    v: ufl.Argument,
    D: ufl.core.expr.Expr,
    phi: ufl.core.expr.Expr,
    eta: ufl.core.expr.Expr,
    beta_pia: ufl.core.expr.Expr,
    beta_exit: ufl.core.expr.Expr,
    tau: ufl.core.expr.Expr,
    g_influx: ufl.core.expr.Expr,
    u: ufl.Coefficient,
    markers: MarkerLookup,
) -> tuple[ufl.Form, ufl.Form]:
    """Return the residual form of the advection-diffusion problem

    Args:
        c: Trial-function for the concentration function space
        c_n: The concentration field at the previous time-step
        v: Test-function for concentration function space
        D: Diffusion coefficient
        phi: Fluid-solid fraction
        eta: Stability parameter for interior penalty/DG method
        beta_pia: Permeability coefficient at the pial surface
        beta_exit: Permeability coefficient at the upper arachnoid and dura surface
        tau: Time-step
        g_influx: Influx at the spinal subarachnoid space
        markers: Lookup table for boundary and interface markers
        u: The velocity field)
    """
    F = advection_diffusion_residual(
        c, c_n, v, D, phi, eta, beta_pia, beta_exit, tau, g_influx, u, markers
    )
    return ufl.system(F)


if __name__ == "__main__":
    import basix.ufl

    # Symbolic variational form (used with C++ or compile_form)
    c_el = basix.ufl.element("Lagrange", basix.CellType.tetrahedron, 1, shape=(3,))
    domain = ufl.Mesh(c_el)
    el = basix.ufl.element(
        "Lagrange", domain.ufl_coordinate_element().cell_type, 2, discontinuous=True
    )
    velocity_el = basix.ufl.element("BDM", domain.ufl_coordinate_element().cell_type, 1)

    V = ufl.FunctionSpace(domain, el)
    Q = ufl.FunctionSpace(domain, velocity_el)
    u = ufl.Coefficient(Q)

    c = ufl.Coefficient(V)
    c_n = ufl.Coefficient(V)
    v = ufl.TestFunction(V)

    D = ufl.Constant(domain)
    phi = ufl.Constant(domain)
    eta = ufl.Constant(domain)
    beta_pia = ufl.Constant(domain)
    beta_exit = ufl.Constant(domain)
    tau = ufl.Constant(domain)
    g_influx = ufl.Constant(domain)

    markers = MarkerLookup(pia=(1,), lv=(2,), am_u=(3,), am_l=(4,), ssas=(5, 8))
    F = advection_diffusion_residual(
        c, c_n, v, D, phi, eta, beta_pia, beta_exit, tau, g_influx, u, markers
    )
    J = ufl.derivative(F, c)

    c_ = ufl.TrialFunction(V)
    a, L = advection_diffusion_system(
        c_, c_n, v, D, phi, eta, beta_pia, beta_exit, tau, g_influx, u, markers
    )
    # If executed with FFCx
    form = [F, J, a, L]

    from mpi4py import MPI
    import dolfinx

    # Compile symbolic form
    F_compiled = dolfinx.fem.compile_form(MPI.COMM_WORLD, F)
    J_compiled = dolfinx.fem.compile_form(MPI.COMM_WORLD, J)
    a_compiled = dolfinx.fem.compile_form(MPI.COMM_WORLD, a)
    L_compiled = dolfinx.fem.compile_form(MPI.COMM_WORLD, L)

    # Create mesh and function spaces
    ct_form = dolfinx.mesh.to_type(domain.ufl_coordinate_element().cell_type.name)
    mesh = dolfinx.mesh.create_unit_cube(
        MPI.COMM_WORLD,
        10,
        10,
        10,
        cell_type=ct_form,
        ghost_mode=dolfinx.mesh.GhostMode.shared_facet,
    )
    Vh = dolfinx.fem.functionspace(mesh, V.ufl_element())
    ch = dolfinx.fem.Function(Vh)
    ch_n = dolfinx.fem.Function(Vh)

    Qh = dolfinx.fem.functionspace(mesh, Q.ufl_element())
    uh = dolfinx.fem.Function(Qh)

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

    F_compiled = dolfinx.fem.create_form

    vh = ufl.TestFunction(Vh)
    Fh = advection_diffusion_residual(
        ch, ch_n, vh, Dh, phih, etah, beta_piah, beta_exith, tauh, gh, uh, markers
    )

    Jh = ufl.derivative(Fh, ch)
    J_compiled = dolfinx.fem.form(Jh)
    F_compiled = dolfinx.fem.form(Fh)

    # Note to self:
    # Check idealized 3D geometry and refine and create boundary markers

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.