Panagiotaki et al. (2014) proposed a multi-compartment model called VERDICT to characterize the composition of tumorous tissues. VERDICT models the diffusion in tumor cells, the extra-cellular space and surrounding bloodvessels as a restricted Sphere, an isotropic Gaussian Ball and a Stick compartment, respectively. VERDICT's design is as follows: \begin{equation} E_{\textrm{VERDICT}}= \underbrace{f_{\textrm{Tumor}}\overbrace{E(D|\lambda_{\textrm{intra}})}^{\textrm{Sphere}}}_{\textrm{Tumor Cells}} + \underbrace{f_{\textrm{extra}}\overbrace{E(\cdot|\lambda_{iso})}^{\textrm{Ball}}}_{\textrm{Hindered Extra-Cellular}}+\underbrace{f_{blood}\overbrace{E(\lambda_\parallel, \boldsymbol{\mu})}^{\textrm{Stick}}}_{\textrm{Vascular}} \end{equation} where $D$ is the sphere's diameter. VERDICT uses the Gaussian Phase approximation to model the sphere (Balinov et al. 1993), which accounts to changes in gradient pulse duration $\delta$ and separation $\Delta$. Furthermore, some particular parameter constraints are imposed:
We can define the VERDICT model in the following lines of code:
We first load the necessary modules to import the Sphere, Ball, and Stick models.
from dmipy.signal_models import sphere_models, cylinder_models, gaussian_models
/home/rutger/anaconda2/lib/python2.7/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`. from ._conv import register_converters as _register_converters
The Sphere model must be initiated with VERDICT's setting for intra-spherical diffusivity, while the Ball and Stick model are regularly imported.
sphere = sphere_models.S4SphereGaussianPhaseApproximation(diffusion_constant=0.9e-9)
ball = gaussian_models.G1Ball()
stick = cylinder_models.C1Stick()
We then assemble the three models into a multi-compartment model:
from dmipy.core.modeling_framework import MultiCompartmentModel
verdict_mod = MultiCompartmentModel(models=[sphere, ball, stick])
verdict_mod.parameter_names
['G1Ball_1_lambda_iso', 'S4SphereGaussianPhaseApproximation_1_diameter', 'C1Stick_1_mu', 'C1Stick_1_lambda_par', 'partial_volume_0', 'partial_volume_1', 'partial_volume_2']
Visualize the model:
from IPython.display import Image
verdict_mod.visualize_model_setup(view=False, cleanup=False, with_parameters=True)
Image('Model Setup.png')
We then fix the Ball's diffusivity, and adjust the Stick's optimization range for $\lambda_\parallel$ for the optimization.
verdict_mod.set_fixed_parameter('G1Ball_1_lambda_iso', 0.9e-9)
verdict_mod.set_parameter_optimization_bounds(
'C1Stick_1_lambda_par', [3.05e-9, 10e-9])
We are now ready to fit the data.
As illustration data we use the voxel of VERDICT data that is freely available at the UCL website at http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.VERDICTcol. We prepared the script to easily import the acquisition scheme and example data from there:
from dmipy.data import saved_data
scheme, data = saved_data.panagiotaki_verdict()
scheme.print_acquisition_info
Acquisition scheme summary total number of measurements: 220 number of b0 measurements: 49 number of DWI shells: 41 shell_index |# of DWIs |bvalue [s/mm^2] |gradient strength [mT/m] |delta [ms] |Delta[ms] |TE[ms] 0 |19 |20 |0 |3.0 |10.0 |17.0 1 |3 |148 |160 |3.0 |10.0 |17.0 2 |3 |231 |200 |3.0 |10.0 |17.0 3 |3 |333 |240 |3.0 |10.0 |17.0 4 |3 |454 |280 |3.0 |10.0 |17.0 5 |3 |593 |320 |3.0 |10.0 |17.0 6 |3 |751 |360 |3.0 |10.0 |17.0 7 |3 |927 |400 |3.0 |10.0 |17.0 8 |12 |4 |0 |3.0 |20.0 |27.0 9 |3 |78 |80 |3.0 |20.0 |27.0 10 |3 |176 |120 |3.0 |20.0 |27.0 11 |3 |313 |160 |3.0 |20.0 |27.0 12 |3 |489 |200 |3.0 |20.0 |27.0 13 |3 |704 |240 |3.0 |20.0 |27.0 14 |3 |959 |280 |3.0 |20.0 |27.0 15 |3 |1253 |320 |3.0 |20.0 |27.0 16 |3 |1585 |360 |3.0 |20.0 |27.0 17 |9 |0 |0 |3.0 |30.0 |37.0 18 |3 |119 |80 |3.0 |30.0 |37.0 19 |3 |268 |120 |3.0 |30.0 |37.0 20 |3 |478 |160 |3.0 |30.0 |37.0 21 |3 |747 |200 |3.0 |30.0 |37.0 22 |3 |1075 |240 |3.0 |30.0 |37.0 23 |3 |1464 |280 |3.0 |30.0 |37.0 24 |3 |1912 |320 |3.0 |30.0 |37.0 25 |3 |2420 |360 |3.0 |30.0 |37.0 26 |3 |2988 |400 |3.0 |30.0 |37.0 27 |13 |9 |0 |3.0 |40.0 |47.0 28 |3 |160 |80 |3.0 |40.0 |47.0 29 |3 |361 |120 |3.0 |40.0 |47.0 30 |3 |643 |160 |3.0 |40.0 |47.0 31 |3 |1004 |200 |3.0 |40.0 |47.0 32 |3 |1446 |240 |3.0 |40.0 |47.0 33 |3 |1969 |280 |3.0 |40.0 |47.0 34 |3 |2572 |320 |3.0 |40.0 |47.0 35 |3 |3255 |360 |3.0 |40.0 |47.0 36 |3 |4018 |400 |3.0 |40.0 |47.0 37 |6 |0 |0 |3.5 |8.0 |15.0 38 |42 |958 |400 |3.5 |8.0 |15.0 39 |2 |0 |0 |10.0 |30.0 |44.0 40 |3 |305 |40 |10.0 |30.0 |44.0 41 |3 |2748 |120 |10.0 |30.0 |44.0 42 |3 |0 |0 |10.0 |40.0 |47.0 43 |3 |419 |40 |10.0 |40.0 |47.0 44 |3 |1679 |80 |10.0 |40.0 |47.0 45 |3 |3778 |120 |10.0 |40.0 |47.0
We can see that the VERDICT acquisition scheme is very particular, having many shells with different diffusion times and often only 3 perpendicular measurements per shell, with the exception of one DTI shell with 42 measurements.
To fit the data we use the MIX algorithm (Farooq et al. 2016), which is efficient for finding the global minimum in models with many compartments. We set parallel processing to False since there is only one voxel.
verdict_fit = verdict_mod.fit(scheme, data, solver='mix', use_parallel_processing=False)
Setup MIX optimizer in 2.90870666504e-05 seconds Fitting of 1 voxels complete in 3.41024303436 seconds. Average of 3.41024303436 seconds per voxel.
To illustrate that Dmipy's VERDICT implementation is correct, we show that we can produce a very similar signal fitting plots as the one shown on the UCL website. On the left we show their original graph, and on the right we show Dmipy's predicted signal together with the measured signal attenuation.
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import matplotlib.image as mpimg
img=mpimg.imread("http://camino.cs.ucl.ac.uk/uploads/Tutorials/VCLSsynth.png")
mask_nonzero = scheme.gradient_strengths>0.
G_nonzero = scheme.gradient_strengths[mask_nonzero]
Delta_nonzero = scheme.Delta[mask_nonzero]
delta_nonzero = scheme.delta[mask_nonzero]
predicted_data = verdict_fit.predict()[0]
predicted_data_nonzero = predicted_data[mask_nonzero]
data_nonzero = data[mask_nonzero]
fig, axs = plt.subplots(1, 2, figsize=[30, 10])
axs = axs.ravel()
axs[0].imshow(img)
axs[0].set_title('UCL Verdict Results', fontsize=30)
axs[0].axis('off')
for delta_, Delta_ in np.unique(np.c_[scheme.shell_delta, scheme.shell_Delta], axis=0):
mask = np.all([Delta_nonzero == Delta_, delta_nonzero == delta_], axis=0)
axs[1].plot(G_nonzero[mask], predicted_data_nonzero[mask])
axs[1].scatter(G_nonzero[mask], data_nonzero[mask], s=3., marker='o')
axs[1].set_title('Dmipy Verdict Results', fontsize=30)
axs[1].set_xlabel('G(T/m)', fontsize=20)
axs[1].set_ylabel('S', fontsize=20);
While we didn't exactly duplicate their graph style, we can see that both the plots and signal fitting are very similar between Dmipy's and UCL's implementation.
Showing the fitting parameters, we can also see that Dmipy's estimate of the tumor cell diameter is 1.54e-5m, i.e. 15.4$\mu m$, which falls exactly within the ranges that Panagiotaki et al. reports.
verdict_fit.fitted_parameters
{'C1Stick_1_lambda_par': array([6.58332016e-09]), 'C1Stick_1_mu': array([[ 2.12873298, -2.53619312]]), 'S0': array([[1., 1., 1., 1., 1., 1.]]), 'S4SphereGaussianPhaseApproximation_1_diameter': array([1.54025803e-05]), 'partial_volume_0': array([0.76130392]), 'partial_volume_1': array([0.09744434]), 'partial_volume_2': array([0.14125174])}