Edit on GitHub

Solve BTPDE

We start by loading SpinDoctor and a Makie plotting backend.

using SpinDoctor
using LinearAlgebra
using Random

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

Set random seed for reproducibility

Random.seed!(123)
Random.TaskLocalRNG()

The built in geometry recipes allow for making various cell configurations. We here consider the case of three twisted axons immersed in an extracellular space (ECS). The axons are extruded from a 2D ground setup consisting of three random disks.

setup = CylinderSetup(;
    ncell = 3,
    ecs = ConvexHullECS(; margin = 2.0),
    height = 40.0,
    bend = 0.0,
    twist = π / 4,
)
CylinderSetup{Float64, ConvexHullECS{Float64}}
  groundsetup: DiskSetup{Float64, ConvexHullECS{Float64}}
  refinement: Float64 Inf
  height: Float64 40.0
  bend: Float64 0.0
  twist: Float64 0.7853981633974483

We also define parameters for the cell compartments and ECS. The required parameters are:

  • 3D Diffusion tensors D (we here use an isotropic coefficient)
  • T2-relaxation times T₂
  • Initial spin densities ρ
  • Permeabilities κ

The parameters D, T₂, and ρ are defined for each compartment, while κ is defined for each interface and outer boundary. The DiskSetup allows for multiple layers in each cell (here nlayer = 1), so cell and cell_boundaries are arrays of length nlayer = 1, while cell_interfaces is of length nlayer - 1 = 0.

coeffs = coefficients(
    setup;
    D = (; cell = [0.002 * I(3)], ecs = 0.002 * I(3)),
    T₂ = (; cell = [Inf], out = Inf, ecs = Inf),
    ρ = (; cell = [1.0], out = 1.0, ecs = 1.0),
    κ = (; cell_interfaces = zeros(0), cell_boundaries = [0.0], cell_ecs = 1e-4, ecs = 0.0),
    γ = 2.67513e-4,
)
(D = [[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.0, 0.0, 0.0001, 0.0, 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)

The following line creates a random cell configuration for our cylinder setup, generates a surface triangulation and calls TetGen to create a tetrahedral finite element mesh. The compartments and boundaries will be ordered in the same way as coeffs.

mesh, surfaces, cells = create_geometry(setup)
(SpinDoctor.FEMesh{Float64, 3}([[1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  1061, 1062, 1064, 1066, 1102, 1106, 1107, 1111, 1113, 1118], [13, 14, 15, 16, 17, 18, 19, 20, 21, 22  …  1050, 1051, 1097, 1098, 1099, 1100, 1103, 1104, 1105, 1109], [30, 31, 32, 33, 34, 35, 36, 37, 38, 39  …  1038, 1039, 1094, 1095, 1096, 1101, 1108, 1110, 1112, 1114], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  1089, 1090, 1091, 1092, 1093, 1115, 1116, 1117, 1119, 1120]], [[15.082718079328439 15.352979394778533 … 8.951502240400457 10.805043881484655; -0.877699437546676 1.1751390607938879 … 7.048908264966403 3.6467824945986953; -20.0 -20.0 … 11.25 -6.25], [5.566461744061371 5.98781998613544 … 0.43441965610525357 -0.02563698116508828; 2.0644025213185593 4.0618578866720005 … 4.190163956819523 4.204390404839309; -20.0 -20.0 … -20.0 -14.414328232678379], [-6.038998579123932 -5.671742422014201 … -9.097074573563216 -10.516569589158443; -2.5745701021829923 -0.5930331991229953 … -3.5277112181337125 -0.7198994633381002; -20.0 -20.0 … 2.5 -20.0], [15.082718079328439 15.352979394778533 … 5.509007211397983 3.954806529363884; -0.877699437546676 1.1751390607938879 … -4.4931455537572536 -5.88647187996765; -20.0 -20.0 … -1.7739273663246702 -16.90175538468711]], [Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0) … Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0); Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0) … Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0); Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0) … [12 9 … 15 14; 13 10 … 1 15; 248 248 … 248 248] Matrix{Int64}(undef, 3, 0); Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0) … Matrix{Int64}(undef, 3, 0) [107 757 … 1090 598; 108 455 … 719 592; 825 460 … 1093 1098]], [[95 117 … 53 179; 166 66 … 55 56; 91 120 … 115 180; 198 196 … 197 197], [215 152 … 260 107; 33 151 … 258 256; 34 253 … 279 258; 277 275 … 280 280], [82 142 … 172 168; 146 64 … 49 111; 90 144 … 170 170; 242 242 … 246 246], [153 409 … 211 711; 161 408 … 719 712; 332 410 … 712 719; 342 1024 … 1093 1093]]), (points = [14.27876664852878 13.742868263666535 … -7.419914275384653 -15.047497514861513; 4.932903284952767 6.932903284952767 … 8.252873729197479 -0.8063742161941281; -20.0 -20.0 … 20.0 20.0], facets = [2 3 … 62 86; 1 2 … 63 88; 62 63 … 65 85], facetmarkers = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3  …  7, 10, 10, 10, 10, 10, 7, 10, 7, 8], regions = [10.27876664852878 -1.1938838093179314 -9.432278711012088 -15.151006859813176; 4.932903284952767 4.00933742966369 -4.717731819722998 -3.5021786272882673; 0.0 0.0 0.0 0.0]), (radii = [4.0 5.554882327575999 4.846487937516692], centers = [10.27876664852878 -1.1938838093179314 -9.432278711012088; 4.932903284952767 4.00933742966369 -4.717731819722998]))

