11.9. Computing CRLB for a 2D Diffusion–\(T_2\) Model and comparing two encoding schemes:

In addition to reconstructing spectra, it is often useful to understand how well a given acquisition can estimate the underlying model parameters. The Cramer–Rao lower bound (CRLB) provides a theoretical lower bound on the variance of any unbiased estimator. For a single voxel with signal model \(g(\boldsymbol\theta,\boldsymbol\gamma_p)\) measured at acquisition setting \(\boldsymbol\gamma_p\), additive white Gaussian noise with variance \(\sigma^2\), and parameter vector \(\boldsymbol\theta\), the Fisher information matrix is

\[[\mathbf F]_{ij} \;=\; \frac{1}{\sigma^2} \sum_{p=1}^P \frac{\partial g(\boldsymbol\theta,\boldsymbol\gamma_p)}{\partial \theta_i} \frac{\partial g(\boldsymbol\theta,\boldsymbol\gamma_p)}{\partial \theta_j}, \qquad \mathrm{Cov}(\hat{\boldsymbol\theta}) \succeq \mathbf F^{-1}.\]

Our CRLB Evaluation Tool automates this computation for a single voxel and a user-specified signal model.

In this tutorial, we will walk through a complete CRLB calculation for a 2-compartment diffusion–\(T_2\) model of the form

\[g(b,TE;\boldsymbol\theta) = \sum_{s=1}^S w_s \bigl(1-e^{-b D_s}\bigr) e^{-TE/T_{2,s}},\]

where each component \(s\) has weight \(w_s\), diffusivity \(D_s\), and \(T_{2,s}\) value.

Citation

If you use this tool, please cite:

Kim et al., Improved efficiency for microstructure imaging using high-dimensional MR correlation spectroscopic imaging, 2017. DOI: 10.1109/ACSSC.2017.8335555

Kim et al., Multidimensional correlation spectroscopic imaging of exponential decays: From theoretical principles to in vivo human applications, 2020. DOI: 10.1002/nbm.4244

11.9.1. Goal

This tutorial will guide you through computing the CRLB for the above 2-compartment model at a set of \((b,TE)\) values. Specifically we will:

• Create a funcfile that encodes the model, spectral parameters, acquisition parameters, acquisition grid, and component locations.

• Run the CRLB_save tool to compute the CRLB.

• Interpret the CRLB output in terms of parameter uncertainty.

11.9.2. Relevant links

• Before getting started with this software, the user can download the following file func_expression_diffT2_1.txt from this link here. This file contains informations required for successful CRLB estimation for a diffusion-T2 experiments.

• Full guide to install DRSuite software tools can be found here.

• For matlab based implementation, please download the folders namely, 'utilities', 'solver', 'ini2struct' and the function 'create_phantom.m' provided in the github here (add url).

11.9.3. Get started

1. Create a working folder

   For this tutorial we assume a folder structure like:

   /Documents/CRLB_2D/
   ├── func_expression_diffT2_1.txt
   ├── ini2struct (For MATLAB only)
   └── CRLB_save.m (For MATLAB only)
   

   You can choose any directory, but we will use /Documents/CRLB_2D/ in the examples below.

