Calculating a Multi-Tissue ICVF from NODDI

[TMNir]

@matteofrigo @rutgerfick In your example for NODDI Watson (https://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_noddi_watson.ipynb) at the very end you explain that to get the ICVF you multiply the partial_volume_1 * SD1WatsonDistributed_1_partial_volume_0

I am following your multi-tissue (MT) NODDI example (https://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_multi_tissue_noddi.ipynb)
Where the last thing you show are the partial_volume_1 and MT_partial_volume_1 outputs.

To get the ICVF for the non multi-tissue model, as before, I can multiply:
partial_volume_1 * SD1WatsonDistributed_1_partial_volume_0
To get the ICVF for the multi-tissue model, I assume I multiply:
MT_partial_volume_1 * SD1WatsonDistributed_1_partial_volume_0
as there is no MT_SD1WatsonDistributed_1_partial_volume_0 output in NODDI_fit.fitted_multi_tissue_fractions_normalized. Is that correct? I just want to confirm that there shouldn’t be an MT normalized volume fraction of the Stick

Thanks so much for your help!!
Best
Talia

[rutgerfick]

Hi Talia!

Yes that is completely correct :-)
the SD1WatsonDistributed_1_partial_volume_0 is exactly just the normalized fraction of the stick inside the Watson Distributed model.
So to get the ICVF of the MT-model, you are correct that you need to do MT_partial_volume_1 * SD1WatsonDistributed_1_partial_volume_0.

Very happy to help you on your way!
Rutger

[TMNir]

Amazing! Thanks so much for the response! So the reason I asked actually is that we were comparing the two versions of multi-compartment SMT you have in the examples:
1 is bundled:
https://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_smt_noddi.ipynb
https://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_multi_compartment_spherical_mean_technique.ipynb

1 is not:
https://nbviewer.jupyter.org/github/AthenaEPI/dmipy/blob/master/examples/example_spherical_mean_models.ipynb

The non-bundled version outputs the ICVF while for the bundled version we have to carry out the multiplication we noted above.

We were outputting both the multi tissue and standard outputs for each as we are mainly interested in the MT outputs:
mcmdi_fit.fitted_parameters
mcmdi_fit.fitted_multi_tissue_fractions_normalized

For the standard outputs--> unbundled ICVF (partial_volume_0) and the bundled ICVF (BundleModel_1_partial_volume_0 * partial_volume_1) are essentially identical

For the MT outputs-->unbundled MT-ICVF (MT_partial_volume_0) and the bundled MT-ICVF (BundleModel_1_partial_volume_0 * MT_partial_volume_1) are quite different, with the unbundled version looking very pixelated.

Any ideas why?

Thanks again!

[matteofrigo]

Hi @TMNir !

Using BundleModel forces the orientation of the compartment to be equal, but the diffusivities and the other parameters remain untouched. Please @rutgerfick correct me if I'm wrong on this.

The two models that you listed are not equivalent as, despite being both defined as the combination of a stick and a zeppelin, they have different sets of constraints. In Dmipy you can setup three versions of that model.

1. Without constraints
2. Forcing the constraints parameter-by-parameter
3. Using BundleModel

Here's an example of what I mean. I will not use the multi-tissue features, as they are not related to the problem. You can obtain the same results specifying the S0 of each compartment as usual. Also, I will consider spherical-mean models, hence the orientation is not an issue. If you use standard multi-compartment models, you should set the orientation of the stick and the zeppelin equal.

# First load the example data
from dmipy.data import saved_data
scheme_hcp, data_hcp = saved_data.wu_minn_hcp_coronal_slice()

Model 1: Not bundled

from dmipy.signal_models import cylinder_models, gaussian_models
from dmipy.core.modeling_framework import MultiCompartmentSphericalMeanModel
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()
mcmdi_not_bundled = MultiCompartmentSphericalMeanModel(
    models=[stick, zeppelin])
print('Free parameters:', mcmdi_not_bundled.parameter_names)

Free parameters: ['C1Stick_1_lambda_par', 'G2Zeppelin_1_lambda_par', 'G2Zeppelin_1_lambda_perp', 'partial_volume_0', 'partial_volume_1']

Notice that here we have two independent parallel diffusivities and one perpendicular diffusivity to fit.

Model 2: force parameters to be equal at top level

stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()
mcmdi_forced_bundled = MultiCompartmentSphericalMeanModel(
    models=[stick, zeppelin])
mcmdi_forced_bundled.set_equal_parameter('C1Stick_1_lambda_par', 'G2Zeppelin_1_lambda_par')
mcmdi_forced_bundled.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
    'C1Stick_1_lambda_par', 'partial_volume_0', 'partial_volume_1')
