Edit on GitHub

Compare ADCs

The apparent diffusion coefficient (ADC) may be sufficient to describe the signal attenuation for small $b$-values. SpinDoctor comes with multiple approaches for computing or estimating the ADC for a ScalarGradient $\vec{g}(t) = f(t) g \vec{d}$:

  • Using the free diffusion coefficient $\vec{d}^\mathsf{T} \mathbf{D} \vec{d}$, which represents unrestricted diffusion in the absence of boundaries;
  • Computing the short diffusion time approximation for the ADC
  • Fitting the log-signal obtained by solving the BTPDE for different $b$-values
  • Solving a homogenized model (HADC) assuming negligible permeability
  • Using the matrix formalism effective diffusion tensor

In this example we will compare the different approaches.

Building a biological model

We start by loading SpinDoctor and a Makie plotting backend.

using SpinDoctor
using LinearAlgebra

if haskey(ENV, "GITHUB_ACTIONS")
    using CairoMakie
else
    using GLMakie
end

Here we create a recipe for five stacked plates with isotropic diffusion tensors. They should allow for free diffusion in the horizontal direction, but a rather restricted vertical diffusion with the permeable membranes.

ncell = 5
setup = SlabSetup(;
    depth = 50.0,
    widths = fill(5.0, ncell),
    height = 50.0,
    bend = 0.0,
    twist = π / 6,
    refinement = 10.0,
)
d = fill(0.002, ncell)
d[2] = 0.005
coeffs = coefficients(
    setup;
    D = [d * I(3) for d ∈ d],
    T₂ = fill(Inf, ncell),
    ρ = fill(1.0, ncell),
    κ = (; interfaces = fill(1e-4, ncell - 1), boundaries = fill(0.0, ncell)),
    γ = 2.67513e-4,
)
(D = [[0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002], [0.005 0.0 0.0; 0.0 0.005 0.0; 0.0 0.0 0.005], [0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002], [0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002], [0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002]], T₂ = [Inf, Inf, Inf, Inf, Inf], κ = [0.0001, 0.0, 0.0, 0.0, 0.0001, 0.0, 0.0, 0.0001, 0.0, 0.0001, 0.0, 0.0, 0.0, 0.0, 0.0], ρ = ComplexF64[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im], γ = 0.000267513)

We then proceed to build the geometry and finite element mesh.

mesh, = create_geometry(setup)
plot_mesh(mesh)

The mesh looks good, so we may then proceed to assemble the biological model and the associated finite element matrices.

model = Model(; mesh, coeffs...)
matrices = assemble_matrices(model);

We may also compute some useful quantities, including a scalar diffusion coefficient from the diffusion tensors.

volumes = get_cmpt_volumes(model.mesh)
D_avg = 1 / 3 * tr.(model.D)' * volumes / sum(volumes)
ncompartment = length(model.mesh.points)
5

The gradient pulse sequence will be a PGSE with both x and z components. This allows for both restricted horizontal diffusion and almost unrestricted vertical diffusion. The different approaches should hopefully confirm this behaviour.

dir = [1.0, 0.0, 1.0]
profile = PGSE(2500.0, 4000.0)
b = 1000
g = √(b / int_F²(profile)) / model.γ
gradient = ScalarGradient(dir, profile, g)
ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865475, 0.0, 0.7071067811865475], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)

Computing the ADC using different methods

Short term approximation

A simple model for the ADC is given by the short term approximation.

adc_sta_cmpts = compute_adc_sta(model, gradient)
adc_sta = volumes'adc_sta_cmpts / sum(volumes)
0.0008978720553849736

Fitting the ADC

A more robust approach is to directly fit the BTPDE signal to a series of b-values. This is however more computationally expensive. We start by building a set of gradients.

