Signal to Concentration#

Hide code cell source
import os
import numpy as np
from matplotlib import pyplot as plt
import csv
import seaborn as sns
import pandas as pd
from plotting_results_nb import plot_bland_altman, bland_altman_statistics
import json
from pathlib import Path

Background#

DCE-MRI data are typically acquired with spoiled gradient echo sequences (SPGR). The measured signal intensities can be converted to concentration values using the spoiled gradient echo equation as described in Schabel et al. (Phys Med Biol 2008):

\[ S(T_1,T_2^*) = M_0 \cdot \frac{\mathrm{sin}(\alpha)(1 - e^{-\mathrm{TR}/T_1}) \cdot e^{-\mathrm{TE}/T_2^*}}{1 - e^{-\mathrm{TR}/T_1} \cdot \mathrm{cos}(\alpha)} \]

with \(S\) the measured signal intensity for spoiled gradient echo sequences, \(T_1\), the longitudinal relaxation time, \(T_2^*\) the transverse effective relaxation time, \(\alpha\) the flip angle, TR, the repitition time, and TE the echo time.

By making several simplifications (e.g. \(T_2^*\) is negligible) the relative signal enhancement (\(\Xi\)) can be expressed as (eq 5 in Schabel et al):

\[ \Xi \approx \frac{(E_1 - E_{1,0})(\mathrm{cos}(\alpha)-1)}{(E_{1,0} - 1)(E_1 \mathrm{cos}(\alpha) -1)} \]

where \(E_{1,0} = e^{-\mathrm{TR}R_{1,0}}\) and \(E_{1} = e^{-\mathrm{TR}R_{1}}\), \(R_1 = 1/T_1\)

In this way R_1 can be solved to obtain a nonlinear approximation as (eq 6 in Schabel et al):

\[ R_1^{nl} \approx - \frac{1}{\mathrm{TR}} \cdot \mathrm{log}[\frac{\Xi(E_{1,0} - 1) + E_{1,0}(1- \mathrm{cos}(\alpha))}{1+\mathrm{cos}(\alpha)(\Xi(E_{1,0} - 1)-1)}] \]

This leads to a nonlinear approximation of the the concentration \(C\) (eq 7 in Schabel et al):

\[ C^{nl} \approx \frac{1}{r_1}(R_1^{nl} - R_{1,0}) \]

with \(r_1\) is the relaxivity of the contrast agent being used.

A linear approximation can also be used (eq (8) in Schabel et al.):

\[ C^l \approx \frac{1}{r_1}(R_{1,0}\Xi) \]

In the current code contributions only the nonlinear approximation has been used.

Test data#

Signal intensity curves randomly selected from one patient with DCE-MRI of the uterus. They were converted to concentration using code from the University of Edinburgh but with various flip angles, baseline T1 values etc rather than the actual values used in the acquisition, to test a wider range of possibilites. For a visualization of the test data, see curves below.

Tolerances: absolute: 0.00001 mM/ relative 0.00001

Source: University of Edinburgh, Mechanistic substudy of UCON https://www.birmingham.ac.uk/research/bctu/trials/womens/ucon/ucon-home.aspx used with permission.
Reference to code: Reavey, J.J., Walker, C., Nicol, M., Murray, A.A., Critchley, H.O.D., Kershaw, L.E., Maybin, J.A., 2021. Markers of human endometrial hypoxia can be detected in vivo and ex vivo during physiological menstruation. Hum. Reprod. 36, 941–950.

Import data#

# Load the meta data
meta = json.load(open("../test/results-meta.json"))
# Loop over each entry and collect the dataframe
df = []
for entry in meta:
    if entry['category'] == 'SI_to_Conc':
        fpath, fname, category, method, author = entry.values()
        df_entry = pd.read_csv(Path(fpath, fname)).assign(author=author)
        df.append(df_entry)
    
# Concat all entries
df = pd.concat(df)
# label data source
df['source']=''
df.loc[df['label'].str.contains('original'),'source']='original'

author_list = df.author.unique()
no_authors = len(author_list)

voxel_list = df.label.unique()
no_voxels = len(voxel_list)

# calculate error between measured and reference R_1 values
df['error_conc'] = df['conc_meas'] - df['conc_ref']

tolerances = { 'conc': {'atol' : 0.00001, 'rtol': 0.00001 }}
dt = 2.5 # in s

Visualization of test data#

