High angular resolution diffusion imaging

In this example we will compute the signal using the same gradient sequence in many different directions.

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 vertical direction, but a rather restricted vertical diffusion with the permeable membranes.

ncell = 5
setup = PlateSetup(; depth = 50.0, widths = fill(5.0, ncell), refinement = 1.0)
coeffs = coefficients(
    setup;
    D = [0.002 * I(2) for d ∈ 1:ncell],
    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.002], [0.002 0.0; 0.0 0.002], [0.002 0.0; 0.0 0.002], [0.002 0.0; 0.0 0.002], [0.002 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, surfaces, cells = create_geometry(setup; recreate = true)
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 / 2 * tr.(model.D)' * volumes / sum(volumes)
ncompartment = length(model.mesh.points)

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

directions = unitcircle(100)[1:2, :]
profile = PGSE(2500.0, 4000.0)
b = 1000
g = √(b / int_F²(profile)) / model.γ
gradients = [ScalarGradient(d, profile, g) for d ∈ eachcol(directions)]

100-element Vector{ScalarGradient{Float64, PGSE{Float64}}}:
 ScalarGradient{Float64, PGSE{Float64}}([1.0, 0.0], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9980267284282716, 0.06279051952931337], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9921147013144779, 0.12533323356430426], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9822872507286887, 0.1873813145857246], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9685831611286311, 0.2486898871648548], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9510565162951536, 0.30901699437494745], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9297764858882515, 0.3681245526846779], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9048270524660196, 0.42577929156507266], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.8763066800438636, 0.4817536741017153], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.8443279255020151, 0.5358267949789967], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ⋮
 ScalarGradient{Float64, PGSE{Float64}}([0.8443279255020147, -0.5358267949789971], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.8763066800438631, -0.4817536741017161], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9048270524660194, -0.425779291565073], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9297764858882515, -0.36812455268467786], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9510565162951535, -0.3090169943749476], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.968583161128631, -0.24868988716485535], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9822872507286887, -0.18738131458572468], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9921147013144778, -0.12533323356430465], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)
 ScalarGradient{Float64, PGSE{Float64}}([0.9980267284282716, -0.06279051952931326], PGSE{Float64}(2500.0, 4000.0), 0.8402604538082508)

The signals are computed from the magnetization field through quadrature.

ρ = initial_conditions(model)
S₀ = abs(compute_signal(matrices.M, ρ))

1250.0

We may solve the BTPDE for each gradient.

btpde = BTPDE(; model, matrices)
solver = IntervalConstantSolver(; timestep = 10.0)
signals = map(gradients) do grad
    @show grad.dir
    ξ = solve(btpde, grad, solver)
    abs(compute_signal(matrices.M, ξ))
end

100-element Vector{Float64}:
 659.8325751201447
 656.64760747521
 647.2207604628692
 631.9798267036172
 611.597730765619
 586.9372558523728
 558.9837208152493
 528.7735991782711
 497.32684053762375
 465.5893004202092
   ⋮
 465.5139504207484
 497.24928418848685
 528.6957751245288
 558.908103422209
 586.8668075279813
 611.5357395658543
 631.9296247727546
 647.1853542507154
 656.6292850001244

We may plot the directionalized signal attenuations.

attenuations = signals ./ S₀
plot_hardi(directions, attenuations)

The signal attenuates the most in the vertical direction, as that is where diffusion is restricted the least.

This page was generated using Literate.jl.