API reference

SpinDoctor Module

SpinDoctor

Diffusion MRI simulation toolbox: finite element simulations of the Bloch-Torrey equation and related models in geometries of cells and extracellular space.

The main functions and types are exported. Lower level building blocks (meshes, assembly, geometry files, gradient sequence internals, callback hooks) are public but not exported, and are used as SpinDoctor.name.

source

The main functions and types are exported. The building blocks listed under Public, not exported are used with the module name, e.g. SpinDoctor.FEMesh.

Setup recipes and geometry

SpinDoctor comes with recipes for creating finite element meshes of plates, slabs, disks, cylinders, spheres, polyhedra and neurons, possibly immersed in an extracellular space (ECS).

SpinDoctor.PlateSetup Type

PlateSetup(; depth = 1.0, widths = [1.0, 1.0, 1.0, 1.0, 1.0], refinement = Inf)

Setup recipe for stacked plates in 2D, with the given widths along the x-axis and extent depth along the y-axis. The triangles have an area of at most refinement.

source

SpinDoctor.SlabSetup Type

SlabSetup(;
    widths = [1.0, 1.0, 1.0, 1.0, 1.0],
    depth = 5.0,
    height = 1.0,
    bend = 0.0,
    twist = 0.0,
    refinement = Inf,
)

Setup recipe for stacked slabs: the extrusion of a PlateSetup.

source

SpinDoctor.DiskSetup Type

DiskSetup(;
    ncell = 3,
    layersizes = [1.0],
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.6,
    nsidewall = 12,
    ecs = NoECS(),
    refinement = Inf,
)

Setup recipe for ncell disks in 2D, possibly immersed in an extracellular space ecs.

• The radii are between rmin and rmax, and the gaps between neighboring cells between dmin and dmax times the mean radius (see create_cells).

• Each cell has layers with the relative radii layersizes (increasing, the last one being one), e.g. [0.6, 1.0] for an axon with a myelin sheath.

• A cell of mean radius has nsidewall points on its boundary.

• The triangles have an area of at most refinement.

source

SpinDoctor.CylinderSetup Type

CylinderSetup(;
    ncell = 1,
    layersizes = [1.0],
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.6,
    nsidewall = 12,
    ecs = NoECS(),
    height = 1,
    bend = 0.0,
    twist = 0.0,
    refinement = Inf,
)

Setup recipe for parallel cylinders: the extrusion of a DiskSetup.

source

SpinDoctor.SphereSetup Type

SphereSetup(;
    ncell = 3,
    layersizes = [1.0],
    rmin = 2.0,
    rmax = 6.0,
    dmin = 0.2,
    dmax = 0.6,
    nsidewall = 200,
    ecs = NoECS(),
    refinement = Inf,
)

Setup recipe for ncell spheres, possibly immersed in an extracellular space ecs. The parameters are those of DiskSetup, with nsidewall points on the surface of a cell of mean radius, and a maximum tetrahedron volume refinement.

source

SpinDoctor.ExtrusionSetup Type

ExtrusionSetup(; groundsetup, height = 1.0, bend = 0.0, twist = 0.0, refinement = Inf)

Setup recipe for the extrusion of the 2D setup groundsetup along the z-axis over height, possibly bent and twisted (see deform_domain!). The tetrahedra have a volume of at most refinement.

source

SpinDoctor.NeuronSetup Type

NeuronSetup(; meshfile, ecs = NoECS(), refinement = Inf)

Setup recipe for a neuron, possibly wrapped in a ECS. The neuron is given by the volume mesh meshfile in the Gmsh format, e.g. from the RealNeuronMeshes repository.

source

SpinDoctor.VoronoiSetup Type

VoronoiSetup(;
    ncell = 20,
    boxsize = 20.0,
    fraction = 0.8,
    lloyd = 10,
    meshsize = boxsize / ncell^(1 / 3) / 8,
    refinement = Inf,
)

Setup recipe for ncell polyhedral cells in a periodic cubic box of side boxsize: the Voronoi cells of random seed points (made more regular by lloyd Lloyd iterations), separated by a uniform gap such that the cells fill the volume fraction fraction. The cells are cut by the faces of the box, their periodic images enter on the opposite faces, and the extracellular space fills the gaps (as for a PeriodicBoxECS).

Unlike packings of spheres, which cannot fill more than about 64 % of the volume, polyhedral cells can fill any fraction below one. The surfaces are triangulated with triangles of size meshsize, such that the opposite faces of the box have the same triangles, and the tetrahedra have a volume of at most refinement. The mesh is periodic, so that periodic boundary conditions may be used with assemble_matrices(model; periodic = true). The coefficients are given as for a SphereSetup with one layer.

source

SpinDoctor.NoECS Type

NoECS()

No extracellular space: the cells are the only compartments.

source