The resulting mesh can be plotted in 3D provided a Makie backend is loaded.

plot_mesh(mesh)

The mesh looks good, so we can proceed with the assembly our biological model and the associated finite element matrices.

model = Model(; mesh, coeffs...)
matrices = assemble_matrices(model)
(M = sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  1145, 1150, 1318, 1323, 1324, 1768, 1772, 1782, 1821, 1824], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824], [1.3325303041367822, 0.16666666666666655, 0.16666666666666696, 0.3309413020824443, 0.3329318187350576, 0.3353238499859468, 0.6662651520683911, 0.16666666666666655, 0.9952159374982228, 0.16666666666666682  …  2.670540222315831, 3.3953556696452445, 1.612635361912306, 2.9084018282843815, 1.488906763332342, 2.135633032553417, 2.468178934200736, 2.7650446245303173, 2.1515222503620994, 20.608948102657422], 1824, 1824), S = sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  1145, 1150, 1318, 1323, 1324, 1768, 1772, 1782, 1821, 1824], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824], [0.006346996905349302, -0.001977075047162113, -0.0012664657614322714, 0.0010338376704796385, -0.002492918837158337, -0.0004164220941403156, -0.0012279528359359026, -0.001977075047162113, 0.006988755547910822, -0.001925260523849295  …  -0.0014067398416301178, -0.0015334181696860757, -0.0029915845959078945, -0.007377062330059798, 0.0019523830483469686, -0.0029824290085358656, -0.006122695920547249, -0.0003349779169020211, -0.0015176138930550226, 0.03059436456985027], 1824, 1824), R = sparse(Int64[], Int64[], Float64[], 1824, 1824), Mx = SparseArrays.SparseMatrixCSC{Float64, Int64}[sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  1145, 1150, 1318, 1323, 1324, 1768, 1772, 1782, 1821, 1824], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824], [19.11782469042412, 2.4315511144888333, 2.3090781458085727, 4.817320599469516, 4.792182159037534, 4.687881386917625, 9.148546520397721, 2.4315511144888333, 14.448326688036559, 2.377954315993194  …  10.835120767317001, 20.585324462159416, 11.206866010760761, 10.83063544590817, 9.898208990797107, 3.7215582199460324, 9.24950870139131, 16.445344590012972, 2.8820119886971813, 82.40474566171758], 1824, 1824), sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  1145, 1150, 1318, 1323, 1324, 1768, 1772, 1782, 1821, 1824], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824], [-0.5238461789627057, 0.08809193431713092, -0.22183063570571326, 0.1801774333145023, -0.07586995498830093, -0.3415365863268203, -0.09194399930883262, 0.08809193431713092, 1.446710561093754, 0.36040626920666974  …  -7.766840671243365, -11.826034216416915, -7.2548123002154785, -13.604974041151529, -6.642380978030072, -8.20066419191971, -7.120506400629633, -9.083852490187384, -11.293568690947072, -91.53780159658153], 1824, 1824), sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  1145, 1150, 1318, 1323, 1324, 1768, 1772, 1782, 1821, 1824], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824, 1824], [-25.817941940399432, -3.2638888888888866, -3.263888888888895, -6.274593858490276, -6.242471601282329, -6.356766631680946, -12.908970970199716, -3.2638888888888866, -19.35274991316038, -3.263888888888892  …  -42.11702453218053, -52.25326282429954, -26.27097249339122, -43.172832032195394, -21.47194304168049, -36.984529280208314, -42.99524947260839, -47.27191317004768, -35.26826469011216, -344.4675572255478], 1824, 1824)], Q = sparse([1, 2, 12, 190, 191, 192, 727, 728, 738, 1788  …  192, 727, 737, 738, 1167, 1171, 1173, 1787, 1790, 1792], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1792, 1792, 1792, 1792, 1792, 1792, 1792, 1792, 1792, 1792], [0.00017356674032594226, 2.1619686615169436e-5, 2.1752261921656375e-5, 4.3190057830009714e-5, 4.341142162614532e-5, 4.3593312332961416e-5, -0.00017356674032594226, -2.1619686615169436e-5, -2.1752261921656375e-5, -4.3190057830009714e-5  …  -0.000350737092314327, 4.3593312332961416e-5, 4.408887426484085e-5, 4.3902953741119305e-5, 4.4082861744488555e-5, 4.3603497814872674e-5, 4.383546592103628e-5, 4.390578196160991e-5, 4.3724344533398e-5, 0.000350737092314327], 1824, 1824), M_cmpts = Any[sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  161, 162, 163, 164, 165, 166, 167, 168, 194, 198], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  198, 198, 198, 198, 198, 198, 198, 198, 198, 198], [1.3325303041367822, 0.16666666666666655, 0.16666666666666696, 0.3309413020824443, 0.3329318187350576, 0.3353238499859468, 0.6662651520683911, 0.16666666666666655, 0.9952159374982228, 0.16666666666666682  …  0.7547365945236499, 1.232536120202826, 2.966019587718394, 2.582828728340802, 0.4004034766974629, 0.7342348779683352, 0.9179639098339009, 1.2665265473867844, 3.9999999999999996, 37.985545474462064], 198, 198), sparse([1, 2, 17, 269, 271, 272, 279, 1, 2, 3  …  266, 267, 268, 269, 270, 271, 272, 276, 279, 280], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  280, 280, 280, 280, 280, 280, 280, 280, 280, 280], [1.856669790307353, 0.23222358552789912, 0.23222358552789926, 0.466297480790817, 0.46388772409787815, 0.4620374143628595, 0.9283348951536765, 0.23222358552789912, 1.388521999781517, 0.23222358552789926  …  1.5026399029263158, 1.9648101036019774, 1.9670432323798288, 1.5006428584523546, 1.7390521894018514, 1.9656014042348289, 1.5015281820682624, 8.944656160236484, 8.820448132647824, 67.09359502099363], 280, 280), sparse([1, 2, 15, 229, 231, 232, 248, 1, 2, 3  …  233, 234, 235, 236, 237, 238, 239, 240, 244, 248], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  248, 248, 248, 248, 248, 248, 248, 248, 248, 248], [1.5913095887190005, 0.19903357091627982, 0.19903357091627996, 0.39990324888555523, 0.3975876525269405, 0.39575154547394503, 0.7956547943595003, 0.19903357091627982, 1.1895702327804494, 0.19903357091627982  …  1.5915259938517545, 1.5839633254194139, 1.591443329928485, 1.5848577906725785, 1.5907025817551501, 1.5862328712407718, 1.5878508032015364, 1.5894318314204066, 7.463758909360495, 26.816670996585767], 248, 248), sparse([1, 2, 12, 53, 1062, 1064, 1066, 1068, 1, 2  …  419, 424, 592, 597, 598, 1042, 1046, 1056, 1095, 1098], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098], [0.9995984854017231, 0.17204447289647445, 0.16108810313772065, 0.16666666666666652, 0.17204447289647445, 0.16666666666666646, 0.16108810313772062, 0.49979924270086157, 0.17204447289647445, 2.321373177078888  …  2.670540222315831, 3.3953556696452445, 1.612635361912306, 2.9084018282843815, 1.488906763332342, 2.135633032553417, 2.468178934200736, 2.7650446245303173, 2.1515222503620994, 20.608948102657422], 1098, 1098)], S_cmpts = Any[sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  161, 162, 163, 164, 165, 166, 167, 168, 194, 198], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  198, 198, 198, 198, 198, 198, 198, 198, 198, 198], [0.006346996905349302, -0.001977075047162113, -0.0012664657614322714, 0.0010338376704796385, -0.002492918837158337, -0.0004164220941403156, -0.0012279528359359026, -0.001977075047162113, 0.006988755547910822, -0.001925260523849295  …  0.0017495754748501349, -0.0007469995388923757, -0.004603240360876724, -0.004645195832269418, 0.002520491675485374, 0.0007315803964535931, 0.0007739644902694654, -0.0005355655805505135, -0.004340678201397818, 0.04860733979632886], 198, 198), sparse([1, 2, 17, 269, 271, 272, 279, 1, 2, 3  …  266, 267, 268, 269, 270, 271, 272, 276, 279, 280], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  280, 280, 280, 280, 280, 280, 280, 280, 280, 280], [0.00804173575255135, -0.0023778856978062665, -0.002057433198595067, -0.00013552793005669348, -0.0031809922333407206, 0.0005542470420892754, -0.0008441437348418784, -0.0023778856978062665, 0.008141834554975774, -0.002447425349101389  …  0.00025261292866535364, -0.0003733555716873339, -0.0003348683808151821, 0.00032045889182093576, 3.283365778210914e-5, -0.0002818141338107674, 0.0003725996480014941, -0.008457233992183353, -0.00837254901151355, 0.05139591384927838], 280, 280), sparse([1, 2, 15, 229, 231, 232, 248, 1, 2, 3  …  233, 234, 235, 236, 237, 238, 239, 240, 244, 248], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  248, 248, 248, 248, 248, 248, 248, 248, 248, 248], [0.007623699661914772, -0.0017121457057828649, -0.001932770481331576, 0.0003414295743444949, -0.003021449198091399, -0.00011132791627018505, -0.0011874359347832418, -0.0017121457057828649, 0.006687998426107217, -0.0017344285103734613  …  -0.0003979889800367365, -0.0007135636309145903, -0.0006307732611122894, -0.0009346644872916034, -0.0008676877801289255, -0.0010894313015867045, -0.001133326763902845, -0.001051885310977223, -0.0043853320925054656, 0.02344190991931171], 248, 248), sparse([1, 2, 12, 53, 1062, 1064, 1066, 1068, 1, 2  …  419, 424, 592, 597, 598, 1042, 1046, 1056, 1095, 1098], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098], [0.010242327857958476, -0.003355306710475011, -0.0033943189510243707, -0.0021682054572781098, 0.0013278427244182197, -0.0009468765864350765, 0.0008371884826219335, -0.0025426513597860624, -0.003355306710475011, 0.013296522866762226  …  -0.0014067398416301178, -0.0015334181696860757, -0.0029915845959078945, -0.007377062330059798, 0.0019523830483469686, -0.0029824290085358656, -0.006122695920547249, -0.0003349779169020211, -0.0015176138930550226, 0.03059436456985027], 1098, 1098)], Mx_cmpts = Vector{Any}[[sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  161, 162, 163, 164, 165, 166, 167, 168, 194, 198], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  198, 198, 198, 198, 198, 198, 198, 198, 198, 198], [19.11782469042412, 2.4315511144888333, 2.3090781458085727, 4.817320599469516, 4.792182159037534, 4.687881386917625, 9.148546520397721, 2.4315511144888333, 14.448326688036559, 2.377954315993194  …  7.007419414814614, 12.553223682350703, 32.10685287663737, 28.747533276483406, 4.825369118793366, 9.367994258122518, 11.738111927021297, 15.948581427710725, 43.95184979899075, 410.53266455483634], 198, 198), sparse([1, 2, 17, 269, 271, 272, 279, 1, 2, 3  …  266, 267, 268, 269, 270, 271, 272, 276, 279, 280], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  280, 280, 280, 280, 280, 280, 280, 280, 280, 280], [8.533169268326766, 1.1349813603772396, 0.963766866555322, 1.9727597261201197, 2.132692793055948, 2.2220547691058554, 3.472542341888884, 1.1349813603772396, 6.841727009077427, 1.1256256092444106  …  2.1192932804716764, 3.210698636344423, 4.681001878260988, 3.6559670376193596, 5.246214895335581, 6.094992883665842, 4.7390224555183496, -2.0873031022928594, 1.8746404429431225, -0.4401188362468291], 280, 280), sparse([1, 2, 15, 229, 231, 232, 248, 1, 2, 3  …  233, 234, 235, 236, 237, 238, 239, 240, 244, 248], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  248, 248, 248, 248, 248, 248, 248, 248, 248, 248], [-10.901057170477387, -1.3127406084861089, -1.4568951306035158, -2.8700730572313238, -2.6976942568995215, -2.615457723015599, -6.044295392108538, -1.3127406084861089, -7.828137790057244, -1.3383645043011105  …  -20.751341182493345, -12.832742466604351, -20.815415635318733, -13.916190619145512, -20.15487840120918, -15.470279189236418, -17.227818064757, -18.885241144434886, -77.97690231739713, -281.2264624305779], 248, 248), sparse([1, 2, 12, 53, 1062, 1064, 1066, 1068, 1, 2  …  419, 424, 592, 597, 598, 1042, 1046, 1056, 1095, 1098], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098], [15.395721402791153, 2.6936585388704195, 2.4239220604158946, 2.6556980158027166, 2.66512340554267, 2.5548917274920298, 2.4029871772534968, 7.856830806286128, 2.6936585388704195, 36.72042994086794  …  10.835120767317001, 20.585324462159416, 11.206866010760761, 10.83063544590817, 9.898208990797107, 3.7215582199460324, 9.24950870139131, 16.445344590012972, 2.8820119886971813, 82.40474566171758], 1098, 1098)], [sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  161, 162, 163, 164, 165, 166, 167, 168, 194, 198], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  198, 198, 198, 198, 198, 198, 198, 198, 198, 198], [-0.5238461789627057, 0.08809193431713092, -0.22183063570571326, 0.1801774333145023, -0.07586995498830093, -0.3415365863268203, -0.09194399930883262, 0.08809193431713092, 1.446710561093754, 0.36040626920666974  …  4.644070572409429, 2.3271125075379713, 16.150474264924107, 8.327024401521001, 2.4944674701569705, 1.9452492864944464, 4.837896151303015, 5.295350347824675, 12.097611172671005, 137.30184703861548], 198, 198), sparse([1, 2, 17, 269, 271, 272, 279, 1, 2, 3  …  266, 267, 268, 269, 270, 271, 272, 276, 279, 280], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  280, 280, 280, 280, 280, 280, 280, 280, 280, 280], [4.673265681138138, 0.8046468253320737, 0.37138194637761246, 0.8395230019003259, 1.1923575908210409, 1.5566291725349846, 2.6655359271269883, 0.8046468253320737, 6.066322520795689, 1.1924876435629637  …  10.410417948199248, 2.8401616897023176, 12.581329157356825, 3.376903059462519, 9.268427358869912, 6.237338737421015, 6.44937007037304, 37.48480347551342, 37.04718770163679, 281.8955624245004], 280, 280), sparse([1, 2, 15, 229, 231, 232, 248, 1, 2, 3  …  233, 234, 235, 236, 237, 238, 239, 240, 244, 248], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  248, 248, 248, 248, 248, 248, 248, 248, 248, 248], [-3.6669842844513165, -0.26303568805926875, -0.6278242651343988, -1.2391640017779517, -0.9516668917860182, -0.6162627989952668, -1.5875458779332023, -0.26303568805926875, -0.5511170388513328, 0.06259752452119742  …  -2.783455814084385, 0.09707940339877752, -1.0356446298192223, 1.467820574423612, 0.5874049073440297, 2.2839519628409204, 2.403424913716571, 1.8038237764398675, -9.878914413531936, -26.68098100225745], 248, 248), sparse([1, 2, 12, 53, 1062, 1064, 1066, 1068, 1, 2  …  419, 424, 592, 597, 598, 1042, 1046, 1056, 1095, 1098], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098], [-0.7299485773910855, -0.000665974639365581, -0.25101909725502464, -0.181337832009083, 0.05205478126986286, -0.07804341389543414, -0.20254309843983243, -0.35966900560849335, -0.000665974639365581, 3.6740370963571367  …  -7.766840671243365, -11.826034216416915, -7.2548123002154785, -13.604974041151529, -6.642380978030072, -8.20066419191971, -7.120506400629633, -9.083852490187384, -11.293568690947072, -91.53780159658153], 1098, 1098)], [sparse([1, 2, 12, 190, 191, 192, 195, 1, 2, 3  …  161, 162, 163, 164, 165, 166, 167, 168, 194, 198], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  198, 198, 198, 198, 198, 198, 198, 198, 198, 198], [-25.817941940399432, -3.2638888888888866, -3.263888888888895, -6.274593858490276, -6.242471601282329, -6.356766631680946, -12.908970970199716, -3.2638888888888866, -19.35274991316038, -3.263888888888892  …  -3.1095397180724413, -5.038137845474553, -11.539612094814075, -9.824286368086675, -1.6404364453611868, -3.1705706964164184, -3.8600101889903584, -4.8884386549898435, -36.666666666666664, -242.8117173890429], 198, 198), sparse([1, 2, 17, 269, 271, 272, 279, 1, 2, 3  …  266, 267, 268, 269, 270, 271, 272, 276, 279, 280], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  280, 280, 280, 280, 280, 280, 280, 280, 280, 280], [-35.973210290104106, -4.547711883254691, -4.547711883254694, -8.839837592131111, -8.697894826835213, -8.759961346606907, -17.986605145052046, -4.547711883254691, -27.000377638358913, -4.547711883254694  …  -24.990190821574572, -32.21496266548419, -32.25227287587231, -24.95662356094601, -28.688661125109448, -32.228183576737365, -24.971505677608917, -106.33505107570753, -152.63587983919805, -982.6065894965417], 280, 280), sparse([1, 2, 15, 229, 231, 232, 248, 1, 2, 3  …  233, 234, 235, 236, 237, 238, 239, 240, 244, 248], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  248, 248, 248, 248, 248, 248, 248, 248, 248, 248], [-30.831823068641313, -3.897740763777146, -3.8977407637771493, -7.58111657115261, -7.454768484880133, -7.503272132184919, -15.415911534320657, -3.897740763777146, -23.13181874648575, -3.897740763777146  …  -28.883359212372863, -28.74943695888349, -28.88189537206496, -28.765276447741627, -28.868777956495492, -28.78962683280338, -28.81827771127525, -28.84627508598441, -127.50588136824177, -500.9246160098465], 248, 248), sparse([1, 2, 12, 53, 1062, 1064, 1066, 1068, 1, 2  …  419, 424, 592, 597, 598, 1042, 1046, 1056, 1095, 1098], [1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098, 1098], [-19.15897097019969, -3.3322412860713153, -3.1222441990285343, -3.263888888888885, -3.1911116445866736, -3.055555555555552, -2.9856797116100395, -9.371235800641154, -3.3322412860713153, -44.56243033845646  …  -42.11702453218053, -52.25326282429954, -26.27097249339122, -43.172832032195394, -21.47194304168049, -36.984529280208314, -42.99524947260839, -47.27191317004768, -35.26826469011216, -344.4675572255478], 1098, 1098)]], G = Any[[3.1090393294905123 -1.1992888095955179 -2.3443790760237686; 2.4323069210533927 0.5864083237345679 -2.7334364461705283; … ; 0.0 0.0 0.0; 0.0 0.0 0.0], [3.1193845762270946 -1.2032794116203933 -3.4528414016351734; 2.5035289780877705 0.13934379968952898 -3.4910072151694207; … ; 0.0 0.0 -31.582407631794318; 0.0 0.0 0.0], [3.0643379264105572 -1.1820455756555337 -3.489184425021361; 2.4480388841794634 0.27885297738705117 -3.155548435753369; … ; 0.0 0.0 0.0; 0.0 0.0 -23.884028509953588], [3.1090393294905123 -1.1992888095955179 -1.011045742690433; 2.4323069210533927 0.5864083237345679 -3.400103112837189; … ; 5.874081955536858 -11.889833332045157 0.7682874663343249; 3.52827637147994 -20.04284146264382 1.3974288381186803]], volumes = [1918.940233608673, 3787.456852888135, 2864.169516788822, 10699.297136081648])