print('Free parameters:', mcmdi_forced_bundled.parameter_names)

Free parameters: ['C1Stick_1_lambda_par', 'partial_volume_0', 'partial_volume_1']

Here we have only one parallel diffusivity to fit. The other one is set equal to the first one, and the perpendicular diffusivity is

Model 3: automatic bundle

from dmipy.distributions.distribute_models import BundleModel
bundle = BundleModel([stick, zeppelin])
bundle.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
    'C1Stick_1_lambda_par','partial_volume_0')
bundle.set_equal_parameter('G2Zeppelin_1_lambda_par', 'C1Stick_1_lambda_par')

mcmdi_auto_bundled = MultiCompartmentSphericalMeanModel(
    models=[bundle])

print('Free parameters:', mcmdi_forced_bundled.parameter_names)

Free parameters: ['BundleModel_1_G2Zeppelin_1_lambda_par', 'BundleModel_1_partial_volume_0']

Notice that these free parameters correspond to the ones of model number 2 (where we forced parameters to be equal at the top level).

Experiment

We fit the three models, expecting to see this:

1. The first model is a degenerate version of the other two.
2. The second and third model are equivalent

fit_not_bundled = mcmdi_not_bundled.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)
fit_forced_bundled = mcmdi_forced_bundled.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)
fit_auto_bundled = mcmdi_auto_bundled.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)

Results

Let's plot the intra-cellular volume fraction obtained with each model.

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=[15, 15])

ax = axes[0]
ax.set_title('ICvf NON BUNDLED')
data = fit_not_bundled.fitted_parameters['partial_volume_0'].squeeze()
ax.imshow(data.T, origin=True)

ax = axes[1]
ax.set_title('ICvf FORCED BUNDLED')
data = fit_forced_bundled.fitted_parameters['partial_volume_0'].squeeze()
ax.imshow(data.T, origin=True)

ax = axes[2]
ax.set_title('ICvf AUTO BUNDLED')
data = fit_auto_bundled.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze()
ax.imshow(data.T, origin=True)

plt.show()

The mean squared difference between the results obtained with the second and the third models is:

import numpy as np
np.mean(np.abs(fit_forced_bundled.fitted_parameters['partial_volume_0'].squeeze() -
fit_auto_bundled.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze()) ** 2)

2.024555072199328e-05

Conclusion

I hope this convinced you that the problem was the set of constraints applied to each model. the smooth solution is obtained when the parallel diffusivity is equal for the stick and the zeppelin compartment, and the perpendicular diffusivity is set with the tortuosity constraint.

If you employ the standard multi-compartment formulation (i.e., without spherical mean), remember to set the orientation of the stick and the zeppelin to be equal.

In this context, setting the S0 of each compartment can be done as usual and will not change the message of this reply, as the problem was in the constraints, not in the S0 responses.

I hope this helps!

Cheers

[TMNir]

Thanks! Yes, we did in fact compare something analogous to your Model 2 and Model 3 and for the non MT outputs, as you also showed, the results were identical. Again, for the non MT outputs, the results are the same, which is a nice sanity check, but the difference is only in the MT outputs.

These are the analogous two models:

Model 2:

ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

mcmdi_model = MultiCompartmentSphericalMeanModel(models=[stick, zeppelin, ball], 
                              S0_tissue_responses=[S0_gm, S0_gm, S0_csf])

mcmdi_model.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
                                   'C1Stick_1_lambda_par',
                                    'partial_volume_0',
                                   'partial_volume_1')
mcmdi_model.set_equal_parameter('C1Stick_1_lambda_par',
                                 'G2Zeppelin_1_lambda_par')
