Defining new sequences
Examples of sequence definitions can be found here. These are the actual sequence definitions used for the pre-implemented sequences in MRIBuilder. In this example, we will look at the file for DiffusionSpinEcho
as an example.
We can see that there are several steps. First we define the sequence as a specific sub-type of Sequence
:
const DiffusionSpinEcho = Sequence{:DiffusionSpinEcho}
Note the duplication of the sequence name (DiffusionSpinEcho
). This name will have to be unique.
The next step is to define the sequence constructor:
function DiffusionSpinEcho(; scanner=DefaultScanner, parameters..., vars...)
build_sequence(scanner) do
(g1, g2) = dwi_gradients(...)
seq = Sequence([
:excitation => excitation_pulse(...),
nothing,
:gradient => g1,
nothing,
:refocus => refocus_pulse(...),
nothing,
g2,
nothing,
:readout => readout_event(...),
nothing,
], name=:DiffusionSpinEcho, vars...)
add_cost_function!(seq[6].duration + seq[7].duration)
return seq
end
end
The crucial part here are the individual parts used to build the sequence, defined as a vector:
[
:excitation => excitation_pulse(...),
nothing,
:gradient => g1,
nothing,
:refocus => refocus_pulse(...),
nothing,
:gradient2 => g2,
nothing,
:readout => readout_event(...),
nothing
]
We can see that this sequence is built in order by an excitation pulse, some unknown dead time (indicated by nothing
), a gradient, more dead time, a refocus pulse, more dead time, another gradient, more dead time, and finally a readout.
Some of these components have been given specific names (e.g., :excitation => ...
). This is optional, but can be useful to refer to individual components. There are helper functions available to create these components.
After creating the sequence object, we can now add secondary objectives to the cost function (using add_cost_function!
). In this example, we have:
add_cost_function!(seq[6].duration + seq[7].duration)
If we check the order of the sequence component, we see that this minimises the sum of the duration of the second gradient and the wait block before this gradient. This cost function has been added to maximise the time between the second gradient and the readout (and hence minimise the effect of eddy currents on the readout). Note that this is a secondary cost function that will only take effect if it does not interfere with any user-defined constraints and cost functions (see sequence optimisation). Some secondary cost functions will be automatically defined for you within the individual components (e.g., a trapezoidal gradient has a secondary cost function to maximise the slew rate). There is even a tertiary cost function, which minimises the total sequence duration.
The next step is to define summary variables that the user can constrain when setting up a specific instance of this sequence:
@defvar begin
diffusion_time(ge::DiffusionSpinEcho) = start_time(ge, :gradient2) - start_time(ge, :gradient)
echo_time(ge::DiffusionSpinEcho) = 2 * (variables.effective_time(ge, :refocus) - variables.effective_time(ge, :excitation))
end
For this sequence, we define the variables.diffusion_time
as the time between the start of the first and second gradient pulse, and the variables.echo_time
as twice the time between the refocus and excitation pulses. These variables need to be defined within a @defvar
block.
In addition to these sequence-specific summary variables, there are also a lot of variables already pre-defined on individual components, such as the variables.slice_thickness
of the RF pulse or the variables.gradient_strength
of the gradient pulses. To access these summary variables on a sequence-level, we need to tell MRIBuilder for which RF pulses/gradients we are interested in computing these variables:
get_pulse(seq::DiffusionSpinEcho) = (excitation=seq[:excitation], refocus=seq[:refocus])
get_gradient(seq::DiffusionSpinEcho) = (gradient=seq[:gradient], gradient2=seq[:gradient2])
get_readout(seq::DiffusionSpinEcho) = seq.readout
Note that we can indicate we are interested in multiple RF pulses/gradients by supplying them as a named tuple ((excitation=..., refocus=...)
).
Setting this allows us to get RF pulse or gradient-specific properties by calling the sequence, for example:
using MRIBuilder
sequence = DiffusionSpinEcho(bval=1., TE=:min, slice_thickness=2.)
variables.flip_angle(sequence)
(excitation = 89.99999999999999, refocus = 180.0000000000496)
Here we can see that we get the variables.flip_angle
for each of the two RF pulses defined using get_pulse
above.
The final component to defining summary variables is to define one or more default coherence pathways using get_pathway
:
get_pathway(seq::DiffusionSpinEcho) = Pathway(seq, [90, 180], 1, group=:diffusion)
Here the coherence Pathway
sets out what a specific set of spins might experience during the sequence. In this case the sequence experiences two RF pulses and is excited by the first pulse and flipped by the second ([90, 180]
). It is then observed during the first readout (1
). For such a Pathway
we can compute:
- the time that the spin spends in each longitudinal and transverse direction, which is particularly useful in the transverse direction to compute the amount of T2-weighting (
variables.duration_transverse
) and the amount of time spent dephasing to compute the amount of T2*-weighting (variables.duration_dephase
). - the diffusion weighting experienced (
variables.bval
,variables.bmat
, andvariables.net_dephasing
)
By defining a default pathway for the sequence, the user can now put constraints on any or all of these variables.
Component helper functions
MRIBuilder.Parts.HelperFunctions.dwi_gradients
— Methoddwi_gradients(; type, optimise=false, refocus=true, orientation=[1, 0, 0], group=:diffusion, variables...)
Returns two diffusion-weighting gradients that are guaranteed to cancel each other out.
Parameters
type
: A symbol describing the type of gradient. One of::trapezoid
: Pulsed trapezoidal gradient profile. SeeTrapezoid
for the relevantvariables
.:instant
: instantaneous gradient (e.g., to test short-pulse approximations). SeeInstantGradient
for the relevantvariables
.
optimise
: Whether to optimise this readout event in isolation from the rest of the sequence. Use this with caution. It can speed up the optimisation (and for very complicated sequences make it more robust), however the resulting parameters might not represent the optimal solution of any external constraints (which are ignored if the readout is optimised in isolation).scanner
: Used for testing. Do not set this parameter at this level (instead set it for the total sequence usingbuild_sequence
).
MRIBuilder.Parts.HelperFunctions.excitation_pulse
— Methodexcitation_pulse(; parameters...)
Create an excitation RF pulse.
By default there is no slice-selective gradient. To enable slice selection min_slice_thickness
has to be set to a number or to :min. If min_slice_thickness
is not set or is set to :min
, then either bandwidth
or duration
should be set, otherwise the optimisation might be unconstrained (ignore this for shape=:instant
).
Parameters
For an InstantPulse
(i.e., shape=:instant
), only the flip_angle
, phase
, and group
will be used. All other parameters are ignored.
optimise
: set to true to optimise this RF pulse separately from the embedding sequence.
Pulse parameters
shape
: The shape of the RF pulse. One of:sinc
(forSincPulse
),:constant
/:hard
(forConstantPulse
), or:instant
(forInstantPulse
). Default is :sinc for slice-selective pulses or :instant otherwise.flip_angle
: size of the flip due to the RF pulse in degrees (default: 90).phase
: angle of the RF pulse in the x-y plane in degrees (default: 0).frequency
: frequency of the RF pulse relative to the Larmor frequency in kHz (default: 0).bandwidth
: width of the RF pulse in Fourier space in kHz (default: free variable).duration
: duration of the RF pulse in ms (default: free variable).Nzeros
: sets the number of zero crossings for aSincPulse
(default: 3). Can be set to a tuple of two numbers to set a different number of zero crossings before and after the pulse maximum.group
: name of the group of which the RF pulse is part. This is used to add transformations after the sequence is optimised.
Slice selection
slice_thickness
: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to:min
to indicate that this should be minimised given the scanner constraints and user values forbandwidth
orduration
.rephase
: set to false to disable the spin rephasing after the RF pulse.rotate_grad
: name of the parameter with which the slice selection gradient will be rotated after sequence optimisation (default::FOV
).scanner
: overrides theglobal_scanner
for this part of the sequence. Recommended to set only if not part of a larger sequence.
MRIBuilder.Parts.HelperFunctions.gradient_spoiler
— Methodgradient_spoiler(; optimise=false, orientation=[0, 0, 1], rotate=:FOV, scale=:spoiler, spoiler=1., duration=:min, variables...)
Returns two DWI gradients that are guaranteed to cancel each other out.
Parameters
orientation
: Orientation of the gradient (default: slice-select direction).rotate
: in which coordinate system is theorientation
defined (default: :FOV).scale
: variable controlling how the gradients should be scaled after optimsation (default: :spoiler).optimise
: Whether to optimise this readout event in isolation from the rest of the sequence. Use this with caution. It can speed up the optimisation (and for very complicated sequences make it more robust), however the resulting parameters might not represent the optimal solution of any external constraints (which are ignored if the readout is optimised in isolation).scanner
: Used for testing. Do not set this parameter at this level (instead set it for the total sequence usingbuild_sequence
).
Variables
variables.spoiler
: maximum spoiler scale (before applying any reductions due toscale
).- Any other parameters expected by
Trapezoid
.
MRIBuilder.Parts.HelperFunctions.interpret_image_size
— Methodinterpret_image_size(fov, resolution, voxel_size, slice_thickness)
Combines the user-provided information of the image size to produce a tuple with:
slice_thickness
: if not set explicitly, will be set to the third element ofvoxel_size
resolution_z
: number of voxels in the slice-select direction.- keywords parameters for the
readout_event
functions with thefov
,resolution
, andvoxel_size
in the x- and y- dimensions.
MRIBuilder.Parts.HelperFunctions.readout_event
— Methodreadout_event(; type, optimise=false, variables...)
Adds a readout event to the sequence.
Parameters
type
: A symbol describing the type of readout. It will default to:epi
if a resolution has been set and:instant
otherwise. Can be set to one of the following::epi
: EPI readout. SeeEPIReadout
for the relevantvariables
.:instant
: single isolated readout eventSingleReadout
(e.g., for NMR). Does not expect anyvariables
.
optimise
: Whether to optimise this readout event in isolation from the rest of the sequence. Use this with caution. It can speed up the optimisation (and for very complicated sequences make it more robust), however the resulting parameters might not represent the optimal solution of any external constraints (which are ignored if the readout is optimised in isolation).scanner
: Used for testing. Do not set this parameter at this level (instead set it for the total sequence usingbuild_sequence
).
MRIBuilder.Parts.HelperFunctions.refocus_pulse
— Methodrefocus_pulse(; parameters...)
Create an excitation RF pulse.
By default there is no slice-selective gradient. To enable slice selection slice_thickness
has to be set to a number or to :min. If slice_thickness
is not set or is set to :min
, then either bandwidth
or duration
should be set, otherwise the optimisation might be unconstrained (ignore this for shape=:instant
).
Parameters
optimise
: set to true to optimise this RF pulse separately from the embedding sequence.
Pulse parameters
For an InstantPulse
(i.e., shape=:instant
), only the flip_angle
, phase
, and group
will be used. All other parameters are ignored.
shape
: The shape of the RF pulse. One of:sinc
(forSincPulse
),:constant
/:hard
(forConstantPulse
), or:instant
(forInstantPulse
). Default is :sinc for slice-selective pulses or :instant otherwise.flip_angle
: size of the flip due to the RF pulse in degrees (default: 180).phase
: angle of the RF pulse in the x-y plane in degrees (default: 0).frequency
: frequency of the RF pulse relative to the Larmor frequency in kHz (default: 0).bandwidth
: width of the RF pulse in Fourier space in kHz (default: free variable).duration
: duration of the RF pulse in ms (default: free variable).Nzeros
: sets the number of zero crossings for aSincPulse
(default: 3). Can be set to a tuple of two numbers to set a different number of zero crossings before and after the pulse maximum.group
: name of the group of which the RF pulse is part. This is used to add transformations after the sequence is optimised.
Slice selection and spoilers
slice_thickness
: minimum slice thickness that should be possible without adjusting the sequence timings in um (not mm!) (default: no slice selection). Can be set to:min
to indicate that this should be minimised given the scanner constraints and user values forbandwidth
orduration
.spoiler
: set to the spatial scale on which the spins should be dephased in mm. For rotating spoilers, this does include the contribution from the slice select gradient as well.rotate_grad
: name of the parameter with which the slice selection and spoiler gradient will be rotated after sequence optimisation (default::FOV
).scanner
: overrides theglobal_scanner
for this part of the sequence. Recommended to set only if not part of a larger sequence.orientation
: vector with orientation of slice select gradient and pulses (defaults: z-direction).
MRIBuilder.Parts.HelperFunctions.saturation_pulse
— Methodsaturation_pulse(; pulse=(), binomial_order=1, variables...)
Creates a saturation pulse consisting of repeating instances of pulse_type
.
Parameters
pulse
: Parameters/variables used to define the basic RF pulse used in the saturation pulse. These can be:type
:ConstantPulse
(default),SincPulse
,InstantPulse
, etc.- Variables describing that type (e.g.,
flip_angle
,duration
,frequency
)
binomial_order
: How many repeats ofpulse_type
there should be in a single binomial pulse.
Variables
nrepeat
: number of repeats of the binomial pulse. This should typically be set by the user.interblock_delay
: time between the start of each pulse. Should be short compared with the T1, but long enough to limit SAR.duration
: total duration of the saturation pulse.
Optimisation helper functions
MRIBuilder.BuildSequences.build_sequence
— MethodWrapper to build a sequence.
Use as
build_sequence(scanner;) do
...
end
Within the code block you can create one or more sequences, e.g.
seq = Sequence(
SincPulse(flip_angle=90, phase=0, duration=2., bandwidth=:max)
nothing.,
SingleReadout
)
You can also add any arbitrary constraints or objectives using one of:
set_simple_constraints!
apply_simple_constraint!
add_cost_function!
As soon as the code block ends the sequence is optimised (if optimise=true
) and returned.
Parameters
scanner
: Set to aScanner
to limit the gradient strength and slew rate. When this call tobuild_sequence
is embedded in another, this parameter can be set tonothing
to indicate that the same scanner should be used.optimise
: Whether to optimise and fix the sequence as soon as it is returned. This defaults totrue
if a scanner is provided andfalse
if no scanner is provided.n_attempts
: How many times to restart the optimiser (default: 20)? Decrease if you want to quickly check feasibility of a simple sequence. Increase if convergence fails for a complex sequence.kwargs...
: Other keywords are passed on as attributes to theoptimiser_constructor
(e.g., setprint_level=3
to make the Ipopt optimiser quieter).
MRIBuilder.BuildSequences.fixed
— Methodfixed(building_block)
Return an equiavalent BuildingBlock
with all free variables replaced by numbers.
This will only work after calling optimize!
(@ref)(global_model
()). It is used internally by build_sequence
.
MRIBuilder.BuildSequences.get_optimiser
— Methodget_optimiser(model)
Returns a JuMP solver (https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) appropriate for this model.
MRIBuilder.BuildSequences.global_scanner
— MethodMRIBuilder.BuildSequences.iterate_cost
— Methoditerate_cost()
Return a sequence of all the cost functions that should be optimised in order.
MRIBuilder.BuildSequences.model_has_integers
— Methodmodel_has_integers(model)
Returns true if the model contains integer variables.
Pathways API
MRIBuilder.Pathways.GradientTracker
— TypeHelper structure for PathwayWalker
, which is itself a helper for Pathway
. You are deep down the rabit hole now...
For documentation, see that structure and walk_pathway!
and update_walker_gradient!
.
MRIBuilder.Pathways.Pathway
— TypePathway(sequence::Sequence, pulse_effects::Vector{:Symbol/Number}, readout_index=1; group=nothing)
Describes how a specific spin/isochromat might experience the sequence.
Only a single pathway through the RF pulses is considered, so that at every point in time the spins are in one of the following four states:
- +longitudinal: initial relaxed state
- +transverse: excited state. During this time gradients will affect the
variables.area_under_curve
(orvariables.qval
) andvariables.bval
. - -longitudinal: inverse state
- -transverse: inverse excited state. During this time all gradients will have the inverse effect compared with +transverse.
The RF pulses cause mappings between these different states as described below.
Parameters
sequence
: MRISequence
to be considered.pulse_effects
: How each RF pulse affects the spins. This can be one of the following::skip
/:ignore
/0: This RF pulse leaves the spins unaffected.:refocus
/:invert
/180: Flips the sign of the spin state (i.e., +longitudinal <-> -longitudinal, +transverse <-> -transverse):excite
/90: Takes spin state one step along the following sequence +longitudinal -> +transverse -> -longitudinal -> -transverse -> +longitudinal:neg_excite
/270/-90: Inverse step compared with:excite
.
readout_index
: After encountering the number of pulses as defined inpulse_effects
, continue thePathway
until the readout given byindex
is reached. If set to 0 thePathway
is terminated immediately after the last RF pulse.group
: which gradient grouping to consider for thenet_dephasing
andbmat
. If not set, all gradients will be considered (using their current alignment).
Attributes
Over the pathway the following values are computed. Each can be accessed by calling the appropriate function:
Timings
variables.duration_state
: The total amount of time spent in a specific state in this pathway (+longitudinal, +transverse, -longitudinal, or -transverse)variables.duration_transverse
: The total amount of time the spins spent in the transverse plane in ms. This can be used to quantify the expected effect of T2-decay.variables.duration_dephase
: The total amount of time the spins spent in the +transverse relative to -transverse state in ms. The absolute value of this can be used to quantify the expected effect of T2'-decay.
Effect of gradients
The area under curve, q-values, and b-values are computed separately for each group of gradients (depending on the group
keyword set during construction).
variables.net_dephasing
: Net displacement vector in k-space/q-space.variables.bmat
: Net diffusion weighting due to gradients along thePathway
in matrix form.variables.bval
: Net diffusion weighting due to gradients along thePathway
as a single number.
MRIBuilder.Pathways.PathwayWalker
— TypeHelper structure for Pathway
.
For documentation, see that structure and walk_pathway!
.
MRIBuilder.Pathways.duration_state_index
— Methodduration_state_index(transverse, positive)
Returns the index of a specific state in the Pathway.duration_state
vector.
This function is used internally to access that vector.
MRIBuilder.Pathways.interpret_pulse_effects
— Methodinterpret_pulse_effects(number_or_symbol)
Interpret the various numbers and symbols that can be passed on to a Pathway
.
The result will be one of:
- :ignore (if input is 0, :ignore, or :skip).
- :excite (if input is 90 or :excite).
- :refocus (if input is 180, :refocus, or :excite).
- :neg_excite (if input is -90, 270, or :negexcite).
MRIBuilder.Pathways.update_gradient_tracker_till_time!
— Methodupdate_gradient_tracker_till_time!(walker::PathwayWalker, key, new_time)
update_gradient_tracker_till_time!(tracker::GradientTracker, new_time)
Update the bmat
for any time passed since the last update (assuming there will no gradients during that period).
The bmat
is updated with the outer produce of qvec
with itself multiplied by the time since the last update.
When called with the first signature the tracker will be created from scratch if a tracker with that key
does not exist.
MRIBuilder.Pathways.update_walker_gradient!
— Methodupdate_walker_gradient!(gradient_block::GradientWaveform, walker::PathwayWalker, gradient_start_time::VariableType; overlapping_pulses=[], overlapping_readouts=[])
Update the walker's qvec
and bmat
based on the given gradient_block
.
The following steps will be taken:
- Do nothing if
walker.is_transverse
is false - increase the appropriate
walker.bmat
by the outer product ofwalker.qvec
with itself multiplied by the time since the last gradient - update the appropriate
walker.qvec
andwalker.bmat
based on the gradient waveform. This will require appropriateqvec
/bmat
functions to be defined for the gradient building block. - update
walker.last_gradient_time
to the time at the end of the gradient.
This requires variables.bmat_gradient
and variables.qvec
to be implemented for the GradientWaveform
.
MRIBuilder.Pathways.update_walker_instant_gradient!
— Methodupdate_walker_instant_gradient!(gradient, walker)
MRIBuilder.Pathways.update_walker_pulse!
— Methodupdate_walker_pulse!(walker::PathwayWalker, pulse_effects::Vector, pulse_time)
Apply the first element of pulse_effects
to the walker
at the given pulse_time
.
The following steps will be taken if the first pulse_effect
is not :ignore
- if
walker.is_transverse
is true before the pulse, increase thewalker.bmat
by the outer product ofwalker.qvec
with itself multiplied by the time since the last gradient - update
walker.duration_states
with time since last pulse. - update
walker.last_pulse_time
- update
walker.is_transverse
, andwalker.is_positive
based on the first value inpulse_effects
. - if
walker.is_positive
changed in the previous step than thewalker.qvec
needs to be flipped. - remove the first element from
pulse_effects
.
If the first element is :ignore
the only effect is that the first element is removed from pulse_effects
.
MRIBuilder.Pathways.update_walker_till_time!
— Functionupdate_walker_till_time!(walker::PathwayWalker, new_time[, gradient_group])
Updates all parts of a PathwayWalker
up to the given time.
This updates the walker.duration_states
and the bmat
for each gradient tracker. If gradient_group
are provided, then only the gradient tracker matching that group will be updated. If that gradient tracker does not exist, it will be created.
This function is used to get the walker
up to date till the start of a gradient, pulse, or final readout.
MRIBuilder.Pathways.walk_pathway!
— Methodwalk_pathway!(bb::Sequence/BuildingBlock, walker::PathwayWalker, pulse_effects::Vector, nreadout::Ref{Int}, start_time)
Computes the effect of a specific ContainerBlock
(starting at start_time
) on the PathwayWalker
.
For individual pulses and gradients, the following behaviour is implemented:
- If a pulse is encountered, call
update_walker_pulse!
(walker, pulse_effects, pulse_effective_time)
- If a gradient is encountered, call
update_walker_gradient!
(gradient, walker, gradientstarttime)
For overlapping gradients/pulses, one should first encounter the part of the gradient before the variables.effective_time
of the pulse, then apply the pulse, and then the rest of the gradient. This effective time can be passed on to update_walker_gradient!
to allow part of the gradient waveform to be applied.
The function should return true
if the Pathway
has reached its end (i.e., the final readout) and false
otherwise.
MRIBuilder.Variables.get_pathway
— Functionget_pathway(sequence)
Gets the main Pathway
that spins are expected to experience in the sequence.
Multiple pathways might be returned as an array or (named)tuple.
MRIBuilder.Variables.variables.area_under_curve
— Functionarea_under_curve(pathway::Pathway)
Return net displacement in k-space (i.e., spoiling) experienced by the spins following a specific Pathway
.
Only gradients active while the spins are in the transverse plane are considered.
Returns a NamedTuple with the area_under_curve
for all gradient groups.
MRIBuilder.Variables.variables.bmat
— Functionbmat(pathway::Pathway)
Return 3x3 diffusion-weighted matrix experienced by the spins following a specific Pathway
in rad^2 ms/um^2.
Only gradients active while the spins are in the transverse plane are considered.
Returns a NamedTuple with the bmat
for all gradient groups.
MRIBuilder.Variables.variables.bval
— Functionbval(pathway::Pathway)
Return size of diffusion-weighting experienced by the spins following a specific Pathway
in rad^2 ms/um^2.
Only gradients active while the spins are in the transverse plane will contribute to the diffusion weighting.
Returns a NamedTuple with the bval
for all gradient groups.
MRIBuilder.Variables.variables.duration_dephase
— Functionduration_dephase(pathway::Pathway)
Returns the net time that spins following the given Pathway
spent in the +transverse versus the -transverse state. This determines the amount of T2'-weighting as $e^{t/T_2'}$, where $t$ is the duration_dephase
.
Also see variables.duration_transverse
for T2-weighting.
MRIBuilder.Variables.variables.duration_state
— Functionduration_state(pathway::Pathway, transverse::Bool, positive::Bool)
Returns how long the Pathway
spent in a specific state.
The requested state can be set using transverse
and positive
as follows:
transverse=false
,positive=true
: +longitudinaltransverse=true
,positive=true
: +transversetransverse=false
,positive=false
: -longitudinaltransverse=true
,positive=false
: -transverse
MRIBuilder.Variables.variables.duration_transverse
— Functionduration_transverse(pathway::Pathway)
Returns the total amount of time that spins following the given Pathway
spent in the transverse plane. This determines the amount of T2-weighting as $e^{t/T_2}$, where $t$ is the duration_transverse
.
Also see variables.duration_dephase
for T2'-weighting.
MRIBuilder.Variables.variables.net_dephasing
— Functionnet_dephasing(pathway::Pathway)
Return net displacement vector in k-space/q-space experienced by the spins following a specific Pathway
in kHz/um.
Only gradients active while the spins are in the transverse plane are considered.
Returns a NamedTuple with the qvec
for all gradient groups.