SpinDoctor.BoxECS Type

BoxECS(; margin = 1.0)

Extracellular space filling the axis-aligned bounding box of the cells, extended by margin on each side. The cells do not touch the faces of the box, so the ECS contains a free layer of thickness margin along the faces.

source

SpinDoctor.ConvexHullECS Type

ConvexHullECS(; margin = 1.0, threshold = pi / 20, maxiter = 10)

Extracellular space filling the convex hull of the cells, moved outwards by margin. Facets of the hull with angles below threshold are refined for up to maxiter iterations (3D only).

source

SpinDoctor.PeriodicBoxECS Type

PeriodicBoxECS(; fraction = 0.5)

Extracellular space of a sample of a periodic tissue. The box of the cell configuration (cells.box, e.g. from periodic_packing) is filled with the cells, cut by its faces, and their periodic images, and the ECS fills the rest. Unlike BoxECS, there is no free layer of ECS along the faces, and the geometry is periodic. The mesh is periodic too (the opposite faces of the box have the same points), so that periodic boundary conditions may be used with assemble_matrices(model; periodic = true); otherwise, the box faces are impermeable boundaries.

If no cell configuration is given to create_geometry, a packing with the cell volume fraction fraction is created (see create_cells). Available for disks and cylinders. Fractions up to about 0.8 (counting the gaps as part of the cells) require many cells, e.g. twenty or more. As in the other recipes, the disks are polygons inscribed in the circles (with nsidewall sides for the mean radius), so that their area is slightly smaller than that of the circles (by 1.1 % for 24 sides).

source

SpinDoctor.periodic_packing Function

periodic_packing(
    radii;
    dim = 2,
    fraction,
    gap = 0.0,
    rng = Random.default_rng(),
    maxiter = 200_000,
)

Pack disks (dim = 2) or spheres (dim = 3) with the given radii in a periodic box, without overlaps and with gaps of at least gap between neighbors. The box size is chosen so that the cells fill the volume fraction fraction of the box.

The positions are found by a collective rearrangement: starting from random positions, the radii grow towards their final values while the overlaps are removed by moving overlapping pairs apart. Random packings of polydisperse disks jam at fractions of about 0.84 and spheres at about 0.64 (counting the gaps as part of the cells); higher fractions cannot be reached. If the packing jams (or takes more than maxiter relaxation sweeps) within one percent of the final radii, the radii are reduced accordingly (with a warning); otherwise an error is thrown.

Returns a cell configuration (; centers, radii, box), with centers of size (dim, ncell) inside box = (lower, upper), for create_geometry with a PeriodicBoxECS.

source

SpinDoctor.create_geometry Function

mesh, surfaces, cells = create_geometry(
    setup;
    cells = nothing,
    rng = Random.default_rng(),
    cache = nothing,
    savedir = nothing,
    recreate = false,
)

Create the cells, the surface triangulation and the finite element mesh of setup.

Random cell configurations are drawn from rng. For disks, cylinders and spheres, a cell configuration (; centers, radii) (as returned by create_cells) may be given instead, with centers of size (dim, ncell) and radii of length ncell, e.g. from a packing algorithm.

By default, no files are written. The geometry may be stored in one of two ways:

• cache: a JLD2 file (requires using JLD2, see save_geometry). If the file exists and was created for the same setup (and cells), the geometry is loaded from it. Otherwise, it is created and saved.

• savedir: a directory with text files in the TetGen format, compatible with the MATLAB version of SpinDoctor (cells, mesh.node, mesh.poly, mesh.1.node, mesh.1.face, mesh.1.ele), and a description of the setup (setup). Existing files are reused if they were created for the same setup.

With recreate = true, the geometry is always created, and the files are overwritten.

source

SpinDoctor.save_geometry Function

save_geometry(file, setup, mesh, surfaces, cells)

Save a geometry created by create_geometry to the JLD2 file file, together with a description of setup. Only plain arrays are stored, so the file does not depend on the SpinDoctor version. Requires using JLD2.

source

SpinDoctor.load_geometry Function

(; setup, mesh, surfaces, cells) = load_geometry(file)

Load a geometry saved with save_geometry. The setup is returned as the description string that was saved. Requires using JLD2.

source

SpinDoctor.compartment_mesh Function

compartment_mesh(mesh, icmpt)

Extract compartment icmpt of mesh as a mesh of its own, with one compartment. The boundary numbering is kept (facets has one row and the same columns as mesh.facets), and point_map still maps to the points of the global mesh.

A model of the compartment alone can then be built with coefficients for one compartment, for example to solve a problem in the extracellular space only.

source

SpinDoctor.box_faces Function

box_faces(mesh, icmpt = length(mesh.points); atol = 1e-8 * diameter)