mcmdi_model.set_fixed_parameter('C1Stick_1_lambda_par', 1.1e-9)
mcmdi_model.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9) 

Model 3:

ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

bundle = BundleModel(models=[stick, zeppelin])
bundle.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
               'C1Stick_1_lambda_par','partial_volume_0')

mcmdi_Bmodel = MultiCompartmentSphericalMeanModel(models=[ball,bundle],
                            S0_tissue_responses=[S0_csf,S0_gm])

mcmdi_Bmodel.set_equal_parameter('BundleModel_1_C1Stick_1_lambda_par',
                        'BundleModel_1_G2Zeppelin_1_lambda_par')
mcmdi_Bmodel.set_fixed_parameter('BundleModel_1_C1Stick_1_lambda_par', 1.1e-9)
mcmdi_Bmodel.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9)

[matteofrigo]

Let's define the models as follows (i.e., with your code and two fake S0s):

S0_gm, S0_csf = [3000, 5000]

# model 2
ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()
mcmdi_model = MultiCompartmentSphericalMeanModel(models=[stick, zeppelin, ball], 
                              S0_tissue_responses=[S0_gm, S0_gm, S0_csf])
mcmdi_model.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
                                   'C1Stick_1_lambda_par',
                                    'partial_volume_0',
                                   'partial_volume_1')
mcmdi_model.set_equal_parameter('C1Stick_1_lambda_par',
                                 'G2Zeppelin_1_lambda_par')
mcmdi_model.set_fixed_parameter('C1Stick_1_lambda_par', 1.1e-9)
mcmdi_model.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9) 

# model 3
ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()
bundle = BundleModel(models=[stick, zeppelin])
bundle.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
               'C1Stick_1_lambda_par','partial_volume_0')
mcmdi_Bmodel = MultiCompartmentSphericalMeanModel(models=[ball,bundle],
                            S0_tissue_responses=[S0_csf,S0_gm])
mcmdi_Bmodel.set_equal_parameter('BundleModel_1_C1Stick_1_lambda_par',
                        'BundleModel_1_G2Zeppelin_1_lambda_par')
mcmdi_Bmodel.set_fixed_parameter('BundleModel_1_C1Stick_1_lambda_par', 1.1e-9)
mcmdi_Bmodel.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9)

The two models with have these free parameters:

mcmdi_model.parameter_names:
['partial_volume_0', 'partial_volume_1', 'partial_volume_2']

mcmdi_Bmodel.parameter_names:
['partial_volume_0', 'partial_volume_1', 'partial_volume_2']

Running this optimization

fit1 = mcmdi_model.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)
fit2 = mcmdi_Bmodel.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)

I get a mean squared error of 9.556214660115517e-05 . Visually, they look the same too.

Notice that the ICvf is BundleModel_1_partial_volume_0 * partial_volume_1 for model 3, while for model 2 it is sufficient to look at partial_volume_0. The reason is that the bundled stick and zeppelin have a partial fraction within the bundle that needs to be rescaled by the fraction of the signal explained by the whole bundle.

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=[15, 15])

ax = axes[0]
ax.set_title('Model 2')
data = fit1.fitted_parameters['partial_volume_0'].squeeze()
ax.imshow(data.T, origin=True)

ax = axes[1]
ax.set_title('Model 3')
data = fit2.fitted_parameters['partial_volume_1'].squeeze() * fit2.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze()
ax.imshow(data.T, origin=True)
plt.show()

[TMNir]

Yes. I get the same as well. But the images I pasted are from when I compare:
fit1.fitted_multi_tissue_fractions_normalized['partial_volume_0']
and
fit2.fitted_multi_tissue_fractions_normalized['partial_volume_1'] * it2.fitted_parameters['BundleModel_1_partial_volume_0']

Aren't those the MT ICVF outputs?

[matteofrigo]

Yes, they are. I had misunderstood the question yet another time. My bad, sorry.

The patches are due to the same effect that we discussed in another issue , namely to the amplification of the differences between IC/EC and CSF due to the high difference between the S0 of the GM and the one of CSF. It's less smooth because multi-tissue has more contrast.

I also did some experiments and the short conclusion is: use the Bundled model.

Here's what I looked at:

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import matplotlib.pyplot as plt

