Bug when running multi-tissue MultiCompartmentSphericalMeanModel

[villalonreina]

@rutgerfick We are trying to run a multi-tissue, multi-compartment SMT model, with the classic stick and zeppelin bundled together, and tortuosity and equality constraints included. We have precalculated a tissue response separately for the whole population (S0_gm) which is fed into the framework. See below a snippet of the code and below that the error message:

# First we call the Stick and the Zeppelin:
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

# We must create a BundleModel with Stick and zeppelin
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')

# Combine stick and zeppelin in a multi-compartment, multi-tissue model:
mcmdi_mod = modeling_framework.MultiCompartmentSphericalMeanModel(
        models=[bundle], S0_tissue_responses=[S0_gm])
mcmdi_mod.parameter_names

# Fitting MC-MDI with and without multi-tissue
print('Grey matter multi-tissue MC-MDI fitting starting')
t = time()
mcmdi_fit = mcmdi_mod.fit(scheme, masked_data, mask=datamask)

This is the error message. It fails in the secondary multi-tissue optimization:

Setup brute2fine optimizer in 0.01388406753540039 seconds
Fitting of 170694 voxels complete in 2132.7913370132446 seconds.
Average of 0.01249482311629726 seconds per voxel.
Starting secondary multi-tissue optimization.
Traceback (most recent call last):
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 2712, in safe_execfile
    self.compile if shell_futures else None)
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/IPython/utils/py3compat.py", line 168, in execfile
    exec(compiler(f.read(), fname, 'exec'), glob, loc)
  File "/Users/jvillalo/Documents/IGC/Python_tests/Python_codes/mcmdi_multitissue_S0_GM.py", line 119, in <module>
    mcmdi_fit = mcmdi_mod.fit(scheme, masked_data, mask=datamask)
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/dmipy-0.1.dev0-py3.7.egg/dmipy/core/modeling_framework.py", line 1622, in fit
    mt_fractions[idx] = fit_func(voxel_S, parameters)
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/dmipy-0.1.dev0-py3.7.egg/dmipy/optimizers/multi_tissue_convex_optimizer.py", line 58, in __call__
    maxfun=2000)
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/scipy/optimize/cobyla.py", line 167, in fmin_cobyla
    **opts)
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/scipy/optimize/cobyla.py", line 233, in _minimize_cobyla
    f = c['fun'](x0, *c['args'])
  File "/Users/jvillalo/anaconda3/lib/python3.7/site-packages/dmipy-0.1.dev0-py3.7.egg/dmipy/optimizers/multi_tissue_convex_optimizer.py", line 64, in cobyla_positivity_constraint
    return volume_fractions - 0.001
TypeError: unsupported operand type(s) for -: 'list' and 'float'

[rutgerfick]

Ah I think I know what's happening. I'll fix it shortly.

[villalonreina]

@rutgerfick
I changed my code by adding two tissue responses in the modeling_framework like this:

# First we call the Stick and the Zeppelin:
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

# We must create a BundleModel with Stick and zeppelin
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')

# Combine stick and zeppelin in a multi-compartment, multi-tissue model:
mcmdi_mod = modeling_framework.MultiCompartmentSphericalMeanModel(
        models=[bundle], S0_tissue_responses=[S0_gm, S0_gm])
mcmdi_mod.parameter_names

It is now giving this error:

Traceback (most recent call last):
  File "/ifs/loni/faculty/thompson/four_d/jvillalo/scripts2/Python_scripts/mcmdi_multitissue_S0_GM.py", line 116, in <module>
    models=[bundle], S0_tissue_responses=[S0_gm, S0_gm])
  File "/ifs/loni/faculty/thompson/four_d/jvillalo/Programs/dmipy/dmipy/core/modeling_framework.py", line 1380, in __init__
    msg.format(len(S0_tissue_responses), self.N_models))
ValueError: Number of S0_tissue responses 2 must be same as number of input models 1.

[rutgerfick]

this is happening because models=[bundle] is of length 1 and you're giving two responses S0_tissue_responses=[S0_gm, S0_gm].

if you want to give different S0 responses to the things inside the Bundle, then you need to bring everything to the top level. But, you will not be able to use MIX optimization in this case when you use tortuosity, but you can use brute2fine.

[villalonreina]

Yes, we want to assign different S0 responses to each compartment. To bring everything to the top level basically means not to bundle the compartments at all?

[rutgerfick]

yes indeed, so for your example it would be something like

 First we call the Stick and the Zeppelin:
stick = cylinder_models.C1Stick()
zeppelin = gaussian_models.G2Zeppelin()

# Combine stick and zeppelin in a multi-compartment, multi-tissue model:
mcmdi_mod = modeling_framework.MultiCompartmentSphericalMeanModel(
        models=[stick, zeppelin], S0_tissue_responses=[S0_gm, S0_gm])

