Extending KomaMRI for Bloch–McConnell simulations via SimulationMethod

[Alena795]

Hello,

First of all, thank you very much for developing KomaMRI.

I would like to ask about the feasibility of implementing a two-pool Bloch–McConnell simulator, using the extensibility of “SimulationMethod”. In the paper “KomaMRI.jl: An open-source framework for general MRI simulations with GPU acceleration”, it is stated that “The Bloch–McConnell equations could also be implemented accurately with the operator splitting method” suggesting that a multi-pool exchange model could be possible in Koma.

From reading the Koma documentation, especially the parts mentioning:

• custom spin-state representations,
• overriding excitation/precession behavior by defining custom run_spin_precession and run_spin_excitation externally,
• and possible EPG-style extensions,

it seems that KomaMRI is intentionally designed so that the simulator is not fundamentally tied to the standard Bloch state representation (Mx, My, Mz). My understanding is that Bloch–McConnell dynamics mainly require replacing the magnetization state representation and redefining the rotation/relaxation steps with appropriate operator splitting in the precession/excitation functions. Instead of a single-pool state (Mx,My,Mz), one would evolve a multi-pool state such as: (Mxa,Mya,Mza,Mxb,Myb,Mzb).

So I suppose my question is: could the simulator architecture support a multi-pool Bloch–McConnell implementation purely through external extension using a custom SimulationMethod, or would changes in the source code be necessary?

Thank you very much!

[cncastillo]

Hi @Alena795! Thanks for considering using Koma :)

As you noticed, Koma is designed to be extended. The Bloch–McConnell equation is something we are interested in implementing later this year, after multi-receive and pTx support. We would also be happy to receive contributions from the community.

For Bloch–McConnell specifically, some heavier lifting would be needed. For example, the exchange constants, pool fractions, and other pool-specific properties would need to be added to the Phantom. One possible approach would be to create something like a Phantom2Pools, which should not be too difficult.

Everything related to the simulation itself could likely be extended without touching the core source code. For instance, we could add an Mxy2Pools spin state containing the values for pools “a” and “b”, and then define the simulation around that. My colleague @mcgrathcm has previously implemented Bloch–McConnell simulations using the operator-splitting method in CMRsim, so we do have some experience with this.

One important point is that implementing the simulation may not be the hardest part. A more challenging aspect may be defining a Phantom with meaningful pool constants, as well as defining an experiment whose results we can corroborate. Ideally, we would have some expected MRI result to compare against. We could always compare the implementation to a numerical ODE solution of the Bloch–McConnell equations, but I would prefer tests that make more sense from an MRI perspective.

If you can help us define a relevant experiment with an expected result, that would help a lot. Otherwise, I would be happy to discuss this further in a meeting and guide you through the implementation.

I am also thinking that the exchange constants, pool fractions, and other properties such as T1, T2, and off-resonance could be defined for every position. However, that level of flexibility may not be necessary for an initial implementation. We could start with a simpler version first, possibly one that does not require defining a completely new Phantom.

[Alena795]

Hi @mcgrathcm,

Thank you very much for your answer! :)
I will check the (.) operator and rename the variables, the notation is indeed confusing. For exp(D1[i]) specifically, I think it might be fine because D1 is a vector containing the two eigenvalues of the block A_t1. So exp(D1[i]) would just be the exponential of a scalar eigenvalue, not a matrix exponential

Regarding the different cases, I checked and also believe that assuming T1,T2, kab, kba all strictly positive always leads to the case with two distinct real eigenvalues (covered by case 3), since the discriminant of the characteristic polynomial of the 2*2 block is strictly positive. From what I understand, case 3 also covers the situation where the block is already diagonal with repeated eigenvalue, so discrimiant equal to 0 and two independant eigenvectors. I will try to re-write these functions without the case handling.

Koma is bit over typed at the moment. functions like simulate don't have any multi dispatch for phantom, so there is no real reason to type it. That being said, I think it would be much more clear to dispatch on sim_method, since you already have BlochMcConnell2Pools.

I just had one question concerning this: would it mean that simulate() could be defined as:

function simulate(obj, seq::Sequence, sys::Scanner; sim_params=Dict{String,Any}, ...)

without specifying a type for obj?

[mcgrathcm]

Thanks @Alena795 , and great that the cases do resolve to a single case, makes it much nicer!
For the typing, I think there are two options here. Either like your example where there is no typing for obj, or we make a new abstract type Phantom that both the default phantom and your 2pool phantom inherit, similar to how SimulationMethod is handled. I would lean towards the second, as I expect in the future there will may additional phantom types (specifically time varying relaxation). @cncastillo what do you think?

[cncastillo]

I am wondering if we can just remove the types from simulate (and many other internal functions), it shouldn't have any performance impact. As @mcgrathcm mentions, Koma internals are overtyped (because i wrote it before i knew Julia properly!), and this can generate some problems. We could change stuff to abstract type AbstractPhantom and Phantom <: AbstractPhantom and Phantom2Pool <: AbstractPhantom, but we would still need to change the internals.

Let me see if during the weeking i can help with the overtyping part so its easier to extend, as I would need to do this anyway.

Related #639.

[Alena795]

Thanks @mcgrathcm, @cncastillo ! An abstract type AbstractPhantom like SpinStateRepresentation and SimulationMethod sounds great

[Alena795]

Hi,

I finished implementing the Bloch-McConnell simulator, in this Koma fork: https://github.com/Alena795/KomaMRI.jl.git
I was working on hyperpolarized carbon 13 MRI (specifically hyperpolarized [1-13C]pyruvate to lactate conversion) and the simulation results seem to match the expected Bloch-McConnell exchange dynamics.

The experiment used to validate the Bloch-McConnell simulator consisted of dynamic hyperpolarized MRI simulations involving two metabolite pools: pyruvate and lactate. A two-region phantom was simulated, representing tumor tissue with a higher kPL and healthy tissue with a baseline kPL. The acquisition sequence consisted of repeated frequency-selective sinc excitation pulses followed by spiral readouts, generating a time series of reconstructed images for both metabolites. The pyruvate and lactate signals were then averaged within dedicated ROIs in the tumor and healthy tissue regions to obtain their temporal evolution, providing experimental Mxy(t) curves for subsequent fitting with the Bloch-McConnell kinetic model and estimation of kPL.

The results show the characteristic temporal pyruvate-lactate curves expected from the Bloch-McConnell model, with progressive pyruvate depletion and a corresponding increase, peak, and subsequent decay of the lactate signal.
The simulator is still under development with some improvements needed (particularly in terms of simulation speed), but I just wanted to share these results :).

Results_BLochMcConnell.zip