Find the boundary facets of compartment icmpt (by default the last one, which is the ECS in the setups with an ECS) lying on the faces of its bounding box. This applies to axis-aligned box-shaped compartments, such as a BoxECS that has not been bent or twisted.

Returns a named tuple with the facets (in the local numbering of the compartment) on each face: (; xmin, xmax, ymin, ymax) in 2D and (; xmin, xmax, ymin, ymax, zmin, zmax) in 3D. The points on a face are given by unique(facets).

source

SpinDoctor.compartment_volumes Function

compartment_volumes(mesh)

Volumes (areas in 2D) of the compartments of mesh.

source

Model and matrices

SpinDoctor.coefficients Function

coefficients(setup; D, T2, rho, kappa, gamma)

Distribute the coefficients over the compartments and boundaries of setup, and return them as a named tuple (; D, T2, rho, kappa, gamma), from which a Model may be created. The coefficients are given as follows:

• plates and slabs: D, T2 and rho are vectors with one value per plate, and kappa = (; interfaces, boundaries) contains the permeabilities between consecutive plates and the relaxivities of the outer boundaries of the plates;

• disks, cylinders and spheres: D, T2 and rho are named tuples (; cell, ecs), with one value per layer for cell, and kappa = (; cell_interfaces, cell_boundaries, cell_ecs, ecs) contains the permeabilities between consecutive layers, the relaxivities of the outer boundaries of the layers (e.g. the ends of extruded cylinders), the permeability between the cells and the ECS and the relaxivity of the outer boundary of the ECS;

• neurons: D, T2 and rho are named tuples (; neuron, ecs), and kappa = (; neuron_ecs, neuron, ecs).

The diffusion tensors D are matrices (µm²/µs), T2 times in µs, permeabilities and relaxivities in µm/µs and gamma the gyromagnetic ratio (see Model).

source

SpinDoctor.analytical_coefficients Function

analytical_coefficients(setup, coeffs)

Coefficients of the analytical matrix formalism (in SI units) for plates, slabs, disks, cylinders and spheres with one cell, from the coefficients coeffs of coefficients. Disks, cylinders and spheres have the mean radius.

source

SpinDoctor.Model Type

Model(; mesh, rho, D, T2, kappa, gamma)

Finite element discretized biological model:

• mesh: finite element mesh split into compartments (FEMesh);

• rho: initial spin density of each compartment;

• D: diffusion tensor of each compartment (µm²/µs): a matrix, or a vector with one matrix per element of the compartment (e.g. from diffusion tensor imaging);

• T2: T2 relaxation time of each compartment (µs);

• kappa: permeability of each boundary (µm/µs), for interfaces between compartments and for the outer boundaries (as a boundary relaxivity);

• gamma: gyromagnetic ratio (rad / (µs µm T/m), 2.67513e-4 for water protons).

The coefficients for a setup recipe may be created with coefficients.

source

SpinDoctor.initial_conditions Function

initial_conditions(model)

Initial magnetization (the spin density) at all the points of the mesh.

source

SpinDoctor.assemble_matrices Function

assemble_matrices(model; periodic = false)

Assemble the finite element matrices of model.

Returns a named tuple with

• the global (block diagonal) mass matrix M, stiffness matrix S (with the diffusion tensors), relaxation matrix R (M divided by T2), first order moment matrices Mx (, one per axis d) and flux matrix Q (permeable interfaces and outer boundaries, see assemble_flux_matrix);

• the compartment blocks M_cmpts, S_cmpts and Mx_cmpts;

• the surface integrals G, one matrix per compartment with the integrals   over the boundary of the compartment with outward normal n (one row per point, one column per axis);

• the compartment volumes.

With periodic = true (or one Bool per axis), the mesh is a sample of periodic tissue in its bounding box, such as the meshes of PeriodicBoxECS and VoronoiSetup, and each compartment is periodic along the given axes: the points on opposite faces of the box are identified. The solvers then use periodic boundary conditions instead of impermeable box faces. The matrices are still those of the whole mesh, and the named tuple also contains

• the prolongation P from the periodic points to all the points (see periodic_restriction), its compartment blocks P_cmpts, and the index of one point for each periodic point (masters);

• the matrices Cx (, one per axis d, see assemble_derivative_matrix).

Without periodicity, these fields are nothing.

source

Gradients

SpinDoctor.PGSE Type

PGSE(delta, Delta)

Pulsed gradient spin echo (PGSE) time profile with pulse duration delta and time Delta between the starts of the two pulses.

source

SpinDoctor.DoublePGSE Type

DoublePGSE(delta, Delta, tpause)

Double pulsed gradient spin echo time profile with four pulses of duration delta, separated by pauses of duration Delta - delta, tpause, and Delta - delta respectively.

source

SpinDoctor.CosOGSE Type

CosOGSE(delta, Delta, nperiod)

