Tracer modelling within the brain parenchyma#
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}\]
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