Generalized Multi-Tissue Modeling¶

The measured PGSE diffusion signal depends on Echo Time (TE), gradient strength (G), orientation $\textbf{n}$, pulse separation $\Delta$ and pulse duration $\delta$. The signal representation can be separated in terms of amplitude and the shape:

\begin{equation}\label{eq:separation} S(G, \textbf{n}, \Delta, \delta, TE)=S_0(TE)\cdot E(G,\textbf{n}, \Delta, \delta). \end{equation}

where we can notice that the amplitude only depends on TE and the shape on all the others. In most models (NODDI, SMT, VERDICT etc.) ONLY the signal shape is fitted. Fitting the signal itself has only been explored in Multi-Tissue models like Multi-Tissue CSD.

In Dmipy, we generalize multi-tissue modeling to ANY MC, MC-SM and MC-SH model.

Setting up a Multi-Tissue model : Example Multi-Tissue Ball and Stick¶

Instantiate base models¶

In [1]:

from dmipy.signal_models import gaussian_models, cylinder_models

# setting base models
ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
models = [ball, stick]

# setting arbitrary S0 values for the example
S0_ball = 12.
S0_stick = 8.
S0_responses = [S0_ball, S0_stick]

Instantiate Multi-tissue Ball and Stick for standard MC-model¶

A multi-tissue model is created exactly as a regular one - only the S0 response values need to be given at the upon multi-compartment model instantiation.

In [2]:

from dmipy.core import modeling_framework
mt_BAS_standard = modeling_framework.MultiCompartmentModel(
    models=models,
    S0_tissue_responses=S0_responses)

Simulate some test signal¶

The multi-tissue model only differs from a standard model when fitting it to data. When simulating data it will still generate the signal attenuation.

In [3]:

from dmipy.data.saved_acquisition_schemes import wu_minn_hcp_acquisition_scheme
scheme = wu_minn_hcp_acquisition_scheme()

# generate test data.
params = {
    'G1Ball_1_lambda_iso': 3e-9,
    'C1Stick_1_mu': [0., 0.],
    'C1Stick_1_lambda_par': 1.7e-9,
    'partial_volume_0': 0.5,  # equal volume fractions as SIGNAL fractions
    'partial_volume_1': 0.5
}

# total signal intensity is 10
S0_signal = 10.
S = mt_BAS_standard.simulate_signal(scheme, params) * S0_signal
S.shape

Out[3]:

(288,)

Fit Multi-Tissue model and compare estimated fractions¶

We can fit the model as usual, as the multi-tissue optimization occurs AFTER the standard optimization. It is therefore independent and naturally follows other approaches.

In [4]:

mt_BAS_standard_fit = mt_BAS_standard.fit(scheme, S)

Using parallel processing with 8 workers.
Setup brute2fine optimizer in 0.332586050034 seconds
Fitting of 1 voxels complete in 0.553920984268 seconds.
Average of 0.553920984268 seconds per voxel.
Starting secondary multi-tissue optimization.
Multi-tissue fitting of 1 voxels complete in 0.00136590003967 seconds.

We now have access the the signal fractions based on the signal attenuation, and the non-normalized and normalized volume fractions based on the signal**.

In [5]:

sig_fracts = mt_BAS_standard_fit.fitted_parameters
vol_fracts = mt_BAS_standard_fit.fitted_multi_tissue_fractions
vol_fracts_norm = mt_BAS_standard_fit.fitted_multi_tissue_fractions_normalized

We can see that as we added equal signal fractions of the ball and stick, that indeed the estimated signal fractions are equal to each other.

In [6]:

sig_fracts

Out[6]:

{'C1Stick_1_lambda_par': array([1.69975574e-09]),
 'C1Stick_1_mu': array([[ 3.14156351, -0.65826994]]),
 'G1Ball_1_lambda_iso': array([3.e-09]),
 'partial_volume_0': array([0.50001804]),
 'partial_volume_1': array([0.49998196])}

But the non-normalized volume fractions after the secondary optimization (which does only estimated the linear volume fractions and does not impose unity) are now scaled according to the tissue-specific S0 responses:

In [7]:

vol_fracts

Out[7]:

{'partial_volume_0': array([0.4167901]),
 'partial_volume_1': array([0.62500692])}

But as it's valuable to find the normalized volume fractions, we also easily provide the normalized volume fractions. Here, we can now see that while the signal fraction is indeed equal, the volume fraction, actually in terms of signal production by diffusing particles, is in fact 0.4 to 0.6.

In [8]:

vol_fracts_norm

Out[8]:

{'partial_volume_0': array([0.40006843]),
 'partial_volume_1': array([0.59993157])}

Note here that we can only set the S0 response to a signal value for an MC-model (no voxel varying) and currently data with multiple TE (so multiple S0 responses per model) is not implemented.

Multi-Tissue Spherical Mean Modeling¶

Setting up Multi-Tissue modeling is exactly the same when setting up MC-spherical mean models:

In [9]:

mt_BAS_sm = modeling_framework.MultiCompartmentSphericalMeanModel(
    models=models,
    S0_tissue_responses=S0_responses)
mt_BAS_sm_fit = mt_BAS_sm.fit(scheme, S)

Using parallel processing with 8 workers.
emean 3.425800320368238
emean 2.3937949986066385
emean 1.9601890280398429
Setup brute2fine optimizer in 0.0430829524994 seconds
Fitting of 1 voxels complete in 0.0924401283264 seconds.
Average of 0.0924401283264 seconds per voxel.
Starting secondary multi-tissue optimization.
Multi-tissue fitting of 1 voxels complete in 0.00108790397644 seconds.

In [10]:

mt_BAS_sm_fit.fitted_parameters

Out[10]:

{'C1Stick_1_lambda_par': array([1.55262663e-09]),
 'G1Ball_1_lambda_iso': array([2.99995881e-09]),
 'partial_volume_0': array([0.51901797]),
 'partial_volume_1': array([0.48098203])}

In [11]:

mt_BAS_sm_fit.fitted_multi_tissue_fractions

Out[11]:

{'partial_volume_0': array([0.4325418]),
 'partial_volume_1': array([0.60115111])}

Multi-Tissue Spherical Harmonics Modeling¶

Similarly, MC-SH models are instantiated the same

In [12]:

mt_BAS_sh = modeling_framework.MultiCompartmentSphericalHarmonicsModel(
    models=models,
    S0_tissue_responses=S0_responses)

But the fitting procedure is slightly differently implemented as a 1-step optimization because fitting spherical harmonics is convex already.

Implications¶

Including tissue responses allows to correct a signal fraction estimation to a volume fraction estimation. Fitting an MC-model representing tissues with different true S0 responses without including these estimates will result in biased volume fraction estimated.

As a good example, CSF has a much higher S0 response than white matter. This results in vast overestimation of CSF volume fractions in models such as NODDI. We illustrate in the Multi-Tissue NODDI example.