Oscillating gradient spin echo time profile with two cosine pulses of duration delta, starting at times 0 and Delta, with nperiod periods per pulse.

source

SpinDoctor.SinOGSE Type

SinOGSE(delta, Delta, nperiod)

Oscillating gradient spin echo time profile with two sine pulses of duration delta, starting at times 0 and Delta, with nperiod periods per pulse.

source

SpinDoctor.ScalarGradient Type

ScalarGradient(dir, profile, amplitude)
ScalarGradient(dir, profile; bvalue, gamma)

Gradient sequence g(t) = amplitude * profile(t) * dir, with a constant direction dir (normalized upon construction), a time profile profile and an amplitude amplitude.

Instead of the amplitude, the b-value may be given together with the gyromagnetic ratio gamma.

source

SpinDoctor.GeneralGradient Type

GeneralGradient(; gvec, TE)

General gradient sequence g(t) = gvec(t), whose direction and amplitude may vary in time, with echo time TE.

source

SpinDoctor.bvalue Function

bvalue(gradient, gamma)

b-value of the scalar gradient sequence gradient for the gyromagnetic ratio gamma.

source

SpinDoctor.unitcircle Function

unitcircle(ndirection; half = false, normal = [0, 0, 1])

Create ndirection directions (columns) uniformly distributed on the unit circle in the plane orthogonal to normal.

If half is true, the directions lie on a half-circle instead of the whole circle.

source

SpinDoctor.unitsphere Function

unitsphere(ndirection; half = false, normal = [0, 0, 1])

Create ndirection directions (columns) uniformly distributed on the unit sphere (Fibonacci lattice).

If half is true, the directions are distributed on the hemisphere defined by normal instead.

source

Solvers

SpinDoctor allows for considering a wide range of diffusion MRI problems, each with its own solver function.

SpinDoctor.solve_btpde Function

solve_btpde(model, matrices, gradient, solver = ThetaSolver(); callbacks = [])

Solve the Bloch-Torrey partial differential equation using P1 finite elements in space and the time stepping solver, and return the magnetization at the echo time.

The time stepping is a ThetaSolver. With OrdinaryDiffEq, its ODE algorithms may also be used, e.g. using OrdinaryDiffEqBDF; solve_btpde(model, matrices, gradient, QNDF(); abstol, reltol).

The callbacks (see AbstractCallback) are called after each time step.

With periodic matrices (assemble_matrices(model; periodic = true), only with ThetaSolver), the compartments are periodic in the box of the mesh: the magnetization is  , where m is periodic and   is the wave vector (see wavevector). The time steps then use one factorization per interval only where the gradient is zero.

source

SpinDoctor.ThetaSolver Type

ThetaSolver(; theta = 0.5, timestep = nothing, arraytype = Array)

Theta-scheme time stepping for the BTPDE (theta = 0.5: Crank-Nicolson, second order; theta = 1: implicit Euler, first order), with steps of at most timestep in each characteristic interval of the gradient (by default a thousandth of the echo time). The gradient is evaluated at the middle of each step.

For interval-wise constant scalar gradients (e.g. PGSE, DoublePGSE), the linear systems are solved with one sparse factorization per interval. For other gradients, they are solved with GMRES, preconditioned with a factorization of the diffusion part (one per step size), so that the matrices are not refactorized at every step.

With arraytype = CuArray (after using CUDA, CUDSS), the time stepping runs on an NVIDIA GPU, with the sparse direct solvers of cuDSS.

source

SpinDoctor.solve_hadc Function

solve_hadc(
    model,
    matrices,
    gradient;
    timestep = echotime(gradient) / 1000,
    arraytype = Array,
)

Compute the ADC of each compartment in the direction of the scalar gradient gradient, using the homogenized ADC model (HADC). See solve_hadc_tensors.

source

SpinDoctor.solve_hadc_tensors Function

solve_hadc_tensors(
    model,
    matrices,
    profile;
    timestep = echotime(profile) / 1000,
    arraytype = Array,
)

Compute the apparent diffusion tensor of each compartment for the time profile profile, using the homogenized ADC model (HADC) [14]. The ADC in the direction u is u' * D * u.

The HADC model assumes impermeable compartment boundaries: permeabilities, initial spin densities and T2 relaxation are not taken into account. The ADC of a sample may be estimated by weighting the compartment ADCs with compartment_weights.

The auxiliary diffusion problem is solved with the Crank-Nicolson scheme, with steps of at most timestep in each characteristic interval of the profile. By linearity, the ADC in any direction follows from one solve per axis. With arraytype = CuArray (after using CUDA, CUDSS), it is solved on an NVIDIA GPU.

With periodic matrices (assemble_matrices(model; periodic = true)), the compartments are periodic, and the faces of the box are not boundaries. With one diffusion tensor per element, the tensors are relative to the volume average of the diffusion tensor.