# have to change these parameter names but you get it.
mcmdi_mod.set_tortuous_parameter('G2Zeppelin_1_lambda_perp',
                              'C1Stick_1_lambda_par', 'partial_volume_0', 'partial_volume_1')
mcmdi_mod.set_equal_parameter('G2Zeppelin_1_lambda_par', 'C1Stick_1_lambda_par')

mcmdi_mod.parameter_names

[matteofrigo]

Hi @villalonreina , we just added a new feature that may help you in what you're doing (update your dmipy installation while you read this). It concerns the use of the tortuosity constraint with multi-tissue, which needs to take into account the difference between volume fractions and signal fractions. In the following lines of code you'll find the easy way to encapsulate all the convoluted models in the top level and use the corrected tortuosity constraint.

First, define the multi-compartment model with the S0s associated to each compartment.

from dmipy.signal_models import cylinder_models, gaussian_models 
from dmipy.distributions.distribute_models import SD1WatsonDistributed
from dmipy.core.modeling_framework import MultiCompartmentModel
stick = cylinder_models.C1Stick() 
wstick = SD1WatsonDistributed([stick])
zeppelin = gaussian_models.G2Zeppelin()
wzepp = SD1WatsonDistributed([zeppelin])
ball = gaussian_models.G1Ball() 
s0 = [3000, 4000, 10000]
model = MultiCompartmentModel([wstick, wzepp, ball], S0_tissue_responses=s0)

Now you can add the constraints that will bundle the stick and the zeppelin together and the ones related to tortuosity. Right now the free parameters are

In [3]: model.parameter_names                                                   
Out[3]: 
['SD1WatsonDistributed_1_C1Stick_1_lambda_par',
 'SD1WatsonDistributed_1_SD1Watson_1_mu',
 'SD1WatsonDistributed_1_SD1Watson_1_odi',
 'SD1WatsonDistributed_2_G2Zeppelin_1_lambda_par',
 'SD1WatsonDistributed_2_G2Zeppelin_1_lambda_perp',
 'SD1WatsonDistributed_2_SD1Watson_1_mu',
 'SD1WatsonDistributed_2_SD1Watson_1_odi',
 'G1Ball_1_lambda_iso',
 'partial_volume_0',
 'partial_volume_1',
 'partial_volume_2']

You can see that there are two watson distributions that need to be bundled together. To do that, it is sufficient to set the orientation and dispersion parameters as equal.

model.set_equal_parameter('SD1WatsonDistributed_1_SD1Watson_1_mu',
    'SD1WatsonDistributed_2_SD1Watson_1_mu')
model.set_equal_parameter('SD1WatsonDistributed_1_SD1Watson_1_odi',
    'SD1WatsonDistributed_2_SD1Watson_1_odi')

Now, in order to define the tortuosity constraint, we can set it up in this way.

S0_correction = True
model.set_tortuous_parameter('SD1WatsonDistributed_2_G2Zeppelin_1_lambda_perp',
    'SD1WatsonDistributed_2_G2Zeppelin_1_lambda_par',
    'partial_volume_0', 'partial_volume_1',
    S0_correction=S0_correction)

This should produce the following output.

Employing S0 correction of tortuosity constraint with:
S0_intra: 3000
S0_extra: 4000

The magic that is happening is that the S0 of the intra-axonal and the extra-axonal compartment that you passed when you defined the model are now employed in the calculation of the volume fraction used in the definition of the tortuosity constraint. The math behind it is

lambda_perp = lambda_par * (1 - volume_fraction_intra)

where volume_fraction_intra is defined using the signal fraction and the S0s of the involved compartments.

Now you've still got some free parameters that you may want to fix or link, like the parallel diffusivity of the stick and the zeppelin or the isotropic diffusivity of the ball, but that's up to you and depends on what's the model you want to define.

Hope this helps!

[villalonreina]

Thanks a lot @matteofrigo. Since we are mostly using the spherical mean technique for our grey matter calculations, could this trick also work in the MCMDI framework? Thanks!

[villalonreina]

Hi @rutgerfick @matteofrigo,
We installed the master branch in October 2020, and it still doesn’t have the S0_correction flag for tortuosity in the multi-tissue (MT) models, NODDI or SMT, as shown by Matteo earlier in the thread:

TypeError: set_tortuous_parameter() got an unexpected keyword argument 'S0_correction'

What is the branch we need to install in order to get this feature? Also, I’m confirming that the model as written out by @matteofrigo, runs but the outputs are nonsensical without the extra 'S0_correction' flag.
Since we are mostly using the spherical mean technique for our grey matter calculations, could this trick also work in the MCMDI framework and is it already implemented in the same branch?

Finally, given that our SMT-multitissue and SMT non-multitissue ICVF outputs look very similar without this tortuosity correction, is it safe to assume this is critical to claim that we are using a multitissue approach.

[rutgerfick]

hey @matteofrigo, @villalonreina is right that this feature is not yet merged.

We have a whole PR prepared but you still needed to address the final comments for the code review : #84

can we wrap this one up?