Two-compartment uptake model (2CUM)#

Hide 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 uptake model data

Source: Concentration-time data (n = 27) 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 to 0.1), \(f_p\) (5 to 40 100ml/ml/min) and PS (0.0001 to 0.25 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)

  • \(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/2cum_sd_0.0025_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=9);   
plt.xticks(fontsize=14)
plt.yticks(fontsize=14);
_images/0a232aadf9d093c96e586fd1dec0f1509c128fb4077b12fd5320b3ed15b73d8c.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'] == '2CUM') :
        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_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 tolerance lines are not shown, it means that they are outside the limits of the y-axis.

df.head(n=27)[['label','ps_ref','vp_ref', 'fp_ref']]
label ps_ref vp_ref fp_ref
0 case_1 0.00001 0.02 5
1 case_2 0.01000 0.02 5
2 case_3 0.02500 0.02 5
3 case_4 0.00001 0.02 25
4 case_5 0.01000 0.02 25
5 case_6 0.02500 0.02 25
6 case_7 0.00001 0.02 40
7 case_8 0.01000 0.02 40
8 case_9 0.02500 0.02 40
9 case_10 0.00001 0.05 5
10 case_11 0.01000 0.05 5
11 case_12 0.02500 0.05 5
12 case_13 0.00001 0.05 25
13 case_14 0.01000 0.05 25
14 case_15 0.02500 0.05 25
15 case_16 0.00001 0.05 40
16 case_17 0.01000 0.05 40
17 case_18 0.02500 0.05 40
18 case_19 0.00001 0.10 5
19 case_20 0.01000 0.10 5
20 case_21 0.02500 0.10 5
21 case_22 0.00001 0.10 25
22 case_23 0.01000 0.10 25
23 case_24 0.02500 0.10 25
24 case_25 0.00001 0.10 40
25 case_26 0.01000 0.10 40
26 case_27 0.02500 0.10 40 
# 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":100,
            "height":4,
            "aspect":3
           }

make_catplot(x='label', y="error_ps", data=df[~df['delay']], 
             ylabel="$\Delta$ PS ($min^{-1}$)", **plotopts)
_images/68d173413633479aa96e2e4a420bf05f168c8cb2998a5b62f887bae15e4d4b00.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/caf1f29709cfeea676b5506e15f17e136ecb81df96312d129bf81da307ca21b0.png 
make_catplot(x='label', y="error_fp", data=df[~df['delay']], 
             ylabel="$f_{p,ref}$ (-)", **plotopts)
_images/dc5d4b128c527561fbdc8e6ed47c9b4d34be33b038e04531953e74087e152be8.png

Bias results estimated PS values combined for all voxels

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.000501  0.000590  -0.001658   0.000656
MJT_UoEdinburgh_UK -0.000518  0.000588  -0.001671   0.000635

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.000340  0.000658  -0.000949   0.001630
MJT_UoEdinburgh_UK  0.000487  0.000586  -0.000662   0.001637

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.427252  0.474663  -1.357592   0.503088
MJT_UoEdinburgh_UK -0.096534  0.220072  -0.527876   0.334807

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/8d5ede94943ce7b93211b24908fe822555e2811f7db1b14e364ca5c54be4bfbb.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/eead3fb2ed233cfe873033071a1fe5f4e7e29daf373dfedf054fc8522de5a737.png 
make_catplot(x='label', y="error_fp", data=df[df['delay']], 
             ylabel="$f_{p,ref}$ (-)", **plotopts)
_images/f7d3835f0287c29b4a58663f9272dc0bfd4575a136cc3bfd4a1965591ec2050a.png

Notes#

Additional notes/remarks

References#