2. Create the function expression file

   Open a text editor and create a file named func_expression_diffT2.txt inside CRLB_2D with the following content:

   % Signal model: g(b,TE) = sum_s w_s (1-exp(-b*D_s)) * exp(-TE/T2_s)
   func         = '(1-exp(-b*D))*(exp(-TE/T2))';
   spect_params = 'D,T2';
   acq_params   = 'b,TE';
   
   % Acquisition grid: rows are [b, TE]
   exp_vals = [ ...
       0,   99;
       0,  120;
       0,  160;
       0,  200;
       0,  250;
       0,  300;
       0,  400;
       0.2, 99;
       0.2, 120;
       0.2, 160;
       0.2, 200;
       0.2, 250;
       0.2, 300;
       0.2, 400;
       0.5, 99;
       0.5, 120;
       0.5, 160;
       0.5, 200;
       0.5, 250;
       0.5, 300;
       0.5, 400;
       1.0, 99;
       1.0, 120;
       1.0, 160;
       1.0, 200;
       1.0, 250;
       1.0, 300;
       1.0, 400;
       1.5, 99;
       1.5, 120;
       1.5, 160;
       1.5, 200;
       1.5, 250;
       1.5, 300;
       1.5, 400;
       2.5, 99;
       2.5, 120;
       2.5, 160;
       2.5, 200;
       2.5, 250;
       2.5, 300;
       2.5, 400;
       5.0, 99;
       5.0, 120;
       5.0, 160;
       5.0, 200;
       5.0, 250;
       5.0, 300;
       5.0, 400];
   
   % Component prior locations: [w_s, D_s, T2_s]
   % Here T2 is in ms, TE is also in ms, and b*D must be dimensionless.
   components = [1, 0.5,  70;   % component 1
                 1, 1.3, 100];  % component 2
   
   % Assumed noise standard deviation (in signal units)
   noise_std = 0.05;
   

   Notes:

   • The names of the variables (func, spect_params, acq_params, exp_vals, components, and optional noise_std) must match the required format described in Cramér-Rao Bounds.

   • Ensure the units are consistent. In this example, both TE and T2 are in ms. The product b*D is dimensionless.

3. Run the CRLB tool

   Now you can run the CRLB calculation tool from the commandline or in MATLAB using the command below,

   SH

   CRLBsave.sh -funcfile  Documents/CRLB_2D/func_expression_diffT2.txt -outprefix Documents/CRLB_2D/CRLB_diffT2
   

   MATLAB

   CRLB_save('funcfile',  'Documents/CRLB_1D/func_expression_diffT2.txt','outprefix', 'Documents/CRLB_1D/CRLB_diffT2');
   

4. Interpreting the Output

   The output file CRLB_diffT2.txt is a tab-separated text file with one row per component and three columns:

   CRLB(w_1)    CRLB(D_1)    CRLB(T2_1)
   CRLB(w_2)    CRLB(D_2)    CRLB(T2_2)
   

   For each component \(s\):

   • Column 1 is \(\mathrm{CRLB}(w_s)\) (variance lower bound on \(w_s\)),

   • Column 2 is \(\mathrm{CRLB}(D_s)\),

   • Column 3 is \(\mathrm{CRLB}(T_{2,s})\).

5. To obtain the corresponding minimum achievable standard deviation, take square roots:

   CRB = readmatrix('Documents/CRLB_1D/CRLB_diffT2.txt');
   std_lower_bounds = sqrt(CRB);
   

   If these standard deviations are large compared to the true parameter values, it indicates that the current acquisition and SNR do not provide strong information about those parameters (i.e., the problem is poorly conditioned or the SNR is too low).

   Note

   • The CRLB values scale as \(\text{noise_std}^2\). If you halve noise_std, the CRLB entries decrease by a factor of 4.

   • You can explore the effectof different acquisition designs by changing exp_vals (e.g., adding more b-values or shorter TEs) and recomputing the CRLB to see how the information content changes.

6. Now let’s run the same algorithm with the other input func_expression_diffT2_1.txt. This file is exactly the same as the previous file except the contrast encodings are downsampled,

   SH

   estimate_crlb.sh -funcfile  Documents/CRLB_2D/func_expression_diffT2_1.txt -outprefix Documents/CRLB_2D/CRLB_diffT2_1
   

   MATLAB

   estimate_crlb('funcfile',  'Documents/CRLB_1D/func_expression_diffT2_1.txt','outprefix', 'Documents/CRLB_1D/CRLB_diffT2_1');
   

7. Now if you check the two results you will find out that the CRLB values are lower for the first contrast encoding set. So you may choose to use those contrast encodings in your experiment to estimate the diffusion-T2 spectrum.

In summary, this tutorial demonstrates how to use CRLB_save to quantify the theoretical precision limits for a multi-compartment diffusion–\(T_2\) model, given a specific acquisition and noise level.