# First load the example data
from dmipy.data import saved_data
scheme_hcp, data_hcp = saved_data.wu_minn_hcp_coronal_slice()
means0 = np.mean(data_hcp[..., scheme_hcp.b0_mask], axis=3)

from dmipy.signal_models import cylinder_models, gaussian_models
from dmipy.core.modeling_framework import MultiCompartmentSphericalMeanModel
from dmipy.distributions.distribute_models import BundleModel

S0_gm, S0_csf = [3000, 5000]

Model 1

ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

mcmdi_model = MultiCompartmentSphericalMeanModel(models=[stick, zeppelin, ball], 
                              S0_tissue_responses=[S0_gm, S0_gm, S0_csf])

mcmdi_model.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
                                   'C1Stick_1_lambda_par',
                                    'partial_volume_0',
                                   'partial_volume_1')
mcmdi_model.set_equal_parameter('C1Stick_1_lambda_par',
                                 'G2Zeppelin_1_lambda_par')
mcmdi_model.set_fixed_parameter('C1Stick_1_lambda_par', 1.1e-9)
mcmdi_model.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9) 

mcmdi_model.parameter_names

['partial_volume_0', 'partial_volume_1', 'partial_volume_2']

Model 2

ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

bundle = BundleModel(models=[stick, zeppelin])
bundle.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
               'C1Stick_1_lambda_par','partial_volume_0')

mcmdi_Bmodel = MultiCompartmentSphericalMeanModel(models=[ball,bundle],
                            S0_tissue_responses=[S0_csf,S0_gm])

mcmdi_Bmodel.set_equal_parameter('BundleModel_1_C1Stick_1_lambda_par',
                        'BundleModel_1_G2Zeppelin_1_lambda_par')
mcmdi_Bmodel.set_fixed_parameter('BundleModel_1_C1Stick_1_lambda_par', 1.1e-9)
mcmdi_Bmodel.set_fixed_parameter('G1Ball_1_lambda_iso', 3e-9)

mcmdi_Bmodel.parameter_names

['BundleModel_1_partial_volume_0', 'partial_volume_0', 'partial_volume_1']

Fit

fit1 = mcmdi_model.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)

Using parallel processing with 4 workers.
Setup brute2fine optimizer in 1.6738982200622559 seconds
Fitting of 8181 voxels complete in 46.548293352127075 seconds.
Average of 0.005689804834632328 seconds per voxel.
Starting secondary multi-tissue optimization.
Multi-tissue fitting of 8181 voxels complete in 35.5544228553772 seconds.

fit2 = mcmdi_Bmodel.fit(scheme_hcp, data_hcp, mask=data_hcp[..., 0]>0)

Using parallel processing with 4 workers.
Setup brute2fine optimizer in 0.19896984100341797 seconds
Fitting of 8181 voxels complete in 54.375812292099 seconds.
Average of 0.006646597273206087 seconds per voxel.
Starting secondary multi-tissue optimization.
Multi-tissue fitting of 8181 voxels complete in 4.872658014297485 seconds.

fit1.fitted_multi_tissue_fractions.keys()

dict_keys(['partial_volume_0', 'partial_volume_1', 'partial_volume_2'])

fit2.fitted_multi_tissue_fractions.keys()

dict_keys(['partial_volume_0', 'partial_volume_1'])

Results

Signal Fractions

data1sf = fit1.fitted_parameters['partial_volume_0'].squeeze()
data2sf = fit2.fitted_parameters['partial_volume_1'].squeeze() * fit2.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze()

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=[15, 15])

ax = axes[0]
ax.set_title('m1 mt')
ax.imshow(data1sf.T, origin=True, interpolation='nearest')

ax = axes[1]
ax.set_title('m2 mt')
ax.imshow(data2sf.T, origin=True, interpolation='nearest')

Root mean squared difference between ICsf computed with module 1 and module 2

np.sqrt(np.mean(np.abs(data1sf - data2sf) ** 2))

0.009775589322447787

Rescaling signal fractions to get volume fractions

Volume fractions can also be computed from signal fractions as follows:

where $\bar{S}_0 is$ the mean S0 image, $f_i$ is the volume fraction, $\phi_i$ is the signal fraction and $S_0^i$ is the S0 associated to compartment i.

(sidenote: this may become the default behaviour in dmipy in the future. we'll keep you posted)

data1 = fit1.fitted_parameters['partial_volume_0'].squeeze() * means0.squeeze() / S0_gm

data2 = fit2.fitted_parameters['partial_volume_1'].squeeze() * fit2.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze() * means0.squeeze() / S0_gm

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=[15, 15])

ax = axes[0]
ax.set_title('m1 mt')
ax.imshow(data1.T, origin=True, interpolation='nearest')

ax = axes[1]
ax.set_title('m2 mt')
ax.imshow(data2.T, origin=True, interpolation='nearest')

Using volume fractions fitted by dmipy

data1mt = fit1.fitted_multi_tissue_fractions['partial_volume_0'].squeeze()

data2mt = fit2.fitted_multi_tissue_fractions['partial_volume_1'].squeeze() * fit2.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze()

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=[15, 15])

ax = axes[0]
ax.set_title('m1 mt')
ax.imshow(data1mt.T, origin=True, interpolation='nearest')

ax = axes[1]
ax.set_title('m2 mt')
ax.imshow(data2mt.T, origin=True, interpolation='nearest')

plt.show()

Model 1

Root mean squared difference between ICvf computed rescaling the single-tissue fractions and the multi-tissue one computed by dmipy.

np.sqrt(np.mean(np.abs(data1 - data1mt) ** 2))

0.018687889476120545

Model 2

Root mean squared difference between ICvf computed rescaling the single-tissue fractions and the multi-tissue one computed by dmipy.

np.sqrt(np.mean(np.abs(data2 - data2mt) ** 2))

0.0016351014683166412

Difference between the two using the rescaled fractions

np.sqrt(np.mean(np.abs(data1 - data2) ** 2))

0.014210827418464472

Difference between the two using the dmipy multi tissue fractions

np.sqrt(np.mean(np.abs(data1mt - data2mt) ** 2))

0.02196077634801288

Observations

• The two models yield different signal fractions, with a gap of 9e-3.
• The two models yield different volume fractions.
  • Volume fractions computed by rescaling the signal fractions have a RMS difference of 1.4e-2 between the two models.
  • Volume fractions computed natively by dmipy have a RMS difference of 2.1e-2 between the two models.
• Within a single model, the volume fractions computed natively by dmipy are different with respect to the volume fractions obtained by rescaling the signal fractions.
  • Model 1 - RMS difference: 1.8e-2
  • Model 2 - RMS difference: 1.6e-3

Conclusions

The two models differ are mathematically equivalent, but have some major numerical differences. The technical aspect that distinguishes the two is related to the amout of computations that are necessary to compute the response function of the GM compartments. Model 1 needs more operations than module 2. Considering that these operations are approximated, model 1 is intrinsically less stable (or more approximated) than model 2.
This is partially reflected by the results observed above, where model 2 is shown to yield closer volume fractions when they are computed either by rescale or natively by dmipy. For this reason, model 2 should be preferred.

Also, rescaling the signal fractions is definitely a cheap and trustworth way to get the volume fractions.

[villalonreina]

Hi @matteofrigo @rutgerfick. Now that we know that we should use a bundled model, the question now is how to cascade volume fractions from one bundled model to another bundled model. For instance, if we would like to cascade the ICVF from a bundled model to a second bundled model.
In your example:

data2mt = fit2.fitted_multi_tissue_fractions['partial_volume_1'].squeeze() * fit2.fitted_parameters['BundleModel_1_partial_volume_0'].squeeze()

When cascading the ICVF into another bundled model, we would definitely cascade partial_volume_1, but we are not sure what to do about BundleModel_1_partial_volume_0. Should we cascade both? We do not want to disturb the ECVF (Zeppelin) fraction.
Thanks!

[rutgerfick]

Basically, if you don't cascade a parameter (provide an initial guess), then that parameter will be brute force optimized between the min and max parameter bounds.

So if you want to optimize the second model, while keeping an idea of what the previous model found, you should cascade also the bundle parameter.

Otherwise, it will re-optimize the bundle fraction from scratch