The Bloch-Torrey PDE takes a magnetic field gradient pulse sequence as an input. Here we consider a ScalarGradient with a PGSE time profile.

dir = [1.0, 0.0, 0.0]
profile = PGSE(2000.0, 6000.0)
b = 1000
g = √(b / int_F²(profile)) / coeffs.γ
gradient = ScalarGradient(dir, profile, g)
ScalarGradient{Float64, PGSE{Float64}}([1.0, 0.0, 0.0], PGSE{Float64}(2000.0, 6000.0), 0.8093302043119762)

SpinDoctor provides a solve function, which has the same base signature for all diffusion MRI problems. The BTPDE is one such problem. They generally take a gradient sequence as an input.

btpde = BTPDE(; model, matrices)
ξ = solve(btpde, gradient)
1824-element Vector{ComplexF64}:
 0.42151898276843097 + 0.4225807888202752im
 0.40258366336102885 + 0.45706295788188617im
 0.45577765300614004 + 0.3827143874574771im
  0.5328395420617785 + 0.21163345238394543im
  0.5464390534024666 - 0.017901451694179997im
 0.49007843969439085 - 0.22617349671501216im
 0.41760247883241997 - 0.3524233969269779im
 0.39627074239026594 - 0.37388243410830074im
  0.4446450105061987 - 0.3014344727022567im
  0.5134150228736285 - 0.13921411958240953im
                     ⋮
 0.13786722927088685 - 0.031723613731392014im
  0.2336829373189278 - 0.19591894192687145im
 0.18640987459140992 - 0.04313466408415677im
 0.15552549398558724 - 0.02588167938464815im
 0.11445698256616671 + 0.010184498969924334im
 0.15689950304035777 + 0.07644876228879041im
  0.2539180429885351 + 0.19837903125473585im
 0.18512946984860332 + 0.06023247406419295im
 0.12048062054335919 - 0.018738097056034738im