source

SpinDoctor.compute_adc_sta Function

compute_adc_sta(model, gradient)

Compute the ADC of each compartment in the short diffusion time approximation (STA), for the direction of the scalar gradient gradient.

source

SpinDoctor.solve_karger Function

solve_karger(model, difftensors, gradient; timestep)

Solve the finite pulse Karger model (FPK) using precomputed effective diffusion tensors difftensors (one per compartment), and return the compartment signals.

The compartments exchange magnetization through their interfaces with rates given by the permeabilities, interface areas and volumes. The compartment signals are integrated with a fourth order Magnus method, with steps of at most timestep within each characteristic interval of the gradient sequence.

source

SpinDoctor.solve_laplace Function

solve_laplace(model, matrices; neig_max)

Compute the neig_max smallest Laplace eigenvalues values, the eigenfunctions funcs, the first order moments of products of pairs of eigenfunctions moments and the relaxation matrix massrelax in the eigenfunction basis.

See [21], [4].

source

SpinDoctor.solve_mf Function

solve_mf(model, matrices, lap_eig, gradient; ninterval = 500)

Compute the magnetization at the echo time using the matrix formalism, in the basis of the Laplace eigenfunctions lap_eig (see solve_laplace). Gradients that are not interval-wise constant are approximated by ninterval constant intervals.

See [21], [4].

source

SpinDoctor.eig2length Function

eig2length(lambda, D)

Convert a Laplace eigenvalue to a length scale, for the diffusivity D.

source

SpinDoctor.length2eig Function

length2eig(length, D)

Convert a length scale to a Laplace eigenvalue, for the diffusivity D.

source

SpinDoctor.limit_lengthscale Function

limit_lengthscale(lap_eig, lambda_max)

Keep the Laplace eigenmodes with eigenvalues up to lambda_max (length scales of at least eig2length(lambda_max, D)).

source

SpinDoctor.compute_mf_diffusion_tensor Function

compute_mf_diffusion_tensor(mesh, M, lap_eig, gradient::ScalarGradient)

Compute the effective diffusion tensor of the sample using the matrix formalism, from the mass matrix M and the Laplace eigendecomposition lap_eig.

source

SpinDoctor.solve_analytical_laplace Function

solve_analytical_laplace(; rho, r, D, kappa, T2, gamma, dim, eiglim, eigstep)

Solve the analytical radial Laplace eigenvalue problem of multilayered plates (dim = 1), cylinders (dim = 2) or spheres (dim = 3), with the coefficients of analytical_coefficients, for eigenvalues up to eiglim, searched with steps eigstep. Return the coefficients together with the eigendecomposition.

source

SpinDoctor.solve_analytical_mf Function

solve_analytical_mf(laplace, volumes, gradient)

Compute the signal in multilayered plates, cylinders or spheres using an analytical matrix formalism solution, from the eigendecomposition laplace of solve_analytical_laplace and the compartment volumes. Only implemented for PGSE sequences.

See [22] and [6].

source

Post-processing and output

SpinDoctor.compute_signal Function

compute_signal(M, u)

Compute the signal from the magnetization u, using the mass matrix M for integration.

Given a mesh mesh and the compartment mass matrices M_cmpts, the compartment signals are compute_signal.(M_cmpts, split_field(mesh, u)).

source

SpinDoctor.split_field Function

split_field(mesh, u)

Split the field u (with values at all the points of mesh) into a vector with a view of the values of each compartment.

source

SpinDoctor.compartment_weights Function

compartment_weights(model; TE = 0)

Relative contributions of the compartments to the signal at b = 0: the initial spin density times the volume, attenuated by T2 relaxation at the echo time TE, and normalized to sum to one.

The ADC of the whole sample in the low-b limit is the weighted average of the compartment ADCs without exchange, e.g. compartment_weights(model; TE)' * adcs for compartment ADCs from solve_hadc.

source

SpinDoctor.fit_adc Function

fit_adc(bvalues, signals; degree = length(bvalues) - 1)

Compute the apparent diffusion coefficient (ADC) as the initial slope of the logarithm of the signals with respect to the bvalues, by fitting a polynomial of degree degree (by default interpolating all the points; use degree = 1 for a log-linear least squares fit).

source

SpinDoctor.fit_tensors Function

fit_tensors(directions, adcs)

Fit effective diffusion tensors to the ADCs adcs[i] (one value per compartment) in the directions[:, i]. The components of the symmetric tensors are fitted by least squares.

source

SpinDoctor.save_field Function

save_field(filename, mesh, u; name = "Magnetization")

Save the field u (real and imaginary parts) to the VTK multiblock file filename.vtm, with one block per compartment. It may then be visualized in ParaView.

source

Callbacks

Callbacks are called after every time step when solving the BTPDE.