To get an impression of the test data included so far, the reference concentration time curves are plotted below.

fig, ax = plt.subplots(1, no_voxels, sharex='col', sharey='row', figsize=(12,3))

for current_voxel in range(no_voxels):
    
    subset_data = df[(df['author'] == author_list[0]) & (df['label'] == voxel_list[current_voxel])]
    no_time_points = len(subset_data.conc_meas)
    time_array = np.arange(0, no_time_points*dt,dt)

    ax[current_voxel].plot(time_array, subset_data.conc_ref, color='black', label='ref')

    if current_voxel == 0:
        ax[current_voxel].set_ylabel('C (mM)', fontsize=16)
    
    ax[current_voxel].set_xlabel('time (s)', fontsize=16)
    ax[current_voxel].set_title(voxel_list[current_voxel], fontsize=16)
    
    # set fontsize
    ax[current_voxel].tick_params(axis='x', labelsize=14)
    ax[current_voxel].tick_params(axis='y', labelsize=14)
 

fig.tight_layout()
fig.subplots_adjust(left=0.15, top=0.95)
plt.show()
_images/be94ef75506c8efc5fe49190378b47f6110405871eede8dc7db389d30bdf8052.png

Plot bland altman figure of error between output concentration values and reference values vs the reference values. All datapoints from all test curves are combined into one figure. Six of the seven implementation give near identical results. One implementation showed small deviations (MJT_UoEdinburgh_UK_num, < \(5e^{-7}\) mM), but please note that the limits of the y-axis are extremely small.

# set-up styling for seaborn plots
sns.set(font_scale=1.5)
sns.set_style("ticks")
fig, ax = plt.subplots(1,1, sharey='none', figsize=(10,5))
sns.set(font_scale=1.5)
sns.set_style("whitegrid")
plot_bland_altman(ax, df, tolerances, 'conc', xlim=(0,5), ylim=(-0.000001,0.000001),label_xaxis='$C_{ref}$ (mM)',label_yaxis='$\Delta $C (mM)')
plt.legend(bbox_to_anchor=(1.05, 1),title='ContributionID', loc='upper left', borderaxespad=0.5, fontsize=14)
fig.tight_layout()
fig.subplots_adjust(left=0.15, top=0.95)

                        
_images/4fc476377036fa9ab3e34d2f1da72348262a3e9ce1e542115aa53056d64f2e09.png

Calculate bias and limits of agreement for all contributions wrt the reference values.

The differences between measured and reference values are really small for the current test data. There are slight differences between contributions, that can be attributed due to minor diferences in implementations.

resultsBA = bland_altman_statistics(data=df,par='error_conc',grouptag='author')

print('Bias results estimated concentration values combined for all voxels')
print(resultsBA)
Bias results estimated concentration values combined for all voxels
                                  bias         stdev     LoA lower  \
author                                                               
LCB_BNI_USA              -6.225602e-14  1.002467e-13 -2.587395e-13   
LEK_UoEdinburgh_UK       -7.354648e-14  9.498901e-14 -2.597249e-13   
MB_UoManchester_UK       -7.425111e-14  9.539893e-14 -2.612330e-13   
MJT_UoEdinburgh_UK       -1.570635e-13  3.697846e-12 -7.404842e-12   
MJT_UoEdinburgh_UK_num   -1.759729e-08  1.265300e-07 -2.655961e-07   
OG_MO_AUMC_NL_ICR_RMH_UK -7.354402e-14  9.498492e-14 -2.597145e-13   
ST_USyd_Aus              -5.996202e-14  9.795241e-14 -2.519487e-13   

                             LoA upper  
author                                  
LCB_BNI_USA               1.342275e-13  
LEK_UoEdinburgh_UK        1.126320e-13  
MB_UoManchester_UK        1.127308e-13  
MJT_UoEdinburgh_UK        7.090715e-12  
MJT_UoEdinburgh_UK_num    2.304015e-07  
OG_MO_AUMC_NL_ICR_RMH_UK  1.126264e-13  
ST_USyd_Aus               1.320247e-13  

Notes#

For some contributions the test data needed to be modified. Please check the test code for modification in input data to be able to run the tests.

References#

Schabel MC and Parker DL “Uncertainty and bias in contrast concentration measurements using spoiled gradient echo pulse sequences” Phys Med Biol (2008), doi:10.1088/0031-9155/53/9/010