Matrix Formalism

In this example we will consider the matrix formalism approach for a geometry of cylinders.

using SpinDoctor
using LinearAlgebra

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

setup = CylinderSetup(;
    name = "Slice",
    ncell = 3,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    height = 1.0,
    bend = 0.0,
    twist = 0.0,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.5,
)

CylinderSetup{Float64}("Slice", "meshfiles", 3, 2.0, 6.0, 0.2, 0.3, 1.0, 0.0, 0.0, false, 0.5, :convex_hull, 0.5, Inf)

We also define coefficients for the different cell compartments :in (axon), :out (myelin), and :ecs (ECS).

coeffs = coefficients(
    setup;
    D = (; in = 0.002 * I(3), out = 0.002 * I(3), ecs = 0.002 * I(3)),
    T₂ = (; in = Inf, out = Inf, ecs = Inf),
    ρ = (; in = 1.0, out = 1.0, ecs = 1.0),
    κ = (; in_out = 1e-4, out_ecs = 1e-4, in = 0.0, out = 0.0, ecs = 0.0),
    γ = 2.67513e-4,
)

(D = StaticArrays.SMatrix{3, 3, Float64, 9}[[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], [0.002 0.0 0.0; 0.0 0.002 0.0; 0.0 0.0 0.002]], T₂ = [Inf, Inf, Inf, Inf], κ = [0.0001, 0.0001, 0.0001, 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], γ = 0.000267513)

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

mesh, = 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 / 3 * tr.(model.D)' * volumes / sum(volumes)

0.0019999999999999996

The eigenfunctions of the diffusive part of the Bloch-Torrey operator forms a good basis for the finite element function space. The basis may be truncated at a certain level, thus reducing the number of degrees of freedom. We here opt for 400 eigenfunctions.

laplace = Laplace(; model, matrices, neig_max = 400)
lap_eig = solve(laplace)

(values = [1.6773224737553423e-17, 1.759114435161524e-5, 3.7638862735570044e-5, 6.606071974992256e-5, 7.556880675141598e-5, 9.120631518647759e-5, 0.00011126285453813121, 0.00014045626690437727, 0.00019895354582789173, 0.000233596615963599  …  0.048866033023565784, 0.04888862109002999, 0.04905801018170325, 0.04925786052489779, 0.049450853078994474, 0.04952356708075232, 0.04954963410934925, 0.04979109204567569, 0.04995678278514602, 0.05009861124463911], funcs = [-0.05414740335127055 -0.04434110354044968 … 0.01706784523712072 0.0008928209762885355; -0.054147403351270953 -0.043267017669457945 … -0.02645654110090974 -0.006177585959259527; … ; -0.054147403351271356 -0.01877646243887364 … 0.02621814093149775 -0.049940682018925965; -0.05414740335127108 -0.0176492471641877 … 0.013062983509291342 0.08193287287837207], moments = [[0.15879763282498116 -1.0764757389042474 … 0.0002855292107216134 -0.0016829918791405885; -1.0764757389042472 -0.07354419591814038 … 9.218075523152042e-5 0.0007077429928470354; … ; 0.00028552921072152234 9.21807552314597e-5 … 1.5568551151234813 2.0239241839847604; -0.0016829918791407238 0.0007077429928470009 … 2.023924183984761 0.894526125754191], [-0.030426231116839564 -6.332636513941051 … 0.0004067225641412528 0.00026489817300592855; -6.3326365139410505 0.3919674256426776 … 5.059094961742949e-5 -0.0006726882766968645; … ; 0.0004067225641413291 5.0590949617599495e-5 … 0.5851879230832477 -5.160851555890054; 0.00026489817300598406 -0.0006726882766968506 … -5.160851555890055 1.2585969756272495], [5.377642775528102e-17 5.941285029828058e-8 … 0.001221834266774143 -0.0008244861932605293; 5.941285030012372e-8 4.611212309741687e-6 … -0.0003471989596512768 0.0011960881123455603; … ; 0.001221834266774139 -0.0003471989596512746 … 0.01810220012812577 0.002105773914909766; -0.0008244861932605249 0.001196088112345561 … 0.002105773914909752 0.012284415395733835]], 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])

The resulting eigenvalues may be represented as length scales, describing the wavelength of the eigenfunctions.

length_scales = eig2length.(lap_eig.values, D_avg)

400-element Vector{Float64}:
  3.430493395544739e7
 33.49791870970913
 22.90058750831141
 17.28594870628898
 16.16194364876153
 14.71134618741523
 13.319556719340312
 11.854802138563507
  9.960680931828325
  9.192456788284169
  ⋮
  0.6354201408138161
  0.6343221925815088
  0.6330340900206038
  0.6317976068164046
  0.6313336114837056
  0.6311675239172678
  0.6296352656929225
  0.6285902486434186
  0.6276998529730011

We may also further truncate the eigenfunction basis, if we are satisfied skipping features below a threshold length scale of 3 micrometers.

length_scale = 3
λ_max = length2eig(length_scale, D_avg)
lap_eig = limit_lengthscale(lap_eig, λ_max)