SpinDoctor.Printer Type

Printer(; nupdate = 1, verbosity = 2)

Print time stepping information to the console every nupdate time step (if verbosity >= 2).

source

SpinDoctor.VTKWriter Type

VTKWriter(; nupdate = 1, dir = "output/time_stepping", filename = "solution")

Write the magnetization to VTK files every nupdate time step. The files are collected in the ParaView data collection file dir/filename.pvd, which may be opened in ParaView. The compartments are separate blocks.

source

SpinDoctor.Plotter Type

Plotter(; nupdate = 1)

Plot the evolution of the BTPDE during time stepping. This requires loading a Makie plotting backend (using GLMakie for an interactive window, or using CairoMakie). The plot is updated every nupdate time step. The resulting figure contains a plot of the time profile, total signal attenuation, and magnetization field (complex magnitude and phase shift), and is available as plotter.fig after solving. It is shown in a window during the time stepping with GLMakie.

source

Plotting

The plotting functions require a Makie backend, e.g. using CairoMakie or using GLMakie.

SpinDoctor.plot_surfaces Function

plot_surfaces(surfaces; boundaries = 1:nboundary, title = "", size = (600, 500))

Plot the surface triangulation surfaces (as returned by create_geometry), with one color per boundary. Requires a Makie backend.

source

SpinDoctor.plot_mesh Function

plot_mesh(femesh; title = "", decorations = true, size = (600, 500), kwargs...)
plot_mesh!(ax, femesh; compartments, colors, alpha = 1, strokewidth, strokecolor)

Plot the finite element mesh femesh, with one color per compartment: the elements in 2D, and the surfaces of the compartments in 3D (each interface once). Only the given compartments are shown. colors and alpha (transparency) may be given for each compartment (alpha also as one value). With strokewidth > 0, the edges of the elements or facets are drawn.

plot_mesh creates a figure (without axis decorations if decorations = false), and plot_mesh! plots into the axis ax (e.g. an Axis3 with aspect = :data). Requires a Makie backend.

source

SpinDoctor.plot_mesh! Function

plot_mesh!(ax, femesh; compartments, colors, alpha = 1, strokewidth, strokecolor)

Plot the finite element mesh femesh into the axis ax, with the keywords of plot_mesh.

source

SpinDoctor.plot_field Function

plot_field(femesh, u; part = :abs, title = "", decorations = true, size, kwargs...)
plot_field!(ax, femesh, u; part = :abs, compartments, colormap, colorrange)

Plot the complex field u (e.g. the magnetization) on the finite element mesh: its magnitude (part = :abs), real or imaginary part (:real, :imag) or phase (:phase), on the elements in 2D and on the surfaces of the compartments in 3D. u may be an Observable, for animations.

plot_field creates a figure with a colorbar, and plot_field! plots into the axis ax and returns the plot (e.g. for Colorbar(position, plot)). Requires a Makie backend.

source

SpinDoctor.plot_field! Function

plot_field!(ax, femesh, u; part = :abs, compartments, colormap, colorrange)

Plot the complex field u into the axis ax, with the keywords of plot_field. Returns the plot (e.g. for Colorbar(position, plot)).

source

SpinDoctor.plot_hardi Function

plot_hardi(directions, values; title = "", label = "", decorations = true, size, colormap)
plot_hardi!(ax, directions, values; colormap = :viridis, colorrange)

Plot values (e.g. signals or ADCs) in the given directions (columns of unit vectors) as a surface (3D) or a closed curve (2D) at the distances values from the origin, colored by the values. Requires a Makie backend.

source

SpinDoctor.plot_hardi! Function

plot_hardi!(ax, directions, values; colormap = :viridis, colorrange)

Plot values in the given directions into the axis ax (see plot_hardi).

source

Public, not exported

Meshes and finite elements

SpinDoctor.FEMesh Type

FEMesh(; point_map, points, facets, elements)

Finite element mesh split into compartments, with one P1 (linear) simplex mesh per compartment:

• points[i]: coordinates of the points of compartment i (size dim, npoint);

• elements[i]: point indices of its elements (size dim + 1, nelement);

• facets[i, j]: point indices of the facets of compartment i on boundary j (size dim, nfacet, empty if compartment i does not touch boundary j);

• point_map[i]: indices of the points of compartment i in the global mesh (sorted).

Points on an interface are duplicated, one copy for each compartment.

source

SpinDoctor.split_mesh Function

split_mesh(mesh_all)

Split a global mesh into compartments. The global mesh is a named tuple with points, facets and facetmarkers (boundary number of each facet), and elements and elementmarkers (compartment number of each element), as returned by create_mesh.

source

SpinDoctor.tagged_mesh Function

tagged_mesh(points, elements, elementmarkers)