Here, ξ is a vector containing the complex-valued magnetization at all degrees of freedom at the echo time TE. We may compute the resulting signal as follows:

compute_signal(matrices.M, ξ)
6268.55354366481 - 0.4912644720153221im

The global mass matrix M is used to compute the integral. We may however be interested in the compartment-wise signals. This requires splitting the magnetization field into the respective compartments. The compartment mass matrices are also available.

ξ_cmpts = split_field(mesh, ξ)
compute_signal.(matrices.M_cmpts, ξ_cmpts)
4-element Vector{ComplexF64}:
  952.1587226283762 + 31.195441404487877im
 1243.7482992899832 - 8.093519658681421im
 1166.3676621178745 - 7.144096023692393im
 2906.2788596285764 - 16.44909019412938im

The final magnetization can be visualized using the plot_field function.

plot_field(mesh, ξ)

In this example, we have computed the complex transverse water proton magnetization field using the finite element method. The measured diffusion MRI signal is the integral of this field, and other quantities of interest, such as the apparent diffusion coefficient (ADC), or the effective diffusion tensor, may easily be obtained from this reference field. Directly solving the BTPDE is thus considered to be the "gold standard" for computing these quantities, as arbitrary precision may be obtained.

However, this is also often the most computationally expensive approach. In the following examples, we will consider some other specialized methods provided by SpinDoctor, each having their own domains of validity, use cases, and computational footprints.

This page was generated using Literate.jl.