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
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
: a_tol=0.025, r_tol=0, start=0.01, bounds=(0,1) : a_tol=0.05, r_tol=0, start=0.2, bounds=(0,1) : a_tol=5, r_tol=0.1, start=20, bounds=(0,200), units ml/100ml/minPS: 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);
![_images/b0ade7958e1f7adb437543011aba44276ae2ebce5c19b29ddedbcb8763d1f1fc.png](_images/b0ade7958e1f7adb437543011aba44276ae2ebce5c19b29ddedbcb8763d1f1fc.png)
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)
![_images/2431302a4205859f7a53e85291c59aa66d598bd612431a227f5c66f1980b15cf.png](_images/2431302a4205859f7a53e85291c59aa66d598bd612431a227f5c66f1980b15cf.png)
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)
![_images/b5b6cdb247ffafe20c89a89267bbe5af509f857ac9a04fde21e5fa775d40c1a9.png](_images/b5b6cdb247ffafe20c89a89267bbe5af509f857ac9a04fde21e5fa775d40c1a9.png)
make_catplot(x='label', y="error_ve", data=df[~df['delay']],
ylabel="$v_{e,ref}$ (-)", **plotopts)
![_images/c6578f8e45a8dcf582b3aeef3903740e6cac250b4d856696977a2853fb733a56.png](_images/c6578f8e45a8dcf582b3aeef3903740e6cac250b4d856696977a2853fb733a56.png)
make_catplot(x='label', y="error_fp", data=df[~df['delay']],
ylabel="$f_{p,ref}$ (-)", **plotopts)
![_images/f31d1a03dd3d8f2d5c32a8137d6f27b2450968cdb6f95330ee468fe04309603f.png](_images/f31d1a03dd3d8f2d5c32a8137d6f27b2450968cdb6f95330ee468fe04309603f.png)
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
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
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
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)
![_images/3a8cf79e5cdeff9dd1a528b17ff6e48daefbf7eb9fefa2b70e61adcdb014a2f0.png](_images/3a8cf79e5cdeff9dd1a528b17ff6e48daefbf7eb9fefa2b70e61adcdb014a2f0.png)
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)
![_images/3e6bddec822bf90f4493d9d364dc4b21dc5bb6e8a9b1ca52dbab054785ffae47.png](_images/3e6bddec822bf90f4493d9d364dc4b21dc5bb6e8a9b1ca52dbab054785ffae47.png)
make_catplot(x='label', y="error_ve", data=df[df['delay']],
ylabel="$v_{e,ref}$ (-)", **plotopts)
![_images/d24ed75236b70becbbf657d2d3e4eafd68e0f4b9fcfc64541ada46d9767781af.png](_images/d24ed75236b70becbbf657d2d3e4eafd68e0f4b9fcfc64541ada46d9767781af.png)
make_catplot(x='label', y="error_fp", data=df[df['delay']],
ylabel="$f_{p,ref}$ (-)", **plotopts)
![_images/2977af6c626b41fd47b4cf34953b98f4deea438af48aeb901e639eb4bf454709.png](_images/2977af6c626b41fd47b4cf34953b98f4deea438af48aeb901e639eb4bf454709.png)
Notes#
Additional notes/remarks