Finite element mesh of a simplex mesh created elsewhere, with the compartment number of each element (elementmarkers, from 1 to the number of compartments), e.g. the cell tags of a mesh from Gmsh or FEniCS. The boundaries are found from the elements: the interfaces between pairs of compartments, and the outer boundary of each compartment, numbered as in the setup recipes (see create_surfaces).

source

SpinDoctor.facet_elements Function

facet_elements(elements, facets)

For each facet (column of facets), find an element (column of elements) that has the facet as one of its faces. The result is zero for facets that are not a face of any element.

source

SpinDoctor.element_geometry Function

volumes, centers = element_geometry(points, elements)

Volumes (areas in 2D) and centers of the elements of a mesh with points.

source

SpinDoctor.facet_geometry Function

areas, centers, normals = facet_geometry(points, facets)
areas, centers, normals = facet_geometry(points, facets, elements)

Areas (lengths in 2D), centers and unit normals of the facets of a mesh with points. The orientation of the normals is given by the order of the facet points, unless the elements of a compartment with the facets on its boundary are given: then the normals point out of the compartment.

source

SpinDoctor.periodic_restriction Function

P, masters = periodic_restriction(points, facets, box, dims; atol)

Identify the points on opposite faces of the periodic box = (lower, upper) along the axes dims: the points of the boundary facets on an upper face are images of the points of the facets on the lower face, translated by the box size (up to atol). Points that touch a face without a facet on it (e.g. the vertex of a cell) are not identified.

Returns the prolongation P (size npoint, nreduced), which copies the values of the reduced (periodic) points to all their images, and the index of one point for each reduced point (masters).

An error is thrown if a point on a face has no image on the opposite face.

source

SpinDoctor.assemble_mass_matrix Function

assemble_mass_matrix(elements, points; elementweights = nothing, nodeweights = nothing[, pattern])

Assemble the P1 mass matrix of the simplices elements (columns of point indices) with points (columns of coordinates). The integrand may be weighted by a piecewise constant function (elementweights, one value per element) or by a P1 function (nodeweights, one value per point).

The elements may be of lower dimension than the points, e.g. surface triangles in 3D.

source

SpinDoctor.assemble_stiffness_matrix Function

assemble_stiffness_matrix(elements, points, D = I[, pattern])

Assemble the P1 stiffness matrix    of the simplices elements with points, for a constant diffusion tensor D, or one tensor per element (a vector of matrices).

source

SpinDoctor.assemble_derivative_matrix Function

assemble_derivative_matrix(elements, points, d, D = I[, pattern])

Assemble the P1 matrix   of the simplices elements with points, for the axis d and a constant diffusion tensor D, or one tensor per element (a vector of matrices).

source

SpinDoctor.assemble_flux_matrix Function

assemble_flux_matrix(model; symmetrical = false)

Assemble the global flux matrix Q of the permeable interfaces between compartments and of the relaxing outer boundaries.

At an interface with permeability kappa between compartments with initial spin densities rho1 and rho2, the coefficients on each side are 2 rho2 / (rho1 + rho2) kappa and 2 rho1 / (rho1 + rho2) kappa, so that the equilibrium magnetization has no net flux. With symmetrical = true, kappa is used on both sides.

source

Geometry generation and files

SpinDoctor.AbstractSetup Type

AbstractSetup{T}

Setup recipe describing a geometry, from which create_geometry creates cells, surfaces and a finite element mesh. All lengths are in µm.

source

SpinDoctor.AbstractECS Type

AbstractECS

Extracellular space (ECS) around the cells of a setup recipe.

source

SpinDoctor.create_cells Function

create_cells(setup; rng = Random.default_rng())

Create geometrical configuration of cells. Return nothing or mathematical description (radii, centers). Random cell positions are drawn from rng.

source

SpinDoctor.create_surfaces Function

create_surfaces(setup, cells)

Create the surface triangulation of a setup with the cell configuration cells (for neurons: the file name of the volume mesh), to be meshed by create_mesh.

Returns a named tuple containing:

• points: the points (size dim, npoint);

• facets: the triangles (edges in 2D) (size dim, nfacet);

• facetmarkers: the boundary number of each facet;

• regions: points inside the regions (size dim, npoint), by default one per region;

• regionmarkers (optional): the region number of each of these points, for regions with several parts;

• box and meshsize (samples of periodic tissue): the box (lower, upper), whose opposite faces have the same points, and the size of the facets. The mesh keeps the facets on the box (in 2D) or all the facets (in 3D), so that the mesh is periodic.

The boundaries between n regions are numbered with the interfaces between regions i < j first (row by row: 12, 13, ..., 1n, 23, ...), followed by the outer boundaries of the regions 1, ..., n.

source

SpinDoctor.create_mesh Function

create_mesh(surfaces, refinement)

