Edit on GitHub

Solve BTPDE

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

The built in geometry recipes allow for making various cell configuration. We here consider the case of three twisted axons immersed in an extracellular space (ECS).

setup = CylinderSetup(;
    name = "some-very-real-axons",
    ncell = 3,
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.3,
    height = 40.0,
    bend = 0.0,
    twist = π / 4,
    include_in = false,
    in_ratio = 0.6,
    ecs_shape = :convex_hull,
    ecs_ratio = 0.5,
)
CylinderSetup{Float64}("some-very-real-axons", "meshfiles", 3, 2.0, 6.0, 0.2, 0.3, 40.0, 0.0, 0.7853981633974483, false, 0.6, :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)

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, = create_geometry(setup)
(SpinDoctor.FEMesh{Float64}([[1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  5377, 5402, 5429, 5487, 5522, 5525, 5556, 5560, 5561, 5563], [31, 32, 33, 34, 35, 36, 37, 38, 39, 40  …  5440, 5441, 5443, 5447, 5516, 5523, 5552, 5559, 5562, 5575], [60, 61, 62, 63, 64, 65, 66, 67, 68, 69  …  5456, 5489, 5520, 5548, 5553, 5565, 5568, 5571, 5574, 5579], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  5569, 5570, 5572, 5573, 5576, 5577, 5578, 5580, 5581, 5582]], [[1.2271724573018856 1.464673957512557 … -3.601158865660833 0.4058816632661841; -2.598934435796008 -1.7971428083042387 … -0.3100773841463978 -1.197033205471726; -20.0 -20.0 … -1.0323817355464058 3.9614074838173936], [2.874844919284934 3.108730986825801 … 1.48192204715903 6.876401810925386; -11.217866464179306 -10.417680960996975 … -8.45841952848956 -9.529807312649876; -20.0 -20.0 … -18.274909319517832 14.939764326538537], [5.968809148491234 6.222147523293041 … -0.9049611462682896 -7.839038128627048; 6.290941785212893 7.08174512478809 … 6.193132671067591 8.149341367703233; -20.0 -20.0 … 15.63959852120179 18.49812612990963], [1.2271724573018856 1.464673957512557 … 2.463624949188328 -4.9765350532745165; -2.598934435796008 -1.7971428083042387 … 6.368120355955749 -6.759296306227222; -20.0 -20.0 … 14.27707922725195 -18.64020091779929]], [[320 549 … 360 72; 321 450 … 355 1082; 603 441 … 1081 1079] Matrix{Int64}(undef, 3, 0) … Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0); Matrix{Int64}(undef, 3, 0) [447 434 … 1047 69; 448 549 … 85 1047; 882 547 … 79 1043] … Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0); Matrix{Int64}(undef, 3, 0) Matrix{Int64}(undef, 3, 0) … [18 16 … 829 8; 19 17 … 8 9; 1431 1431 … 1506 1506] Matrix{Int64}(undef, 3, 0); [1614 2758 … 1820 372; 1615 2254 … 1795 4372; 2876 2240 … 4371 4367] [2416 2220 … 4362 325; 2417 2828 … 356 4362; 3974 2826 … 349 4352] … Matrix{Int64}(undef, 3, 0) [18 102 … 3059 4762; 102 17 … 4762 2; 103 101 … 5025 5025]], [[643 117 … 160 1026; 1089 923 … 179 160; 1119 546 … 1084 1084; 1137 922 … 1141 1141], [812 744 … 225 1003; 527 250 … 227 227; 998 742 … 1003 1004; 1051 1065 … 1101 1101], [746 1171 … 55 1308; 748 1173 … 54 53; 994 1357 … 1308 1472; 1340 1456 … 1509 1509], [2088 1562 … 45 24; 3095 2968 … 44 23; 3100 2965 … 4697 4560; 4302 2969 … 5032 5032]]), (points = [2.128207883846267 2.04079828678149 … -8.642374292323556 -8.837956249895907; -1.9339629190351406 -1.1023161557641035 … 10.146286099799731 9.336914288422342; -20.0 -20.0 … 20.0 20.0], facets = [18 102 … 298 247; 102 17 … 297 298; 103 101 … 149 98], facetmarkers = Any[7, 7, 7, 4, 7, 4, 7, 7, 4, 4  …  7, 7, 7, 7, 7, 7, 7, 7, 7, 7], regions = [-1.871792116153733 3.093463592779914 -1.9209403059069228 -7.9364384659005935; -1.9339629190351406 -9.266280858999346 8.093757572126355 8.510091281781285; 0.0 0.0 0.0 0.0]), (radii = [4.0 3.855317646918099 5.02784093379109], centers = [-1.871792116153733 3.093463592779914 -1.9209403059069228; -1.9339629190351406 -9.266280858999346 8.093757572126355]))

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, 30, 638, 640, 642, 1052, 1, 2, 3  …  6757, 6763, 8238, 8287, 8288, 8311, 8448, 8526, 8765, 8783], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783], [0.14411536230946753, 0.017517372544279413, 0.01864903235196686, 0.034850707608734995, 0.035891276258487495, 0.03720697354599878, 0.07205768115473377, 0.017517372544279413, 0.10338396598739119, 0.017134575925347145  …  0.11957015849663083, 0.06683013748219645, 0.17403357420054968, 0.31070704192232207, 0.34244088665562633, 0.30845288522611003, 0.22422859586567614, 0.21244324218152336, 0.19950070446349163, 2.43638189133149], 8783, 8783), S = sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  6757, 6763, 8238, 8287, 8288, 8311, 8448, 8526, 8765, 8783], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783], [0.0037323198717831604, -0.0010291826887019043, -0.0011330474567051072, 0.00024876419612260736, -0.0010373174273221432, -0.00017937604409554706, -0.0006021604510810655, -0.0010291826887019043, 0.0033490520343336613, -0.0009552311946322739  …  -0.004027887294810581, 0.001565970106510206, 0.0009602512294008451, -0.000686171631760155, -0.004229021927761462, -0.004589686456485123, -0.004244516996780964, -0.0007582005242575871, -0.0034799355328381967, 0.024394555578979057], 8783, 8783), R = sparse(Int64[], Int64[], Float64[], 8783, 8783), Mx = SparseArrays.SparseMatrixCSC{Float64, Int64}[sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  6757, 6763, 8238, 8287, 8288, 8311, 8448, 8526, 8765, 8783], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783], [0.13523799910374112, 0.018768303671888593, 0.014163394560854177, 0.03708980165850179, 0.03451960364563238, 0.030193357430884832, 0.04731433558959006, 0.018768303671888593, 0.11870116734309966, 0.020218282081898586  …  -0.519382202873184, -0.31976063435444707, -1.09204243416606, -1.8911483247223062, -1.7719844712650363, -1.6571240484600351, -1.2015559391603396, -1.2905989942973184, -1.017282293526875, -12.084134457437626], 8783, 8783), sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  6757, 6763, 8238, 8287, 8288, 8311, 8448, 8526, 8765, 8783], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783], [-0.34020317438346814, -0.03438936507913983, -0.050986929601869665, -0.07043385190348808, -0.08448465466682933, -0.09922538457373896, -0.1536129740252238, -0.03438936507913983, -0.16553016530978262, -0.021912421899279635  …  -0.8972377760955403, -0.5193357610164643, -1.1900081673440175, -1.9605397716602098, -2.1857268631124462, -1.9093844783984353, -1.6270494788367158, -1.5178861962030448, -1.4776076691889632, -16.03210234599126], 8783, 8783), sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  6757, 6763, 8238, 8287, 8288, 8311, 8448, 8526, 8765, 8783], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783, 8783], [-2.837328513933842, -0.3466979982721967, -0.36909543196601075, -0.6788819125326321, -0.6953934775081954, -0.7247703275270527, -1.4186642569669212, -0.3466979982721967, -2.0390409791154216, -0.3391218151891623  …  -2.2368281745521252, -1.250629222226159, -3.3416777381846074, -5.751194735862685, -6.116168683254481, -5.962048466085882, -4.316242239427791, -3.925309346609161, -3.612695048629814, -45.41484334419949], 8783, 8783)], Q = sparse([1, 2, 30, 638, 640, 642, 3752, 3753, 3781, 6744  …  3642, 3643, 3644, 4168, 4173, 5582, 5601, 8152, 8153, 8154], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  8154, 8154, 8154, 8154, 8154, 8154, 8154, 8154, 8154, 8154], [3.486903462141328e-5, 4.356736284667902e-6, 4.36062873272193e-6, 8.716029087159843e-6, 8.71715229331681e-6, 8.7184882235468e-6, -3.486903462141328e-5, -4.356736284667902e-6, -4.36062873272193e-6, -8.716029087159843e-6  …  -8.656209290291025e-6, -8.650557021946595e-6, -5.192145486085735e-5, 8.650557021946594e-6, 8.653961118191056e-6, 8.656209290291029e-6, 8.65396111819106e-6, 8.656209290291025e-6, 8.650557021946595e-6, 5.192145486085735e-5], 8783, 8783), M_cmpts = Any[sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  800, 801, 913, 1026, 1027, 1084, 1104, 1131, 1138, 1141], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141], [0.14411536230946753, 0.017517372544279413, 0.01864903235196686, 0.034850707608734995, 0.035891276258487495, 0.03720697354599878, 0.07205768115473377, 0.017517372544279413, 0.10338396598739119, 0.017134575925347145  …  0.14748016747567547, 0.07126340980972581, 0.3731518477935909, 0.30338767707625924, 0.09669244003555236, 0.2988093457450153, 0.22539783699052351, 0.2889402905049603, 0.28139490164033537, 2.570693703203142], 1141, 1141), sparse([1, 2, 29, 616, 618, 1062, 1, 2, 3, 614  …  843, 868, 870, 871, 872, 1003, 1004, 1073, 1075, 1101], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2  …  1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101], [0.10233404379340474, 0.01676905595062959, 0.017795948479071688, 0.03337107341763068, 0.03439796594607278, 0.05116702189670237, 0.01676905595062959, 0.09943928199184943, 0.016514356371741413, 0.03295058504529512  …  0.15947878189131334, 0.056848785250546055, 0.11298278849501672, 0.16584372068515776, 0.08428341846159114, 0.2737329180202968, 0.36215901744326134, 0.3969177872424697, 0.32888858454477066, 2.8110679222820476], 1101, 1101), sparse([1, 2, 38, 798, 800, 801, 1477, 1, 2, 3  …  841, 843, 845, 1308, 1310, 1465, 1472, 1480, 1482, 1509], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509], [0.12163929130163595, 0.015691381345758867, 0.01455466666293842, 0.029304032658533596, 0.030573597642120687, 0.03151561299228438, 0.060819645650817974, 0.015691381345758867, 0.199950034462446, 0.01569138134575891  …  0.07006046271328634, 0.1080341764362891, 0.10041080085041197, 0.29313638270209785, 0.2521069055884533, 0.19873237120830342, 0.13622688030723257, 0.3238885041825084, 0.17916597797728068, 1.7624146476026716], 1509, 1509), sparse([1, 2, 30, 2993, 3001, 3009, 5025, 1, 2, 3  …  3006, 3012, 4487, 4536, 4537, 4560, 4697, 4775, 5014, 5032], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032], [0.11601793972821328, 0.014887345931297757, 0.014020262775869503, 0.02983415555699885, 0.029101361156939384, 0.028174814307107793, 0.05800896986410664, 0.014887345931297757, 0.15503655397533617, 0.012230884275906325  …  0.11957015849663083, 0.06683013748219645, 0.17403357420054968, 0.31070704192232207, 0.34244088665562633, 0.30845288522611003, 0.22422859586567614, 0.21244324218152336, 0.19950070446349163, 2.43638189133149], 5032, 5032)], S_cmpts = Any[sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  800, 801, 913, 1026, 1027, 1084, 1104, 1131, 1138, 1141], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141], [0.0037323198717831604, -0.0010291826887019043, -0.0011330474567051072, 0.00024876419612260736, -0.0010373174273221432, -0.00017937604409554706, -0.0006021604510810655, -0.0010291826887019043, 0.0033490520343336613, -0.0009552311946322739  …  -0.004262901908733892, 0.0025274256772752834, -0.008830597834310477, -0.0013024809241955883, 0.0026465166670347376, -0.003911923763319232, 0.0007284076582203906, -0.0044690239666460504, -0.005861454581139328, 0.027624520838851144], 1141, 1141), sparse([1, 2, 29, 616, 618, 1062, 1, 2, 3, 614  …  843, 868, 870, 871, 872, 1003, 1004, 1073, 1075, 1101], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2  …  1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101], [0.0039161614009296615, -0.0010424594041866977, -0.0012891737911371567, 0.00038898949087715236, -0.0014236217213977633, -0.0005498959750851963, -0.0010424594041866977, 0.003557932920766926, -0.0009282002409012727, 0.00018580249357365256  …  -0.0005545721413416811, 0.0014165467928998906, -0.0003824604185940485, 0.00037603274531279216, 1.0391529205492093e-5, -0.0028322250610052285, -0.007513591785730365, -0.004815437651782874, -0.0038603958939315095, 0.02715435659275021], 1101, 1101), sparse([1, 2, 38, 798, 800, 801, 1477, 1, 2, 3  …  841, 843, 845, 1308, 1310, 1465, 1472, 1480, 1482, 1509], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509], [0.003415557936263168, -0.0013281467163810496, -0.0008493687430861293, 0.00012042287544662643, -0.0007781616779873435, 5.880989118460277e-5, -0.0006391135654398747, -0.0013281467163810496, 0.00479959635947959, -0.0011961682101291474  …  -0.0009305682402290665, -0.0021008500131139917, -0.00018500558693026055, -0.010869481465412592, -0.0025094655149133396, -0.0026948489492296325, -0.0008208434542327247, -0.005303224269610718, -0.0029295688361452568, 0.027798718840329383], 1509, 1509), sparse([1, 2, 30, 2993, 3001, 3009, 5025, 1, 2, 3  …  3006, 3012, 4487, 4536, 4537, 4560, 4697, 4775, 5014, 5032], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032], [0.0035418147334346406, -0.0008320817489714874, -0.0010993738988919614, 0.00010541659740629318, -0.0006300586850641889, -0.00017332338529440778, -0.000912393612618888, -0.0008320817489714874, 0.0043030918415647105, -0.001092380773628524  …  -0.004027887294810581, 0.001565970106510206, 0.0009602512294008451, -0.000686171631760155, -0.004229021927761462, -0.004589686456485123, -0.004244516996780964, -0.0007582005242575871, -0.0034799355328381967, 0.024394555578979057], 5032, 5032)], Mx_cmpts = Vector{Any}[[sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  800, 801, 913, 1026, 1027, 1084, 1104, 1131, 1138, 1141], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141], [0.13523799910374112, 0.018768303671888593, 0.014163394560854177, 0.03708980165850179, 0.03451960364563238, 0.030193357430884832, 0.04731433558959006, 0.018768303671888593, 0.11870116734309966, 0.020218282081898586  …  0.21360851861259122, 0.08031077408663592, -0.10905917605344846, -0.10747281872570905, -0.020231774450526367, 0.13656605721454185, 0.08512430482415308, 0.1423154604237319, 0.11966474454042167, 1.213598154704466], 1141, 1141), sparse([1, 2, 29, 616, 618, 1062, 1, 2, 3, 614  …  843, 868, 870, 871, 872, 1003, 1004, 1073, 1075, 1101], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2  …  1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101], [0.2737323216911787, 0.04629710107099697, 0.04473418565970973, 0.09366167764813035, 0.09279585878493765, 0.12287856728390698, 0.04629710107099697, 0.28308454607116496, 0.04729089202707013, 0.09618774870133637  …  1.2210207346825273, 0.40359439030915656, 0.8435246672306673, 1.317247083446249, 0.6579132669504534, 1.9911276893696628, 2.6435563818757872, 2.6616187248074756, 2.127062441291757, 19.218660011581445], 1101, 1101), sparse([1, 2, 38, 798, 800, 801, 1477, 1, 2, 3  …  841, 843, 845, 1308, 1310, 1465, 1472, 1480, 1482, 1509], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509], [0.6806545193213576, 0.09015189370757387, 0.07910326333000278, 0.15914097848936015, 0.17001669537490227, 0.1790639860946528, 0.3208113642776256, 0.09015189370757387, 1.1432069869692463, 0.0939239096168234  …  -0.6266525657585522, -0.9319065252202214, -0.8158778278198413, -2.2488734239386687, -1.7494713199895724, -1.5348262851033485, -1.07342750547924, -2.3154257683958814, -1.4287831864255547, -13.821110489279308], 1509, 1509), sparse([1, 2, 30, 2993, 3001, 3009, 5025, 1, 2, 3  …  3006, 3012, 4487, 4536, 4537, 4560, 4697, 4775, 5014, 5032], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032], [0.1754460526402889, 0.02442227407164004, 0.018626365714435336, 0.04872525398417343, 0.04463868751406572, 0.03890480824804698, 0.10438770565452954, 0.02442227407164004, 0.2857119772492599, 0.02146181208145352  …  -0.519382202873184, -0.31976063435444707, -1.09204243416606, -1.8911483247223062, -1.7719844712650363, -1.6571240484600351, -1.2015559391603396, -1.2905989942973184, -1.017282293526875, -12.084134457437626], 5032, 5032)], [sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  800, 801, 913, 1026, 1027, 1084, 1104, 1131, 1138, 1141], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141], [-0.34020317438346814, -0.03438936507913983, -0.050986929601869665, -0.07043385190348808, -0.08448465466682933, -0.09922538457373896, -0.1536129740252238, -0.03438936507913983, -0.16553016530978262, -0.021912421899279635  …  -0.23958242207818944, -0.1533722838336731, -0.4209593124919463, -0.1826284648622132, -0.2258793871870895, -0.36775353402054417, -0.23927292775562925, -0.5522381927277867, -0.6053172754632452, -2.8974050462702454], 1141, 1141), sparse([1, 2, 29, 616, 618, 1062, 1, 2, 3, 614  …  843, 868, 870, 871, 872, 1003, 1004, 1073, 1075, 1101], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2  …  1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101], [-1.1155223409738944, -0.17762915644424657, -0.19998049462225395, -0.3550579669688583, -0.377472647336507, -0.5469195975689773, -0.17762915644424657, -1.017116419998232, -0.1636552776752972, -0.32815603662437026  …  -1.5693851671818895, -0.6013890764579226, -1.1820056352124075, -1.593172028404115, -0.8522573277237312, -2.470245295132048, -3.1768421097748276, -3.9093126258380586, -3.2730451393332016, -26.955561283726553], 1101, 1101), sparse([1, 2, 38, 798, 800, 801, 1477, 1, 2, 3  …  841, 843, 845, 1308, 1310, 1465, 1472, 1480, 1482, 1509], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509], [0.7699570609888606, 0.1051445941067693, 0.08636318476008575, 0.17671262846465174, 0.19448386927525857, 0.21026474337675738, 0.3843322518076261, 0.1051445941067693, 1.4592983550649934, 0.12247985717799244  …  0.5806915631812638, 0.9487358826041197, 0.9149428947449872, 2.2783518379257295, 1.9499791180209338, 1.6769138315440937, 1.1818835716477374, 2.7688371680604873, 1.3402952447155227, 14.407887323861859], 1509, 1509), sparse([1, 2, 30, 2993, 3001, 3009, 5025, 1, 2, 3  …  3006, 3012, 4487, 4536, 4537, 4560, 4697, 4775, 5014, 5032], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032], [-0.30288828599216666, -0.03312554380642633, -0.04200398253945114, -0.06812350802969447, -0.0758932890991762, -0.08250463666115229, -0.15336410248355342, -0.03312554380642633, -0.25912597427190215, -0.01390852412501704  …  -0.8972377760955403, -0.5193357610164643, -1.1900081673440175, -1.9605397716602098, -2.1857268631124462, -1.9093844783984353, -1.6270494788367158, -1.5178861962030448, -1.4776076691889632, -16.03210234599126], 5032, 5032)], [sparse([1, 2, 30, 638, 640, 642, 1052, 1, 2, 3  …  800, 801, 913, 1026, 1027, 1084, 1104, 1131, 1138, 1141], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141, 1141], [-2.837328513933842, -0.3466979982721967, -0.36909543196601075, -0.6788819125326321, -0.6953934775081954, -0.7247703275270527, -1.4186642569669212, -0.3466979982721967, -2.0390409791154216, -0.3391218151891623  …  0.5657379381384585, 0.2706414081513551, 1.3247405455864443, 1.3206667712079039, 0.3404582184518021, 1.4450467064804953, 0.6042018839542208, 0.9386441511389528, 1.2000772811543627, 10.004956143044973], 1141, 1141), sparse([1, 2, 29, 616, 618, 1062, 1, 2, 3, 614  …  843, 868, 870, 871, 872, 1003, 1004, 1073, 1075, 1101], [1, 1, 1, 1, 1, 1, 2, 2, 2, 2  …  1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101, 1101], [-2.0184437761332186, -0.3318875656895439, -0.3522114803149605, -0.6500581007896422, -0.6701680794716334, -1.009221888066609, -0.3318875656895439, -1.9612206941413723, -0.3268466365240488, -0.6418580761633724  …  2.2687082726347647, 0.8962316223028352, 1.771333690115169, 2.6038449262250487, 1.331387144958804, 3.9257942556593166, 5.559237443396291, 6.26251258814385, 4.636805092241713, 42.08889495136903], 1101, 1101), sparse([1, 2, 38, 798, 800, 801, 1477, 1, 2, 3  …  841, 843, 845, 1308, 1310, 1465, 1472, 1480, 1482, 1509], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509, 1509], [-2.394705307993995, -0.3105585891348109, -0.2880611110373229, -0.5707978549805339, -0.5923634543160883, -0.6138840395058762, -1.1973526539969974, -0.3105585891348109, -3.9509390970543823, -0.3105585891348117  …  1.3135548272811126, 2.01421483402462, 1.8713598281689723, 5.616135287551398, 4.562194163919546, 3.5131155774021625, 2.619426326157715, 5.974751244777006, 3.267196969744559, 32.66255208482103], 1509, 1509), sparse([1, 2, 30, 2993, 3001, 3009, 5025, 1, 2, 3  …  3006, 3012, 4487, 4536, 4537, 4560, 4697, 4775, 5014, 5032], [1, 1, 1, 1, 1, 1, 1, 2, 2, 2  …  5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032, 5032], [-2.2840628233054963, -0.2946453882236014, -0.2774843674390839, -0.5811382943192065, -0.5638388724157006, -0.5488079152785197, -1.1420314116527481, -0.2946453882236014, -3.063200305255279, -0.24206958462731268  …  -2.2368281745521252, -1.250629222226159, -3.3416777381846074, -5.751194735862685, -6.116168683254481, -5.962048466085882, -4.316242239427791, -3.925309346609161, -3.612695048629814, -45.41484334419949], 5032, 5032)]], G = Any[[0.6434424434047774 -0.25731790838979657 -0.6049557970594269; 0.5146771833386863 -0.07255047635654992 -0.5674251721865191; … ; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.4868893493279125 -0.17677048516684346 -0.6466451440753818; 0.5135065759225813 -0.06797087169770381 -0.6300091932789103; … ; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.3326437802088955 -0.09769988098075813 -0.4254341229621356; 0.507692882789164 -0.09845392378710108 -1.2676046701777688; … ; 0.0 0.0 0.0; 0.0 0.0 0.0], [0.6434424434047774 -0.25731790838979657 -0.4888150580341625; 0.5146771833386863 -0.07255047635654992 -1.0523971290224126; … ; 0.0 0.0 0.0; 0.0 0.0 0.0]], volumes = [1995.7625763798055, 1853.0307929680541, 3161.8999333686215, 8537.013937066502])

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}([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)
8783-element Vector{ComplexF64}:
  0.3768999655836339 + 0.3309834594154104im
 0.36358192810329737 + 0.34229683652152737im
 0.35747370221571445 + 0.3471241695722385im
  0.3581520546535996 + 0.34680852022902586im
  0.3789057313288629 + 0.3183342804420766im
  0.4030183490755389 + 0.2898576849647053im
  0.4317375726663847 + 0.2575902341926537im
  0.4668676633333275 + 0.19886293586111095im
 0.49744400075876455 + 0.13282271887008268im
  0.5213405379030064 + 0.05802487080635975im
                     ⋮
  0.5695945534139173 + 0.10134319755603721im
  0.5264862627507649 - 0.06265173251602779im
   0.579906372631599 + 0.24040165988132478im
 0.45086396654981203 - 0.15059733687532043im
  0.4906752100691949 - 0.24109251299010173im
  0.4341579017701663 - 0.10710707832773204im
  0.5870462467041845 - 0.18030398569286865im
 0.49155487369159384 + 0.0706614283827105im
  0.5314925105455783 - 0.1562988278336503im

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, ξ)
7060.181187493257 + 0.08015554231403499im

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}:
  962.0588969196565 + 1.5498878994450136im
  951.8171146808451 - 4.151787654661453im
 1271.9383499938085 - 3.149635661743872im
 3874.3668258989474 + 5.831690959274354im

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.