(values = [1.6773224737553423e-17, 1.759114435161524e-5, 3.7638862735570044e-5, 6.606071974992256e-5, 7.556880675141598e-5, 9.120631518647759e-5, 0.00011126285453813121, 0.00014045626690437727, 0.00019895354582789173, 0.000233596615963599  …  0.0014264796969238225, 0.001436555759237195, 0.0014649571723772676, 0.0015333786940824007, 0.001582394118468939, 0.0016038248306292089, 0.0016978243493671167, 0.0017886138275346926, 0.0019086718579114136, 0.002107690096542497], funcs = [-0.05414740335127055 -0.04434110354044968 … 0.004981805801448434 -0.011666013905094549; -0.054147403351270953 -0.043267017669457945 … 0.0045757322955425825 -0.01592529066457598; … ; -0.054147403351271356 -0.01877646243887364 … 0.1493979765605434 -0.08529633951336474; -0.05414740335127108 -0.0176492471641877 … 0.1316741486057071 -0.08890689529630194], moments = [[0.15879763282498116 -1.0764757389042474 … 0.06971516278016571 0.0231186516795535; -1.0764757389042472 -0.07354419591814038 … 0.03306481125311574 0.02869016142665392; … ; 0.06971516278016573 0.03306481125311582 … -0.8567245342661187 1.8104505542075036; 0.02311865167955366 0.028690161426653923 … 1.8104505542075031 -0.6370199454455916], [-0.030426231116839564 -6.332636513941051 … -0.05221263090123049 -0.020334892072978998; -6.3326365139410505 0.3919674256426776 … -0.08660513680117332 0.007140888658539651; … ; -0.052212630901230406 -0.08660513680117336 … -1.2411870184501883 1.83073776046857; -0.020334892072978984 0.007140888658539429 … 1.8307377604685715 -0.10407820485917213], [5.377642775528102e-17 5.941285029828058e-8 … 1.9536692050002927e-6 7.275806889266764e-5; 5.941285030012372e-8 4.611212309741687e-6 … -5.943934188315687e-6 -1.783637699868533e-5; … ; 1.9536692050046295e-6 -5.94393418831482e-6 … -0.0007718153673362434 -0.00032135327100471663; 7.275806889267848e-5 -1.783637699868311e-5 … -0.0003213532710047151 0.0006471389116635435]], 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])

Each of the resulting eigenfunctions is represented in the same way as the initial magnetization field ρ.

ncompartment, nboundary = size(mesh.facets)
fig = Figure()
for i = 1:3, j = 1:4
    ieig = 6(i - 1) + j
    ϕ_cmpts = split_field(mesh, lap_eig.funcs[:, ieig])
    ax = Axis3(fig[i, j]; title = "n = $ieig, ℓ = $(length_scales[ieig])", aspect = :data)
    nboundary = size(mesh.facets, 2)
    scene = nothing
    first = true
    for icmpt = 1:ncompartment, iboundary = 1:nboundary
        facets = mesh.facets[icmpt, iboundary]
        points = mesh.points[icmpt]
        mesh!(ax, points', facets', color = ϕ_cmpts[icmpt], shading = false)
    end
end
fig

We observe that the first functions have large features, while the higher-index functions have more rapidly varying features. We may now choose a gradient and compute the projection of magnetization field onto the truncated basis.

dir = [1.0, 0.0, 0.0]
profile = CosOGSE(5000.0, 5000.0, 2)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)

ScalarGradient{Float64}([1.0, 0.0, 0.0], CosOGSE{Float64}(5000.0, 5000.0, 2), 4.201554156121445)

The matrix formalism problem is solved in the same way as the BTPDE. The time profile is approximated on 500 points, since it is non-constant.

mf = MatrixFormalism(; model, matrices, lap_eig)
ξ = solve(mf, gradient; ninterval = 500)

612-element Vector{ComplexF64}:
  0.2555561569186843 + 0.03386852374967236im
 0.25336658422734376 + 0.035001350964830785im
 0.24513672838387973 + 0.03289432246639865im
  0.2339777930908032 + 0.028502175302723892im
 0.22174366718554953 + 0.023022282066736444im
 0.20931735735923843 + 0.016864695738340194im
  0.2016035894740579 + 0.009748029751519437im
 0.19788543642817338 + 0.004451767714466655im
 0.19799063471165887 + 0.0002619500235386363im
 0.20041180563372668 - 0.004478761327529738im
                     ⋮
 0.21035523318093638 + 0.04383102262214298im
  0.2538203599292748 - 0.05170018639274829im
  0.2643807414855294 - 0.06091454430440761im
  0.2325791035287399 - 0.025516694829094373im
  0.2140130411614197 + 0.004209094177860869im
 0.23328467181740029 - 0.021772726130303524im
  0.2086884370257897 + 0.011752752710981187im
 0.20556556712361537 + 0.02746024371472481im
  0.2067940590206238 + 0.0318601324822539im

The resulting magnetization field may be plotted.

plot_field(model.mesh, ξ)

This page was generated using Literate.jl.