Create a global finite element mesh of the regions enclosed by the surface triangulation surfaces (see create_surfaces), with TetGen in 3D and Triangle in 2D. The elements have quality constraints, and a maximum volume (area in 2D) refinement if it is finite. For samples of periodic tissue (with a box), the facets on the box (in 2D) or all the facets (in 3D) are kept as they are, so that the mesh is periodic.

Returns a named tuple with points, facets and facetmarkers (boundary numbers), and elements and elementmarkers (region numbers).

source

SpinDoctor.deform_domain! Function

deform_domain!(points, bend, twist)

Bend the points of a 3D mesh in the x-direction (a parabola in z, with a deflection of 30 * bend times the width at the top and bottom), and twist them around the z-axis by a total angle twist over the height.

source

SpinDoctor.radial_dimension Function

radial_dimension(setup)

Dimension of the radial direction of plates (1), disks (2) and spheres (3).

source

SpinDoctor.fibonacci_sphere Function

fibonacci_sphere(npoint; half = false)

Create npoint evenly distributed points on the unit sphere (or on the upper hemisphere z > 0 if half is true).

source

SpinDoctor.read_gmsh Function

points, elements = read_gmsh(filename)

Read the points and tetrahedra of a volume mesh in the Gmsh format (version 2, as in the RealNeuronMeshes repository). Other elements (e.g. surface triangles) are skipped.

source

SpinDoctor.read_cells Function

read_cells(filename)

Read a cell configuration (; centers, radii), with the fields box and gap if present in the file.

source

SpinDoctor.save_cells Function

save_cells(filename, cells)

Save a cell configuration (; centers, radii), and its fields box and gap if present (e.g. from periodic_packing or VoronoiSetup, whose cells have no radii). The file contains the spatial dimension: 2 for disks and cylinders, 3 for spheres.

source

SpinDoctor.read_surfaces Function

read_surfaces(filename)

Read a surface triangulation from the files filename.node and filename.poly.

source

SpinDoctor.save_surfaces Function

save_surfaces(filename, surfaces)

Save a surface triangulation to the files filename.node and filename.poly, which may be passed to TetGen.

source

SpinDoctor.read_mesh Function

read_mesh(filename)

Read a global mesh from the TetGen output files filename.node, filename.face and filename.ele.

source

SpinDoctor.save_mesh Function

save_mesh(filename, mesh_all)

Save a global mesh in the TetGen output format (filename.node, filename.face and filename.ele).

source

Gradient sequences

SpinDoctor.AbstractTimeProfile Type

AbstractTimeProfile{T}

Time profile f(t) of a gradient sequence (normalized to a maximum of one).

A time profile implements a calling method returning its value at time t, and intervals, returning the characteristic times at which the profile may be discontinuous (the last one being the echo time). It can optionally implement integral and int_F2, which are otherwise computed numerically.

source

SpinDoctor.AbstractGradient Type

AbstractGradient{T}

Magnetic field gradient sequence g(t) (a vector).

source

SpinDoctor.integral Function

integral(f, t = echotime(f))

Integral F(t) of the time profile f between 0 and t.

Closed forms are used for PGSE, DoublePGSE, CosOGSE and SinOGSE. Otherwise the integral is computed numerically.

source

SpinDoctor.int_F2 Function

int_F2(f)

Integral   of the squared integral F of the time profile f. It relates the gradient amplitude g to the b-value: b = gamma^2 * g^2 * int_F2(f).

source

SpinDoctor.wavevector Function

wavevector(gradient, t, gamma)

Wave vector    of the gradient sequence gradient at time t, for the gyromagnetic ratio gamma.

source

SpinDoctor.intervals Function

intervals(f)

Characteristic times of the time profile (or gradient) f: the profile is smooth between two consecutive times. The first time is zero and the last one is the echo time.

source

SpinDoctor.echotime Function

echotime(f)

Echo time of the time profile (or gradient) f, which is the end of the last characteristic interval.

source

SpinDoctor.isconstant Function

isconstant(f)

Return true if the time profile f is interval-wise constant.

source

SpinDoctor.diffusion_time Function

diffusion_time(f)

Diffusion time of the time profile f.

source

SpinDoctor.diffusion_time_sta Function

diffusion_time_sta(f)

Diffusion time of the time profile f in the short-time approximation of the ADC (see compute_adc_sta).

source

SpinDoctor.j_integral Function

j_integral(f, lambda)

Compute the integral

for the time profile f and Laplace eigenvalue lambda.

source

Callbacks

SpinDoctor.AbstractCallback Type

AbstractCallback

Callback to be performed during the time stepping of the BTPDE. The callback is initialized before time stepping with initialize!(callback, problem, gradient, u, t), updated after every step with update!(callback, problem, gradient, u, t), and finalized after time stepping with finalize!(callback), where problem = (; model, matrices) and u is the magnetization at time t.

source