bvalues = 0:400:4000
gvalues = map(b -> √(b / int_F²(profile)) / coeffs.γ, bvalues)
gradients = [ScalarGradient(gradient.dir, gradient.profile, g) for g ∈ gvalues]
11-element Vector{ScalarGradient{Float64, PGSE{Float64}}}:
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 0.0)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 0.5314273723601551)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 0.7515517974080282)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 0.9204592094606131)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.0628547447203103)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.1883077297013998)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.3017258976304167)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.4060246671574907)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.5031035948160565)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.5942821170804655)
 ScalarGradient{Float64, PGSE{Float64}}([0.7071067811865476, 0.0, 0.7071067811865476], PGSE{Float64}(2500.0, 4000.0), 1.6805209076165015)

We compute the magnetization for each gradient. Since the gradients are interval-wise constant, we can use a specialized solver.

btpde = BTPDE(; model, matrices)
solver = IntervalConstantSolver(; θ = 0.5, timestep = 5.0)
signals_cmpts = map(gradients) do grad
    @show grad.g
    ξ = solve(btpde, grad, solver)
    compute_signal.(matrices.M_cmpts, split_field(model.mesh, ξ))
end
signals = sum(signals_cmpts)

Fitting the ADC is straightforward.

adc_fit = fit_adc(bvalues, signals)
adc_fit_cmpts =
    [fit_adc(bvalues, [s[icmpt] for s ∈ signals_cmpts]) for icmpt = 1:ncompartment]

Homogenized ADC model

The HADC-model uses homogenization and assumes negligible permeability between the compartments. This does require solving an ODE involving all the degrees of freedom.

hadc = HADC(; model, matrices)
adc_homogenized_cmpts = solve(hadc, gradient)
adc_homogenized = volumes'adc_homogenized_cmpts / sum(volumes)
0.0018624675877559336

Matrix Formalism

If a Laplace eigendecomposition is available, the HADC can be approximated with little additional computational expense.

laplace = Laplace(; model, matrices, neig_max = 400)
lap_eig = solve(laplace)
(values = [-4.781790432054211e-17, 6.472281898550158e-6, 9.2733102574735e-6, 9.559958685849962e-6, 1.596174670129549e-5, 1.6477449715310355e-5, 1.7801179137596793e-5, 2.4454070190484735e-5, 2.5878161999050007e-5, 3.2704519471457126e-5  …  0.0009374097288678095, 0.0009421200928268779, 0.0009428320468513823, 0.0009474386348285552, 0.0009484197416795614, 0.0009503782788641251, 0.000952113696274504, 0.0009538405221474771, 0.000954704625622591, 0.000956088409193153], funcs = [-0.004000158043524799 0.005542254415819706 … -6.262240444875598e-5 -0.02605767648260966; -0.004000158043523761 0.00511913887385148 … -0.0001041168483042911 -0.010840664813136574; … ; -0.004000158043523739 -0.0054646531934355504 … -0.007470110774208972 0.00021060496995797272; -0.004000158043523913 -0.005487169740022619 … -0.0027402041678590418 -4.7780044549801375e-5], moments = [[0.000222082165838966 7.057787537559466 … 0.006305959896651989 0.007452511718661784; 7.057787537559467 0.26691894812817046 … 0.006548810382687481 -0.010033275548614738; … ; 0.006305959896651926 0.006548810382687439 … 4.362310562567342 0.038544066351455485; 0.007452511718661684 -0.010033275548614717 … 0.03854406635145549 -10.20802236623691], [-0.00010093363309052195 -0.0008207669586151045 … -0.0029284004526531557 0.0058737033325518846; -0.0008207669586155764 -0.0011456707505955732 … -0.0035427684863816172 -0.005725625871294609; … ; -0.002928400452653024 -0.0035427684863815964 … -0.42765697033658234 0.019782828453520734; 0.005873703332551855 -0.005725625871294834 … 0.019782828453520744 -0.8800170207559181], [-3.8418825013708524e-5 0.0011205010513714864 … 0.0008169278529724794 9.595864454528666e-5; 0.0011205010513714864 2.102670009218599e-5 … -0.0007404798219458886 0.0001351497116467192; … ; 0.0008169278529725627 -0.0007404798219461661 … -1.859804747304655 -0.025730285388067584; 9.595864454525153e-5 0.0001351497116464475 … -0.025730285388067594 9.068799985860611]], massrelax = [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])

