Two-compartment exchange model (2CXM)#
Show code cell source
# import statements
import os
import numpy as np
from matplotlib import pyplot as plt
import csv
import pandas as pd
import seaborn as sns
from plotting_results_nb import plot_bland_altman, bland_altman_statistics, make_catplot
import json
from pathlib import Path
Test data#
Data summary: simulated two-compartment exchange model data
Source: Concentration-time data (n = 24) generated by M. Thrippleton using Matlab code at mjt320/DCE-functions
Detailed info:
Temporal resolution: 0.5 s
Acquisition time: 300 s
AIF: Parker function, starting at t=10s
Noise: SD = 0.001 mM
Arterial delay: none or 5s Since it is challenging to fit all parameters across a wide area of parameter space, data is generated with high SNR.
Reference values: Reference values are the parameters used to generate the data. All combinations of \(v_p\) (0.02, 0.1), \(v_e\) (0.1, 0.2), \(f_p\) (5 to 40 100ml/ml/min) and PS (0.05 to 0.15 per min) are included. A delayed version of the test data was created by shifting the time curves with 5s. This data is labeled as ‘delayed’ and only used with the models that allow the fitting of a delay.
Citation: Code used in Manning et al., Magnetic Resonance in Medicine, 2021 https://doi.org/10.1002/mrm.28833 and Matlab code: mjt320/DCE-functions
Tolerances
\(v_p\): a_tol=0.025, r_tol=0, start=0.01, bounds=(0,1)
\(v_e\): a_tol=0.05, r_tol=0, start=0.2, bounds=(0,1)
\(f_p\): a_tol=5, r_tol=0.1, start=20, bounds=(0,200), units ml/100ml/min
PS: a_tol=0.005, r_tol=0.1, start=0.6, bounds=(0,5), units [/min]
delay: a_tol=0.5, r_tol=0, start=0, bounds=(-10,10), units s
Visualize test data#
To get an impression of the test data, the concentration time curves of the test data are plotted below.
#plot test data
filename = ('../test/DCEmodels/data/2cxm_sd_0.001_delay_0.csv')
# read from CSV to pandas
df1 = pd.read_csv(filename)
no_voxels = len(df1.label)
fig, ax = plt.subplots(1, 1, sharex='col', sharey='row', figsize=(6,6))
for currentvoxel in range(no_voxels):
labelname = 'case_' + str(currentvoxel+1)
testdata = df1[(df1['label']==labelname)]
t = testdata['t'].to_numpy()
t = np.fromstring(t[0], dtype=float, sep=' ')
c = testdata['C_t'].to_numpy()
c = np.fromstring(c[0], dtype=float, sep=' ')
ax.plot(t, c, label=labelname)
ax.set_ylabel('C (mM)', fontsize=14)
ax.set_xlabel('time (ms)', fontsize=14)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0,fontsize=10);
plt.xticks(fontsize=14)
plt.yticks(fontsize=14);
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'] == 'DCEmodels') & (entry['method'] == '2CXM') :
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)
author_list = df.author.unique()
no_authors = len(author_list)
# split delayed and non-delayed data
df['delay'] = df['label'].str.contains('_delayed')
# calculate error between measured and reference values
df['error_ps'] = df['ps_meas'] - df['ps_ref']
df['error_vp'] = df['vp_meas'] - df['vp_ref']
df['error_ve'] = df['ve_meas'] - df['ve_ref']
df['error_fp'] = df['fp_meas'] - df['fp_ref']
# tolerances
tolerances = { 'ps': {'atol' : 0.005, 'rtol': 0.1 }, 'vp': {'atol':0.025, 'rtol':0}, 've': {'atol':0.05, 'rtol':0}, 'fp': {'atol':5, 'rtol':0.1}}
Results#
Non-delayed data#
Some models allow the fit of a delay. For the tests with non-delayed data, the delay was fixed to 0. If the tolerances are not shown, it means that they are outside the limits of the y-axis.
df.head(n=23)[['label','ps_ref','vp_ref', 've_ref','fp_ref']]
| label | ps_ref | vp_ref | ve_ref | fp_ref | |
|---|---|---|---|---|---|
| 0 | case_1 | 0.05 | 0.02 | 0.1 | 5 |
| 1 | case_2 | 0.15 | 0.02 | 0.1 | 5 |
| 2 | case_3 | 0.05 | 0.02 | 0.1 | 25 |
| 3 | case_4 | 0.15 | 0.02 | 0.1 | 25 |
| 4 | case_5 | 0.05 | 0.02 | 0.1 | 40 |
| 5 | case_6 | 0.15 | 0.02 | 0.1 | 40 |
| 6 | case_7 | 0.05 | 0.02 | 0.2 | 5 |
| 7 | case_8 | 0.15 | 0.02 | 0.2 | 5 |
| 8 | case_9 | 0.05 | 0.02 | 0.2 | 25 |
| 9 | case_10 | 0.15 | 0.02 | 0.2 | 25 |
| 10 | case_11 | 0.05 | 0.02 | 0.2 | 40 |
| 11 | case_12 | 0.15 | 0.02 | 0.2 | 40 |
| 12 | case_13 | 0.05 | 0.10 | 0.1 | 5 |
| 13 | case_14 | 0.15 | 0.10 | 0.1 | 5 |
| 14 | case_15 | 0.05 | 0.10 | 0.1 | 25 |
| 15 | case_16 | 0.15 | 0.10 | 0.1 | 25 |
| 16 | case_17 | 0.05 | 0.10 | 0.1 | 40 |
| 17 | case_18 | 0.15 | 0.10 | 0.1 | 40 |
| 18 | case_19 | 0.05 | 0.10 | 0.2 | 5 |
| 19 | case_20 | 0.15 | 0.10 | 0.2 | 5 |
| 20 | case_21 | 0.05 | 0.10 | 0.2 | 25 |
| 21 | case_22 | 0.15 | 0.10 | 0.2 | 25 |
| 22 | case_23 | 0.05 | 0.10 | 0.2 | 40 |
Show code cell source
# set-up styling for seaborn plots
sns.set(font_scale=1.5)
#sns.set_style("whitegrid")
sns.set_style("ticks")
plotopts = {"hue":"author",
"dodge":True,
"s":80,
"height":4,
"aspect":3
}
# Instead of a regular bland-altman plot we opted for a catplot + swarm for these results
# In this way we can appreciate the results of the different contributions per test case better.
# The downside is that the values of the test cases are not obvious.
# This might be something to improve upon at a later stage
make_catplot(x='label', y="error_ps", data=df[~df['delay']],
ylabel="$\Delta$ PS ($min^{-1}$)", **plotopts)
Bias results estimated PS values combined for all voxels
make_catplot(x='label', y="error_vp", data=df[~df['delay']],
ylabel="$v_{p,ref}$ (-)", **plotopts)
make_catplot(x='label', y="error_ve", data=df[~df['delay']],
ylabel="$v_{e,ref}$ (-)", **plotopts)
make_catplot(x='label', y="error_fp", data=df[~df['delay']],
ylabel="$f_{p,ref}$ (-)", **plotopts)
resultsBA = bland_altman_statistics(data=df[~df['delay']],par='error_ps',grouptag='author')
print(resultsBA)
bias stdev LoA lower LoA upper
author
LEK_UoEdinburgh_UK 0.000150 0.004673 -0.009010 0.009309
MB_QBI_UoManchester 0.000061 0.004847 -0.009440 0.009561
MJT_UoEdinburgh_UK 0.000259 0.004817 -0.009183 0.009701
OGJ_OsloU_Norway 0.000061 0.004847 -0.009440 0.009561
Bias results estimated \(v_p\) values combined for all voxels
resultsBA = bland_altman_statistics(data=df[~df['delay']],par='error_vp',grouptag='author')
print(resultsBA)
bias stdev LoA lower LoA upper
author
LEK_UoEdinburgh_UK 0.001673 0.003827 -0.005828 0.009173
MB_QBI_UoManchester 0.000856 0.003520 -0.006044 0.007756
MJT_UoEdinburgh_UK 0.000606 0.003206 -0.005678 0.006889
OGJ_OsloU_Norway 0.000856 0.003520 -0.006044 0.007756
Bias results estimated \(v_e\) values combined for all voxels
resultsBA = bland_altman_statistics(data=df[~df['delay']],par='error_ve',grouptag='author')
print(resultsBA)
bias stdev LoA lower LoA upper
author
LEK_UoEdinburgh_UK -0.001461 0.003401 -0.008127 0.005204
MB_QBI_UoManchester -0.000666 0.003261 -0.007058 0.005725
MJT_UoEdinburgh_UK -0.000492 0.002916 -0.006207 0.005223
OGJ_OsloU_Norway -0.000666 0.003261 -0.007058 0.005725
Bias results estimated \(f_p\) values combined for all voxels
resultsBA = bland_altman_statistics(data=df[~df['delay']],par='error_fp',grouptag='author')
print(resultsBA)
bias stdev LoA lower LoA upper
author
LEK_UoEdinburgh_UK -0.499689 0.535799 -1.549855 0.550476
MB_QBI_UoManchester -0.028328 0.239519 -0.497785 0.441130
MJT_UoEdinburgh_UK 0.003312 0.124339 -0.240393 0.247016
OGJ_OsloU_Norway -0.028328 0.239519 -0.497785 0.441130
Delayed results#
Some contributions allowed the fitting of a delay. For those, additional tests with a delay were performed.
make_catplot(x='label', y="error_ps", data=df[df['delay']],
ylabel="$\Delta$ PS ($min^{-1}$)", **plotopts)
Bias results estimated PS values combined for all voxels
make_catplot(x='label', y="error_vp", data=df[df['delay']],
ylabel="$v_{p,ref}$ (-)", **plotopts)
make_catplot(x='label', y="error_ve", data=df[df['delay']],
ylabel="$v_{e,ref}$ (-)", **plotopts)
make_catplot(x='label', y="error_fp", data=df[df['delay']],
ylabel="$f_{p,ref}$ (-)", **plotopts)
Notes#
Additional notes/remarks