About the numerical solver of KomaMRI and T2*

[timnotavailable]

Hello KomaMRI team, I have several questions regarding your software:
1.What kind of solver are you using for the Bloch equation? is your solver adaptive like Dopri5 or fixed step?
2.How do you simulate the T2* effectively, meaning that how many spins per voxel do you need to simulate to avoid spurious signal?
3.Can I simulate the adiabatic pulse with your software?

[cncastillo]

Hi! Thanks for you interest :). I will attempt to answer below:

Simulation solver

Our current solver uses the operator splitting method, and for each sub-problem, what we do can be described as an adaptive first-order implicit method for the precession (trapezoidal integration) and adaptive zero-order explicit method for the excitation. But we will upgrade the excitation method to match the precession (see #458), as we are very interested in simulating flow and more advanced motions (check #539 (comment)).

$T_2^*$ decay

We currently do not simulate T2* effects (se issue #62). We store a T2s variable in the phantom that can be used for that purpose, but you would need to create the different isochromats by hand:

function include_T2s(obj, Niso)
	T2prime = (1 ./ obj.T2s .- 1 ./ obj.T2) .^ -1
	prob(Δw) = T2prime ./ (π .* (1 .+ T2prime .^ 2 .* Δw .^ 2 ))
	Δw_lim = 2π * 200 # -200 Hz to 200 Hz
	obj_aux = copy(obj)
    obj_with_T2s = Phantom()
    for  Δw in range(-Δw_lim, Δw_lim, Niso))
		obj_aux.Δw .= Δw
		obj_with_T2s += prob(Δw) * obj_aux
    end
	return obj_with_T2s
end

You can create a simulation method that does this by default by extending the function initialize_spins_state.

struct BlochWithT2s <: SimulationMethod
	Niso::Int
end
function KomaMRICore.initialize_spins_state(
    obj::Phantom{T}, 
    sim_method::BlochWithT2s # <-------
) where {T<:Real}
    obj = include_T2s(obj, sim_method.Niso) # <-------
    Nspins = length(obj)
    Mxy = zeros(T, Nspins)
    Mz = obj.ρ
    Xt = Mag{T}(Mxy, Mz)
    return Xt, obj
end

Both code snippets could have some errors, but that is the overall idea. Then, just set the obj.T2s to the appropriate values and simulate using sim_params["sim_method"] = BlochWithT2s(100) (for Niso = 100, but you can choose any number of isochromats).

It is not incorporated by default because I believe there is a better way to do this, and how exactly the frequencies of the isochromats are chosen could by itself be a research topic, so at the moment it is responsibility of the user to choose some appropriate method for their problem.

Adiabatic pulses

You can simulate adiabatic pulses, in fact we have a few examples that should be in the documentation soon (see #227, or https://github.com/JuliaHealth/KomaMRI.jl/tree/master/examples/3.tutorials/adiabatic_pulses). If you have a cool example you would like us to implement please let us know. Random fact, the cover of my PhD thesis uses an adiabatic pulse simulated with KomaMRI:

[timnotavailable]

Hello, thanks for the answer. Regarding the reply I still have a bit problems:

For your advanced simulation like flow, will you meet the problem that 1st order implicit method does not work anymore? In fact as far as I knew JEMRIS used Adam which utilize higher order term for simulation and this approach is better than Runge-Kutta 4th or 5th order, how do you guarantee that only with 1st order implicit method stiffness problem of Bloch equation can be solved easily?

$T2* decay$
I read the issue provided, however my concern is the T2' the same for every voxels or different voxel by voxel?

For adiabatic pulse, I'm gonna use simulate some pulses with arbitrary waveform and frequency. Predefined BIR / HS/Wurst is interesting but might not so suitable for my case, but thank you for your example of adiabatic pulses.

[cncastillo]

For your advanced simulation like flow, will you meet the problem that 1st order implicit method does not work anymore? In fact as far as I knew JEMRIS used Adam which utilize higher order term for simulation and this approach is better than Runge-Kutta 4th or 5th order, how do you guarantee that only with 1st order implicit method stiffness problem of Bloch equation can be solved easily?

We can easily include higher order integration methods if there is a problem with the current ones, but I haven't seen problems in our preliminary results. Do you have a case where this could be a concern? You can also reduce this parameter sim_params["Δt"] to reduce the maximum gradient step size, in case the automatically determined one doesn't work.

I read the issue provided, however my concern is the T2' the same for every voxels or different voxel by voxel?

T2' would in fact be different voxel by voxel.

For adiabatic pulse, I'm gonna use simulate some pulses with arbitrary waveform and frequency. Predefined BIR / HS/Wurst is interesting but might not so suitable for my case, but thank you for your example of adiabatic pulses.

You should be able to provide arbitrary AM and FM modulation. If you see a problem with the simulation accuracy, try reducing sim_params["Δt_rf"] to reduce the maximum allowable RF step size.

If you have a more specific example I could help you implement it to get started. Feel free to ask as many questions as you like! :)

[timnotavailable]

Hello , thank you for your nice reply! Everything is clear for me now!