An effective diffusion tensor can be computed.

D_mf = compute_mf_diffusion_tensor(model.mesh, matrices.M, lap_eig, gradient)
3×3 Matrix{Float64}:
  0.000420604  -2.78279e-8   2.98397e-8
 -2.78279e-8    0.00225041  -4.95e-6
  2.98397e-8   -4.95e-6      0.00225097

We observe that $D_{x x}$ and $D_{y y}$ are very close to the intrinsic diffusion coefficient $D = 0.002$, which confirms that diffusion in the horizontal direction is almost unrestricted. $D_{z z}$ is significantly smaller, confirming the presence of membranes along the vertical direction.

In particular, we may deduce the MF-ADC in our direction.

adc_mf = dir'D_mf * dir / dir'dir
0.001335817329315722

Comparing results

Compartment ADCs

Here we make a comparison between the compartment ADCs. For MF, only a global ADC was computed.

n = ncompartment
fig = Figure()
ax = Axis(
    fig[1, 1];
    xticks = (1:4, ["STA", "Fit BTPDE", "HADC", "MF"]),
    ylabel = "ADC / D",
    title = "Compartment ADCs",
)
barplot!(ax, fill(1, n), adc_sta_cmpts ./ D_avg; dodge = 1:n)
barplot!(ax, fill(2, n), adc_fit_cmpts ./ D_avg; dodge = 1:n)
barplot!(ax, fill(3, n), adc_homogenized_cmpts ./ D_avg; dodge = 1:n)
barplot!(ax, [4], [adc_mf / D_avg])
fig

The STA and HADC are the same for all compartments, as they consider them separately and all the compartments have the same size. The fitted ADC is larger for the three inner compartments, as they all have permeable membranes both below and above, in contrast to the top and bottom compartments that have hard walls.

Signal attenuation

We may also inspect the resulting signal attenuations.

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "b", yscale = log10, title = "Signal attenuation")
lines!(
    ax,
    [0, bvalues[end]],
    [1, exp(-adc_sta * bvalues[end])];
    linestyle = :dash,
    label = "ADC STA",
)
lines!(
    ax,
    [0, bvalues[end]],
    [1, exp(-adc_fit * bvalues[end])];
    linestyle = :dash,
    label = "ADC Fit",
)
lines!(
    ax,
    [0, bvalues[end]],
    [1, exp(-adc_homogenized * bvalues[end])];
    linestyle = :dash,
    label = "HADC",
)
lines!(
    ax,
    [0, bvalues[end]],
    [1, exp(-adc_mf * bvalues[end])];
    linestyle = :dash,
    label = "ADC MF",
)
lines!(
    ax,
    [0, bvalues[end]],
    [1, exp(-D_avg * bvalues[end])];
    linestyle = :dash,
    label = "Free diffusion",
)
scatterlines!(ax, bvalues, abs.(signals) ./ abs(signals[1]); label = "BTPDE Signal")
axislegend(ax)
fig

We observe the following:

  • The exact signal starts to visually deviate from the log-linear regime after $b = 2000$. For higher $b$-values, the ADC is no longer sufficient to describe the attenuation, as higher order terms can no longer be neglected.
  • The fit-ADC signal coincides with the exact signal for the lowest $b$-values. This makes sense since that is how this ADC was obtained to begin with. This is considered to be the "reference" ADC.
  • The free diffusion signal attenuates more thant the exact signal by many orders of magnitude, which confirms the presence of restrictive membranes and boundaries in the gradient direction.
  • The HADC signal attenuates less than the exact signal, as it assumes a more severe restriction with impermeable membranes.
  • The MF-ADC signal coincides with the HADC signal, as the former is simply an MF approximation of the latter.

This page was generated using Literate.jl.