Brain Clearance of Contrast Agent in Intravenous DCE-MRI Is Measurable and Cannot Be Explained by Clearance Across the BBB Alone.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § Methods › DCE‐MRI › Data Extraction and Model Fitting ↔ functions.py, lines 797–922 · score 0.61 · venous impulse response, impulse response functions, deconvolved, arterial, plug, uptake
[2] § Results ↔ rebuttal_functions.py, lines 563–619 · score 0.53 · fast exchange limit, water exchange, BBB, Patlak, model

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Python · 1,562 lines · 58 KB · no license · 1 match

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 2 13:55:59 2021
@author: Martin Kozár
[email hidden]
"""
"""
Code created using:
Python 3.8.10
Matplotlib 3.3.4
Numpy 1.20.1
Pandas 1.2.4
Scipy 1.6.2
"""
import matplotlib.pyplot as plt
from scipy.integrate import cumulative_trapezoid
from scipy.optimize import least_squares
from scipy.signal import convolve, deconvolve
import numpy as np
import os
from multiprocessing import Pool
#from psutil import cpu_count
from multiprocess_function import pool_fun_wrap
from tqdm import tqdm
import copy
import warnings
warnings.simplefilter(action='ignore', category=UserWarning)
import pandas as pd
# %% Helpful wrappers
from functools import wraps
from time import time
def measure(func):
"""
Measures function execution time.
"""
@wraps(func)
def _time_it(*args, **kwargs):
start = int(round(time() * 1000))
try:
return func(*args, **kwargs)
finally:
end_ = int(round(time() * 1000)) - start
print(f"Exe-time of {func.__name__}: {end_ if end_ > 0 else 0} ms")
return _time_it
def filecheck_csv(func):
"""
Checks if a csv file exists and reads it in if it does.
Adds a first positional argument -> filepath : str
"""
@wraps(func)
def _file_exists(filepaths, *args,
index_col = None, header = 'infer', **kwargs):
if all(os.path.isfile(filepath) for filepath in filepaths):
t = tuple(
pd.read_csv(
filepath,
index_col = index_col,
header = header
)
for filepath in filepaths
)
return t[0] if len(t) == 1 else t
else:
return func(filepaths, *args, **kwargs)
return _file_exists
def deconv(data, irf):
N = int(2**(np.ceil(np.log2(len(irf) + len(data)))))
return np.abs(np.fft.ifft(np.fft.fft(data, N) / np.fft.fft(irf, N))), None
# %% Helpful functions
@filecheck_csv
def get_data(filepath_s, brain_loc, slice_data, undersampled):
"""
Checks if a file with all the data exists, if not, create it.
Parameters
----------
brain_loc : list
A list of brain locations from which the data was taken.
Returns
-------
multiDF : pandas.core.frame.DataFrame
Multiindex dataframe where index 0 are brain locations and index 1 are
time points.
"""
dfs = []
udsstr = "_undersampled" if undersampled else ""
# read in the data from 3 .csv files (different brain_loc, same VIF)
for i in range(len(brain_loc)):
filepath_r = f"PvU_remaster/CSVData/ROI_fits_{brain_loc[i]}_June18_signal_clean{udsstr}.csv"
dfs.append(pd.read_csv(filepath_r, decimal = ','))
# join the data into a single dataframe
multiDF = dfs_into_multiDF(dfs, brain_loc, slice_data)
multiDF.to_csv(make_dir(*filepath_s))
return multiDF
def make_dir(filepath):
"""
Checks if a directory exists and creates it if not. Useful when trying to
save a file to a new directory.
Parameters
----------
filepath : str
File path to save the file.
Returns
-------
filepath : str
"""
head, _ = os.path.split(filepath)
if not os.path.isdir(head):
os.makedirs(head, mode = 0o766)
return filepath
@measure
def wrapper_process(inputs, cores):
"""
A function for running concurrent processing. The tricky part of using
Pool().map() is that you have to supply exactly 2 arguments: the function
and a list of all inputs you want to run. I use pool_model_wrap() to then
call model_optimise() with all the arguments from inputs unpacked.
Parameters
----------
inputs : list
A list of all inputs model_optimise requires.
cores : integer
Number of concurrent processes.
Returns
-------
list
A list containing the outputs of running pool_model_wrap with each input.
"""
with Pool(cores) as p:
r = list(tqdm(p.imap(pool_fun_wrap, inputs), total = len(inputs)))
print()
return r
def remove_files(filepaths, delete = False):
"""
Checks if a file exists and removes it if it does.
Parameters
----------
filepaths : list-like
List of filepaths.
delete : bool, optional
The default is False.
"""
if delete:
for filepath in filepaths:
if os.path.isfile(filepath):
print(f"Removing {filepath}.")
os.remove(filepath)
def check_input(user_input: str, method = "isalnum", boolean = True,
ignore_space = True) -> str:
"""
Checks if the user_input is or isn't a chosen type based on the method.
Parameters
----------
user_input : str
method : str, optional
String method that needs to return 'boolean' to pass the test. Accepts:
isalnum, isalpha, isdecimal, isdigit, isidentifier, islower, isnumeric,
isprintable, isspace, istitle, isupper. The default is "isalpha".
boolean : bool, optional
ignore_space : bool, optional
Returns
-------
user_input : str
"""
# string methods that return True/False
str_methods = ["isalnum", "isalpha", "isdecimal", "isdigit", "isidentifier",
"islower", "isnumeric", "isprintable", "isspace", "istitle",
"isupper"]
# check if function parameters are valid
if not isinstance(method, str):
raise TypeError("The 'method' parameter only accepts strings.")
if method not in str_methods:
raise ValueError("Please select a valid built-in Python method for a "
"string that returns True/False, Options are: "
f"{' '.join(str_methods)}. "
f"You chose: {method}")
temp_input = user_input
if ignore_space:
temp_input = temp_input.replace(" ", "")
# apply method to user_input
if getattr(temp_input, method)() != boolean:
raise ValueError(f"Invalid input! Input must cause str.{method}() to"
f" return {boolean}. Your input was: {user_input}")
return user_input
# %% Loading data and simple data statistics
def get_filters(patient_dict, index_data):
"""
Creates a dictionary where keys are patient descriptions and values are all
column names that correspond to that description.
"""
catg_filters = {}
for description, category in patient_dict.items():
catg_filters[description] = ([idx for idx in index_data
if idx.startswith(category)])
return catg_filters
def dfs_into_multiDF(data, brain_loc, slice_data):
"""
Turns a list of dataframes into a single multiindex dataframe.
data: List[Dataframe], brain_loc: List[string],
return multiDF: Dataframe
"""
loaded_data = []
for i, csv_data in enumerate(data):
temp = (csv_data.dropna(how='all')
.dropna(axis='columns')
.set_index("time")
)
tissue = temp.iloc[:(temp.shape[0]//2)]
VIF = temp.iloc[(temp.shape[0]//2):]
if slice_data:
tissue = tissue.iloc[:int(tissue.shape[0]*slice_data)]
VIF = VIF.iloc[:int(VIF.shape[0]*slice_data)]
temp = pd.concat((tissue, VIF),
keys = ["tissue", "VIF"],
names = ["tissue", "time"])
temp = temp.T
loaded_data.append(temp)
multiDF = pd.concat(loaded_data,
keys = brain_loc,
names = ["brain_loc", "ID"])
return multiDF
@filecheck_csv
def statistics(filepath, data, catg_filters):
"""
Summarises data for each patient_desc as columns of stats.
"""
stat_labels = ['mean','std','sem']
methods = [pd.DataFrame.mean, pd.DataFrame.std, pd.DataFrame.sem]
# for groupby using a dictionary, assign each person a category
# c - category ("control", "stroke", "Parkinson's")
# P - list of IDs for category
# p - individual IDs
groupby_filters = {p: c for c, P in catg_filters.items() for p in P}
# applies each method in methods and joins up the results in a single df
stat_df = pd.concat([(data.groupby(["brain_loc", groupby_filters],
level = "ID")
.apply(method)
.round(6)
) for method in methods], keys = stat_labels)
stat_df.to_csv(make_dir(*filepath))
stat_df.T.to_csv(make_dir(f"{os.path.split(*filepath)[0]}/raw_data_statistics_forViewing.csv"))
return stat_df
@filecheck_csv
def assign_group_VIF(filepath, data_indiv, stat_df, catg_filters, brain_loc):
"""
Assigns group-averaged VIF data to a copy of the dataframe with raw data.
"""
groupby_filters = {p: c for c, P in catg_filters.items() for p in P}
data_group = data_indiv.copy()
data_group["VIF"] = (data_indiv["VIF"].groupby(["brain_loc", groupby_filters],
level = "ID")
.transform(np.mean)
.round(6))
data_group.to_csv(make_dir(*filepath))
return data_group
@filecheck_csv
def join_noVIF_data(filepath, data_i):
"""
Remove VIF data from the dataset and save (used when plotting individual curves).
"""
df = data_i.loc[:,"tissue"]
df.to_csv(make_dir(*filepath))
return df
@filecheck_csv
def shift_raw_data(filepath, data_raw, cwd):
from scipy import interpolate
tissue = data_raw.loc[:,"tissue"]
VIF = data_raw.loc[:,"VIF"]
#the index value at which the bolus arrives
BATtissue = pd.read_csv(os.path.join(cwd, 'PvU_remaster', 'CSVData', 'BATtissue.csv'), header = None)
BATvif = pd.read_csv(os.path.join(cwd, 'PvU_remaster', 'CSVData', 'BATvif.csv'), header = None)
#create an array of time points for each person
tissueTime = np.tile(tissue.columns.values.astype(float), (tissue.shape[0],1))
#difference between vif and tissue BAT (in units of array index)
timeDiff = (BATvif - BATtissue).to_numpy()
#calculate necessary timeshift
tissueTimeSftd = tissueTime - timeDiff * 7.6
dfTisTimeSftd = pd.DataFrame(tissueTimeSftd, index = tissue.index, columns = tissue.columns)
#df where left is tissue signal, right shifted time
df = pd.concat([tissue, dfTisTimeSftd], axis = 1, keys = ["tissue", "shift"])
def apply_interp(row, x):
"""
Interpolate the signal values to a shift. Expects a DataFrame where each
row is separated by a multiindex into tissue signal and shifted timepoints.
Uses scipy.interpolate
Parameters
----------
row : Series
Row of the dataframe supplied by apply function
x : array-like
Original time points when the signal was taken
"""
tck = interpolate.splrep(x, row["tissue"], s = 0)
return interpolate.splev(row["shift"], tck, der = 0)
tissueSftd = df.apply(apply_interp, args = (tissue.columns.values.astype(float), ), axis = 1, result_type='expand')
tissueSftd.columns = tissue.columns
output = pd.concat([tissueSftd, VIF], axis = 1, keys = ["tissue", "VIF"], names = ['tissue', 'time'])
output.to_csv(make_dir(*filepath))
return output
# %% Optimisation (data fitting)
@filecheck_csv
def model_setup(fps_predict, data_indiv_VIF, model_type, time_data,
time_windows, brain_loc, weighted = True):
"""
Sets up the combinations of settings for how to fit model to measurements.
Each setting then gets passed to the wrapper_process for multiprocessing.
The output then gets concatenated into 3 dataframes and returned.
Parameters
----------
fps_predict : str
Filepath to save the prediction values.
data : tuple
Contains the individual-VIF and group-VIF data.
model_type : str
time_data : numpy.ndarray
The timestamps at which the data was measured.
time_window : list
List containing the cut-off time stamp for weighting in seconds.
brain_loc : list
Strings describing the brain locations where measurements were taken.
weighted : bool, optional
Toggle weighting for the residual function. The default is True.
Returns
-------
concatd_data : tuple
Contains (values predicted by model, residual function, constant
values fitted to the model) - each element is a DataFrame.
"""
print("Setting up model.")
if not weighted:
time_windows = [(0, 0)]
model_settings = []
concat_keys = []
for time_window in time_windows:
model_settings.append([model_optimise, data_indiv_VIF, model_type,
time_data, time_window, brain_loc])
concat_keys.append(str(time_window))
print("Done.", end = "\n\n")
# choose number of cores to use
L = len(model_settings)
if L > 6:
cores = 6
elif 6 > L > 1:
for i in range(L, 1, -1):
if L % i == 0:
cores = i
else:
cores = 1
print(f"Fitting model to data using {cores} cores.")
fitted_data_disjoined = wrapper_process(model_settings, cores)
concatd_data = tuple(pd.concat(df,
keys = concat_keys,
names = ["time_window","brain_loc","ID"])
for df in zip(*fitted_data_disjoined))
print("Done.", end = "\n\n")
# calculate v_p and k_trans from 2CUptake predictions
if model_type == "2CUptake":
pred_const_data = concatd_data[2]
pred_const_data["v_p"] = pred_const_data["F_p"] * pred_const_data["T_p"]
pred_const_data["k_trans"] = pred_const_data["F_p"] * pred_const_data["E"]
if model_type == "2CExchange":
pred_const_data = concatd_data[2]
#sourbron perfusion eq. [9]
pred_const_data["T_B"] = np.reciprocal(pred_const_data["K_plus"] - pred_const_data["E_minus"] * (pred_const_data["K_plus"] - pred_const_data["K_minus"]))
#[10]
pred_const_data["T_E"] = np.reciprocal(pred_const_data["T_B"] * pred_const_data["K_plus"] * pred_const_data["K_minus"])
#[11]
pred_const_data["T_P"] = np.reciprocal(pred_const_data["K_plus"] + pred_const_data["K_minus"] - np.reciprocal(pred_const_data["T_E"]))
#next 3 are [12]
pred_const_data["v_p"] = pred_const_data["F_P"] * pred_const_data["T_B"]
pred_const_data["PS"] = pred_const_data["F_P"] * (pred_const_data["T_B"] / pred_const_data["T_P"] - 1)
pred_const_data["v_e"] = pred_const_data["PS"] * pred_const_data["T_E"]
print("Saving model data.")
for data, filepath in zip(concatd_data, fps_predict):
data.to_csv(filepath)
#predictions, resid_df, const_pred, x0_df
return concatd_data
def model_optimise(data, model_type, time_data, time_window, brain_loc):
"""
Sets up empty dataframes and fills them up with the model values.
Parameters
----------
data : pandas.core.frame.DataFrame
model_type : str
time_data : numpy.ndarray
The timestamps at which the data was measured.
time_window : list
List containing the cut-off time stamp for weighting in seconds.
brain_loc : list
Strings describing the brain locations where measurements were taken.
Returns
-------
predictions : pandas.core.frame.DataFrame
Values predicted by model.
resid_df : pandas.core.frame.DataFrame
Residual function.
const_pred : pandas.core.frame.DataFrame
Constant values fitted to the model.
"""
#x0 - array of starting points
const_list, x0, bounds = _set_model_inputs(model_type)
person_list = data.index.get_level_values("ID").unique()
midx = pd.MultiIndex.from_product([brain_loc, person_list],
names = ["brain_loc","ID"])
#prediction curves
predictions = pd.DataFrame(index=midx, columns=time_data, dtype = np.float64)
#fitted parameters
const_pred = pd.DataFrame(index=pd.MultiIndex.from_product([brain_loc, person_list],
names = ["brain_loc", "ID"]),
columns=const_list)
resid_df = pd.DataFrame(index=midx, columns=time_data)
#starting parameters
x0_df = pd.DataFrame(index=midx, columns = const_list)
(predictions[predictions.columns],
resid_df[resid_df.columns],
const_pred[const_pred.columns],
x0_df[x0_df.columns]) = zip(*data.apply(_scipy_fun_wrap,
axis = 1,
args = (residual,
x0,
time_data/60,
time_window,
model_type
),
method = "trf",
bounds = bounds
)
)
return predictions, resid_df, const_pred, x0_df
def residual(x, VIF, measured_data, time_data, time_window,
model_type, integral):
"""
Calculates the difference between model and measured values.
Parameters
----------
x : array_like
Contains the variables to be determined.
VIF : array_like
Contains input data used for prediction.
measured_data : array_like
Contains measured data to compare the prediction against.
time_data : array_like
time_data points at which the data was measured.
time_window : tuple
'weighted = True' sets data in the period time_window to zero.
model_type : string
Determines which model to use.
Returns
-------
array_like
Array-like with differences between predicted and measured values.
"""
if model_type == "Intravascular":
#x -> v_p
temp = measured_data - x[0] * VIF
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "Patlak":
#x -> v_p, k_trans
temp = measured_data - (x[0]*VIF + x[1]*integral)
#temp[60 < temp.index < 280] = 0
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "2CUptake":
#x -> F_p, T_p, E
# the search function in scipy searches values outside of float64
#so just ignore the warnings if necessary
# with np.errstate(over='ignore', invalid='ignore'):
R = np.exp(-time_data/x[1]) + (x[2] * (1 - np.exp(-time_data/x[1])))
convolution = convolve(x[0] * R, VIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
temp = measured_data - (convolution[:hl+1])
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "2CUptakeDisp":
VIF = VIF.to_numpy()
time_data = time_data.to_numpy()
#x -> F_p, T_p, E, T_v
# the search function in scipy searches values outside of float64
#so just ignore the warnings if necessary
# with np.errstate(over='ignore', invalid='ignore'):
# this is a time array for the impulse response function. It needs to be longer than time_data
time_irf = np.arange(0, time_data[1]*(2*len(VIF)-2), time_data[1])
#this is the venous impulse response function
irf_v = np.exp(-time_irf/x[3])/(x[3])
#this is the length we need to pad the VIF with zeroes
lengthtopad = len(time_irf)-len(VIF)
#this is the padded VIF. You will probably need to change how this is done in python.
VIF_pad = np.concatenate((VIF, np.tile(0, lengthtopad)))
#this should provide the full deconvolution
C_p, _ = deconv(VIF_pad, irf_v)
C_p = C_p / time_data[1]
#this crops the full deconvolution to the bit we are interested in
C_p2 = C_p[0:len(VIF)]
#this defines the irf of the cappilary bed
irf_c = np.exp(-time_irf/x[1])/(x[1])
#this bit pads the cappilary concentration with zeros
C_p2_padded = np.concatenate((C_p2, np.tile(0, lengthtopad)))
#this deconvolves the cappilary concentration to the arterial concentration
C_a, _ = deconv(C_p2_padded, irf_c)
C_a = C_a / time_data[1]
#this crops the arterial concentration
AIF = C_a[0:len(VIF)]
R = irf_c[0:len(VIF)] + (x[2] * (1 - np.exp(-time_data/x[1])))
convolution = convolve(x[0] * R, AIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
temp = measured_data - (convolution[:hl+1])
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "2CUptakePlug":
VIF = VIF.to_numpy()
time_data = time_data.to_numpy()
#x -> F_p, T_p, E, T_v, T_l
# the search function in scipy searches values outside of float64
#so just ignore the warnings if necessary
# with np.errstate(over='ignore', invalid='ignore'):
# this is a time array for the impulse response function. It needs to be longer than time_data
time_irf = np.arange(0, time_data[1]*(2*len(VIF)-2), time_data[1])
#this is the venous impulse response function
irf_v = np.exp(-time_irf/x[3])/(x[3])
#this is the length we need to pad the VIF with zeroes
lengthtopad = len(time_irf)-len(VIF)
#this is the padded VIF. You will probably need to change how this is done in python.
VIF_pad = np.concatenate((VIF, np.tile(0, lengthtopad)))
#this should provide the full deconvolution
C_p, _ = deconv(VIF_pad, irf_v)
C_p = C_p / time_data[1]
#this crops the full deconvolution to the bit we are interested in
C_p2 = C_p[0:len(VIF)]
#this defines the irf of the cappilary bed
irf_c = np.heaviside(x[1] - time_irf, 0)*np.exp(-time_irf/x[4])/(x[1])
#this bit pads the cappilary concentration with zeros
C_p2_padded = np.concatenate((C_p2, np.tile(0, lengthtopad)))
#this deconvolves the cappilary concentration to the arterial concentration
C_a, _ = deconv(C_p2_padded, irf_c)
C_a = C_a / time_data[1]
#this crops the arterial concentration
AIF = C_a[0:len(VIF)]
k_trans = x[2] * x[0]
irf_e = (k_trans / x[2]) * (1 - np.exp(-time_data / x[4])) * np.heaviside(x[1] - time_data, 0) + k_trans * np.heaviside(time_data - x[1], 0)
h_p = irf_c[0:len(VIF)]
v_p = x[0] * x[1]
H = h_p * v_p + irf_e
convolution = convolve(H, AIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
temp = measured_data - (convolution[:hl+1])
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "ExTofts":
#x -> v_p, k_trans, v_e
#Extended Tofts is the Tofts model (TM) plus intravascular component
convolution = convolve(np.exp(-time_data*x[1]/x[2]), VIF)*(time_data[2]-time_data[1])
#TM = x[1] * convolution[:round(convolution.size/2)+1]
TM = x[1] * convolution[:len(VIF)]
#concentration curve
C = x[0]*VIF + TM
temp = measured_data - C
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "Clearance":
#x -> v_p, k_trans, k
convolution = convolve(np.exp(-time_data*x[2]), VIF)*(time_data[2]-time_data[1])
TM = x[1] * convolution[:round(convolution.size/2)+1]
#concentration curve
C = x[0]*VIF + TM
temp = measured_data - C
temp[temp.index.astype(float) < time_window] = 0
return temp
elif model_type == "2CExchange":
#x -> F_P, E_minus, K_plus, K_minus
# R = e^(-t*K+) + (E- * (e^(-t*K-) - e^(-t*K+)))
R = np.exp(-time_data*x[2]) + (x[1] * (np.exp(-time_data * x[3]) - np.exp(-time_data * x[2])))
# (F_P * R) (X) VIF
convolution = convolve(x[0] * R, VIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
temp = measured_data - (convolution[:hl+1])
temp[temp.index.astype(float) < time_window] = 0
return temp
def _set_model_inputs(model_type):
"""
Sets initial guess and bounds based on model_type (different number of
parameters) and zero_ktrans (sets k_trans to zero).
Parameters
----------
model_type : str
Returns
-------
const_list : list
List of names of the fitted parameters.
x0 : numpy.ndarray
Initial guess for optimisation function.
bounds : 2-tuple containing lists
Sets the boundaries on the fitted parameters.
"""
#x0 is a list of np arrays
if model_type == "Intravascular":
const_list = ["v_p"]
x0 = [np.array([0.03])]
bounds = ([-np.inf], [np.inf])
elif model_type == "Patlak":
const_list = ["v_p", "k_trans"]
x0 = list(np.array([0.03, x]) for x in np.linspace(0.0001, 0.005, 20))
bounds = ([-np.inf,-np.inf], [np.inf,np.inf])
elif model_type == "2CUptake":
const_list = ["F_p", "T_p", "E"]
x0 = list(np.array([x, 0.1, 0.00125]) for x in np.linspace(0.05, 1, 20))
bounds = ([0,0,0],[1,0.3,1])
elif model_type == "2CUptakeDisp":
const_list = ["F_p", "T_p", "E", "T_v"]
x0 = list(np.array([x, 0.1, 0.00125, 0.1]) for x in np.linspace(0.05, 1, 20))
bounds = ([0,0,0,0],[1,0.3,1,0.3])
elif model_type == "2CUptakePlug":
const_list = ["F_p", "T_p", "E", "T_v", "T_l"]
x0 = list(np.array([x, 0.1, 0.00125, 0.1, 0.1]) for x in np.linspace(0.05, 1, 20))
bounds = ([0,0,0,0,0],[1,0.3,1,0.3,0.3])
#v_e limits based on DOI 10.1002/mrm.25793
elif model_type == "ExTofts":
const_list = ["v_p", "k_trans", "v_e"]
x0 = list(np.array([0.03, x, 0.1]) for x in np.linspace(0.0001, 0.005, 20))
bounds = ([0,-np.inf,0.0], [1,np.inf,1])
elif model_type == "2CExchange":
const_list = ["F_P", "E_minus", "K_plus", "K_minus"]
x0 = list(np.array([x,1,1,0.5]) for x in np.linspace(0.05, 1, 20))
bounds = ([0,0,0,0],[5,100,100,100])
else:
print(f"model_type = {model_type}")
#raise ValueError('model_type invalid input, valid inputs are: "Patlak", "2CUptake"')
return const_list, x0, bounds
def _scipy_fun_wrap(row, fun, x0, time_data, time_window, model_type, **kwargs):
"""
Calls the function scipy.optimize.least_squares to do the model fitting and
uses the results to calculate measurement prediction based on VIF.
Parameters
----------
row : pandas Series
Individual rows given by .apply() function from pandas.
fun : function
The function that calculates the residual error between model and measurement.
x0 : array-like
Initial guess of the fitted constants.
time_data : array-like
Time stamps at which the data was taken.
time_window : integer
If weighted == True, residuals up to this point get weighted to zero.
model_type : string
Which model to use for fitting and predicting.
**kwargs :
Extra params for least_squares().
Returns
-------
predictions: array-like
Array containing predictions after fitting.
residual: array-like
Residuals after fitting is complete.
constants: array-like
Fitted values for constants.
xnot: array-like
Starting values with lowest sum of residuals.
"""
# add a loop over several x0 starting points and choose lowest residual sum
tissue, VIF = row["tissue"], row["VIF"]
integral = None
if model_type == "Patlak":
VIF_copy = VIF.copy()
VIF_copy[VIF_copy.index.astype(float) < time_window] = 0
integral = cumulative_trapezoid(
VIF_copy,
time_data,
initial=0
)
residuals_list = []
#loop over all starting points to find global minimum
for xnot in x0:
res_lsq = least_squares(fun, xnot,
ftol=1e-12, xtol=1e-12, gtol=1e-12,
args = (VIF,
tissue,
time_data,
time_window,
model_type,
integral),
**kwargs
)
residuals_list.append((res_lsq.fun.sum(), res_lsq, xnot))
#get the result with minimum residual sum
res_lsq, xnot = min(residuals_list, key = lambda x:x[0])[1:]
prediction = model_predict(res_lsq.x, VIF, time_data, integral, model_type)
return prediction, res_lsq.fun, res_lsq.x, xnot
def model_predict(x, VIF, time_data, integral, model_type):
"""
Predicts measured value based on VIF and fitted constants.
Parameters
----------
x : array_like
List of constants in the same order as in residual().
VIF : array_like
Contains input data used for prediction.
time_data : array_like
time_data points at which the data was measured.
integral : array_like or None
In Patlak model we do a single integral instead of many convolutions.
model_type : string
Determines which model to use.
Returns
-------
array_like
Array-like with the predictions based on VIF.
"""
if model_type == "Intravascular":
#x -> v_p
return x[0]*VIF
elif model_type == "Patlak":
#x -> v_p, k_trans
return x[0]*VIF + x[1]*integral
elif model_type == "2CUptake":
#x -> F_p, T_p, E
R = np.exp(-time_data/x[1]) + (x[2] * (1 - np.exp(-time_data/x[1])))
convolution = convolve(x[0] * R, VIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
return convolution[:hl+1]
elif model_type == "2CUptakeDisp":
#x -> F_p, T_p, E, T_v
VIF = VIF.to_numpy()
time_data = time_data.to_numpy()
#x -> F_p, T_p, E, T_v
# the search function in scipy searches values outside of float64
#so just ignore the warnings if necessary
# with np.errstate(over='ignore', invalid='ignore'):
# this is a time array for the impulse response function. It needs to be longer than time_data
time_irf = np.arange(0, time_data[1]*(2*len(VIF)-2), time_data[1])
#this is the venous impulse response function
irf_v = np.exp(-time_irf/x[3])/(x[3])
#this is the length we need to pad the VIF with zeroes
lengthtopad = len(time_irf)-len(VIF)
#this is the padded VIF. You will probably need to change how this is done in python.
VIF_pad = np.concatenate((VIF, np.tile(0, lengthtopad)))
#this should provide the full deconvolution
C_p, _ = deconv(VIF_pad, irf_v)
C_p = C_p / time_data[1]
#this crops the full deconvolution to the bit we are interested in
C_p2 = C_p[0:len(VIF)]
#this defines the irf of the cappilary bed
irf_c = np.exp(-time_irf/x[1])/(x[1])
#this bit pads the cappilary concentration with zeros
C_p2_padded = np.concatenate((C_p2, np.tile(0, lengthtopad)))
#this deconvolves the cappilary concentration to the arterial concentration
C_a, _ = deconv(C_p2_padded, irf_c)
C_a = C_a / time_data[1]
#this crops the arterial concentration
AIF = C_a[0:len(VIF)]
R = irf_c[0:len(VIF)] + (x[2] * (1 - np.exp(-time_data/x[1])))
convolution = convolve(x[0] * R, AIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
return convolution[:hl+1]
elif model_type == "2CUptakePlug":
VIF = VIF.to_numpy()
time_data = time_data.to_numpy()
#x -> F_p, T_p, E, T_v, T_l
# the search function in scipy searches values outside of float64
#so just ignore the warnings if necessary
# with np.errstate(over='ignore', invalid='ignore'):
# this is a time array for the impulse response function. It needs to be longer than time_data
time_irf = np.arange(0, time_data[1]*(2*len(VIF)-2), time_data[1])
#this is the venous impulse response function
irf_v = np.exp(-time_irf/x[3])/(x[3])
#this is the length we need to pad the VIF with zeroes
lengthtopad = len(time_irf)-len(VIF)
#this is the padded VIF. You will probably need to change how this is done in python.
VIF_pad = np.concatenate((VIF, np.tile(0, lengthtopad)))
#this should provide the full deconvolution
C_p, _ = deconv(VIF_pad, irf_v)
C_p = C_p / time_data[1]
#this crops the full deconvolution to the bit we are interested in
C_p2 = C_p[0:len(VIF)]
#this defines the irf of the cappilary bed
irf_c = np.heaviside(x[1] - time_irf, 0)*np.exp(-time_irf/x[4])/(x[1])
#this bit pads the cappilary concentration with zeros
C_p2_padded = np.concatenate((C_p2, np.tile(0, lengthtopad)))
#this deconvolves the cappilary concentration to the arterial concentration
C_a, _ = deconv(C_p2_padded, irf_c)
C_a = C_a / time_data[1]
#this crops the arterial concentration
AIF = C_a[0:len(VIF)]
k_trans = x[2] * x[0]
irf_e = (k_trans / x[2]) * (1 - np.exp(-time_data / x[4])) * np.heaviside(x[1] - time_data, 0) + k_trans * np.heaviside(time_data - x[1], 0)
h_p = irf_c[0:len(VIF)]
v_p = x[0] * x[1]
H = h_p * v_p + irf_e
convolution = convolve(H, AIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
return convolution[:hl+1]
elif model_type == "ExTofts":
#x -> v_p, k_trans, v_e
convolution = convolve(np.exp(-time_data*x[1]/x[2]), VIF)*(time_data[2]-time_data[1])
#TM = x[1] * convolution[:round(convolution.size/2)+1]
TM = x[1] * convolution[:len(VIF)]
#concentration curve
return x[0]*VIF + TM
elif model_type == "2CExchange":
#x -> F_P, E_minus, K_plus, K_minus
R = np.exp(-time_data * x[2]) + (x[1] * (np.exp(-time_data * x[3]) - np.exp(-time_data * x[2])))
convolution = convolve(x[0] * R, VIF)*(time_data[2]-time_data[1])
hl = round(convolution.size/2)
return convolution[:hl+1]
# %% Plotting
def pred_plt_wrap(model_type, df_pred, df_meas, dir_plt):
"""
Sets up a dataframe with each row having predicted and the appropriate measured
values next to each other, then feeds each row for plotting.
Parameters
----------
model_type : str
df_pred : pandas.core.frame.DataFrame
Dataframe containing prediction values.
df_meas : pandas.core.frame.DataFrame
Dataframe containing measured values. Expected to have layers
compatible with df_pred for joining.
Returns
-------
None.
"""
df_pred.columns = pd.MultiIndex.from_product([['predictions'],
df_pred.columns],
names = ["data", "time"])
df_meas.columns = pd.MultiIndex.from_product([['measurements'],
df_meas.columns],
names = ["data", "time"])
df_P_M = (df_pred.reset_index(("time_window"), drop = False)
.join(df_meas)
.set_index(["time_window"], append = True)
)
df_P_M = df_P_M.reorder_levels(["time_window", "brain_loc", "ID"])
# it might make sense to use threading instead since this is a lot of saving? haven't tested the speed
with Pool(6) as p:
#remove .head(5) to plot all rows
list(tqdm(p.imap(prediction_plot,
((model_type, dir_plt, x) for x in df_P_M.iterrows())),
total = df_P_M.shape[0]))
print()
def prediction_plot(row):
"""
Creates a lineplot and labels it using data from row, saves and clears figure.
Parameters
----------
row : tuple
Tuple of (model_type, cwd, (index labels and pandas Series)).
Returns
-------
None.
"""
model_type = row[0]
dir_plt = row[1]
labels = row[2][0]
data = row[2][1]
data.rename(index=lambda val: int(float(val)), level = "time", inplace = True)
data["predictions"].plot()
data["measurements"].plot()
dirdesc = '/'.join(str(x) for x in labels[0:-1])
plt.title(f"{model_type} Time window: {labels[0]}, Location: {labels[1]}, ID: {labels[2]}")
plt.xlabel('Time (s)')
plt.ylabel('Contrast agent concentration (mM)')
plt.legend(["prediction", "measurement"])
dirpath = (f"{dir_plt}/graphs/{model_type}/" + dirdesc)
if not os.path.isdir(dirpath):
os.makedirs(dirpath, mode = 0o766)
plt.savefig(f"{dirpath}/pred_vs_meas_{model_type}_{labels}.jpg")
plt.clf()
def boxplot_const_wrap(pred_const_data, model_type, time_windows, dir_plt):
"""
Plots a boxplot for each constant and brain location.
Parameters
----------
pred_const_data : pandas.core.frame.DataFrame
model_type : str
Returns
-------
None.
"""
dpi = 96
res = (1920,1080)
unstack_lvls = ['time_window', 'patient_catg', 'ID']
# with ktrans -> Patlak
plot_data = (pred_const_data.unstack(level = unstack_lvls)
.stack(level = 0)
)
# unique values describing the plotted data
uniq_cols = plot_data.columns.get_level_values("time_window").unique()
if uniq_cols.dtype == int:
uniq_cols = uniq_cols.astype(int)
for row in plot_data.iterrows():
data = row[1]
const_label = row[0][1]
brain_loc = row[0][0]
fig, axs = plt.subplots(1, uniq_cols.size,
num = f'{const_label} {brain_loc}',
sharey = "row",
figsize = (res[0]/dpi, res[1]/dpi),
dpi = dpi)
data = data.unstack(level = "time_window")
data = data.reindex(["control","stroke","Parkinson's"], level = "patient_catg")
categories = data.index.get_level_values("patient_catg").unique()
axs = data.boxplot(by = "patient_catg",
ax = axs,
grid = False)
for i, window in enumerate(uniq_cols):
for w, description in enumerate(categories):
scat_data = data.loc[description, window]
#add a random jiggle to the scatter points
x = np.random.normal(w+1, 0.05, len(scat_data))
axs[i].scatter(x = x,
y = scat_data,
alpha = 0.4)
axs[i].set_title(f"Cutoff: {window}s")
# pandas adds the input of "by=" as xlabel, this removes it
axs[i].set(xlabel = "")
axs[i].grid(axis = "y")
if const_label == "v_e":
axs[i].set_ylim([0, 0.06])
fig.suptitle(f"Boxplots of {model_type} predictions at various time windows for "
f"constant: {const_label} in: {brain_loc}")
if const_label == "k_trans":
axs[0].set_ylabel("$\mathregular{K_{trans}}$ [min⁻¹]")
elif const_label == "v_p":
axs[0].set_ylabel("$\mathregular{v_p}$ [ml⁻¹ tissue]")
elif const_label == "E":
axs[0].set_ylabel("E")
elif const_label == "F_p":
axs[0].set_ylabel("$\mathregular{F_p}$ [ml min⁻¹]")
elif const_label == "T_p":
axs[0].set_ylabel("$\mathregular{T_p}$ [min⁻¹]")
elif const_label == "v_e":
axs[0].set_ylabel("$\mathregular{v_e}$ [ml cc⁻¹ tissue]")
elif const_label == "E_minus":
axs[0].set_ylabel("$\mathregular{E-}$")
elif const_label == "K_minus":
axs[0].set_ylabel("$\mathregular{K-}$ [min⁻¹]")
elif const_label == "K_plus":
axs[0].set_ylabel("$\mathregular{K+}$ [min⁻¹]")
dirpath = f"{dir_plt}/boxplots_constants/{model_type}"
if not os.path.isdir(dirpath):
os.makedirs(dirpath, mode = 0o766)
filepath = f"{dirpath}/box_{model_type}_{const_label}_{brain_loc}.png"
if not os.path.isfile(filepath):
fig.savefig(filepath)
plt.close(fig)
def boxplot_AIC_wrap(AICs, pred_const_data, descriptions, brain_locs, models,
undersampled):
AICs_lst = copy.deepcopy(AICs)
dpi = 96
res = (1920,1080)
uniq_cols = pred_const_data.index.get_level_values("time_window").unique()
#reorder so I get index with patient category labels with the same order as AICs
midx = (pred_const_data.copy()
.unstack("patient_catg")
.reindex(AICs[0].index)
.stack("patient_catg")
).index
for i in range(len(AICs_lst)):
AICs_lst[i].index = midx
data = pd.concat(AICs_lst, keys = tuple(models)).unstack("time_window")
print(data)
temp = list(data.index.names)
temp[0] = "model"
data.index.names = temp
data = data.reorder_levels(["patient_catg","brain_loc","model","ID"])
s = "_undersampled" if undersampled else ""
dirpath = f"{os.getcwd()}/PvUgraphs&tables{s}/boxplots_AIC/"
models_in_order = pd.unique(data.index.get_level_values('model'))
# Convert only that level to a categorical with a frozen order
data.index = data.index.set_levels(
pd.CategoricalIndex(models_in_order, categories=models_in_order, ordered=True),
level='model'
)
for desc in descriptions:
for brain_loc in brain_locs:
fig, axs = plt.subplots(1, uniq_cols.size,
num = f'{desc} {brain_loc}',
sharey = "row",
figsize = (res[0]/dpi, res[1]/dpi),
dpi = dpi)
#the groupby inside boxplot sorts the data which is annoying
axs = data.loc[(desc, brain_loc)].boxplot(by = "model",
ax = axs,
rot = 45,
grid = False)
for j, ax in enumerate(axs):
plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")
ax.set_title(f"Window: {uniq_cols[j]}s")
ax.set(xlabel = "")
ax.grid(axis = "y")
for w, model in enumerate(models):
scat_data = data.loc[(desc, brain_loc, model), uniq_cols[j]]
x = np.random.normal(w+1, 0.05, len(scat_data))
ax.scatter(x = x,
y = scat_data,
alpha = 0.4)
axs[0].set_ylabel("AIC")
fig.suptitle(f"Akaike information criterion for fits at different time windows in {brain_loc} for {desc}")
dirpath_d = f"{dirpath}{desc}/"
filepath = make_dir(dirpath_d) + f"AIC_model_comparison_{brain_loc}.png"
#only closes the figure when saving
#if not os.path.isfile(filepath):
print(f"Saving {filepath}")
plt.savefig(filepath)
plt.close(fig)
# %% Statistics on prediction values
def add_catg_label(pred_const_data, catg_filters, model_type):
"""
Adds a multiindex level for the 3 patient categories to separate IDs by
category easier, useful when doing groupby and statistics per category.
Parameters
----------
pred_const_data : pandas.core.frame.DataFrame
catg_filters : dict
model_type : str
Returns
-------
pred_const_data : pandas.core.frame.DataFrame
Dataframe with the new midx level added.
"""
# get tuples with patient category for each patient ID
midx_arr = [(p, c) for c, P in catg_filters.items() for p in P]
# create multiindex
midx = pd.MultiIndex.from_tuples(midx_arr, names = ["ID", "patient_catg"])
# remove index levels except for patient ID
unstack_labels = pred_const_data.index.names[:-1]
pred_const_data = pred_const_data.unstack(unstack_labels)
# replace index values, essentially adding the new layer
pred_const_data.index = midx
# return previous index levels
pred_const_data = pred_const_data.stack(unstack_labels)
# reorder index levels back to original order
idcs = pred_const_data.index.names
pred_const_data = pred_const_data.reorder_levels(idcs[2:] + idcs[0:2][::-1])
return pred_const_data
@filecheck_csv
def stats_on_constants(filepath, pred_const_data):
"""
Groups patient data and applies a list of aggregate functions. Then adds
extra columns for standard deviation of a agg result divided by mean of
results in the group.
Parameters
----------
filepath : list
pred_const_data : pandas.core.frame.DataFrame
Dataframe containing fitted constant values.
Returns
-------
pred_const_stats : pandas.core.frame.DataFrame
Dataframe containing a summary of fitted constant values.
"""
methods = ["mean", "std", "sem"]
pred_const_stats = (pred_const_data.groupby(pred_const_data.index.names[:-1])
.agg(methods))
pred_const_stats.rename_axis(columns = ["constant","aggregate"], inplace = True)
for val in pred_const_stats.columns.unique(level = 0):
pred_const_stats[(val, "std/mean")] = (pred_const_stats[(val, "std")].divide(pred_const_stats[(val,"mean")])
.fillna(0))
# reorder columns
pretty_order_0 = pred_const_stats.columns.unique(level = 'constant')
pretty_order_1 = ["mean", "std", "sem", "std/mean"]
pred_const_stats = (pred_const_stats.reindex(columns = pretty_order_0,
level = "constant")
.reindex(columns = pretty_order_1,
level = "aggregate"))
pred_const_stats.to_csv(make_dir(*filepath))
return pred_const_stats
@filecheck_csv
def get_AIC(filepath, resids, model_type):
"""
Calculate the Akaike information criterion using ΔAIC = 2*k + n*ln(RSS)
where RSS = residual sum of squares, n is the number of samples and k
number of parameters of the model.
Parameters
----------
filepath : str
resids : pd.DataFrame
Array of residuals.
model_type : str
Returns
-------
resids_AIC : pd.DataFrame
"""
#model_type : k
num_model_params = {"Intravascular": 1,
"Patlak": 2,
"2CUptake": 3,
"ExTofts": 3,
"2CExchange": 4,
"2CUptakeDisp": 4,
"2CUptakePlug": 5}
sum_list = []
for cutoff in resids.index.get_level_values("time_window").unique():
#sum up the residuals corresponding to a window -> Series
window_data_slice = resids.loc[cutoff].pow(2).sum(axis = 1)
df_slice = pd.DataFrame(window_data_slice, columns = ["sum"])
#add the number of elements with index value greater than cutoff
df_slice["num_elements"] = (resids.columns.values.astype(float) >= int(cutoff)).sum()
sum_list.append(df_slice)
#this WILL store AIC information, doesn't have it yet
resids_AIC = pd.concat(sum_list)
#return time_window index information
resids_AIC.index = resids.index
k = num_model_params[model_type]
#2 * number of parameters + number of data points
# https://www.sciencedirect.com/science/article/pii/S2468042719300508
resids_AIC["AIC"] = 2*k + resids_AIC["num_elements"]*np.log(resids_AIC["sum"]/resids_AIC["num_elements"])
resids_AIC.to_csv(make_dir(*filepath))
return resids_AIC
@filecheck_csv
def get_AICc(filepath, resids, model_type):
"""
https://pmc.ncbi.nlm.nih.gov/articles/PMC3291742/
Calculate the Akaike information criterion with small sample correction using ΔAIC = 2*k + n*ln(RSS) + 2*k*(k+1)/(n-k-1)
where RSS = residual sum of squares, n is the number of samples and k
number of parameters of the model.
Parameters
----------
filepath : str
resids : pd.DataFrame
Array of residuals.
model_type : str
Returns
-------
resids_AIC : pd.DataFrame
"""
#model_type : k
num_model_params = {"Intravascular": 1,
"Patlak": 2,
"2CUptake": 3,
"ExTofts": 3,
"2CExchange": 4,
"2CUptakeDisp": 4,
"2CUptakePlug": 5}
sum_list = []
for cutoff in resids.index.get_level_values("time_window").unique():
#sum up the residuals corresponding to a window -> Series
window_data_slice = resids.loc[cutoff].pow(2).sum(axis = 1)
df_slice = pd.DataFrame(window_data_slice, columns = ["sum"])
#add the number of elements with index value greater than cutoff
df_slice["num_elements"] = (resids.columns.values.astype(float) >= int(cutoff)).sum()
sum_list.append(df_slice)
#this WILL store AIC information, doesn't have it yet
resids_AICc = pd.concat(sum_list)
#return time_window index information
resids_AICc.index = resids.index
k = num_model_params[model_type]
#2 * number of parameters + number of data points
# https://www.sciencedirect.com/science/article/pii/S2468042719300508
resids_AICc["AICc"] = 2*k + resids_AICc["num_elements"]*np.log(resids_AICc["sum"]/resids_AICc["num_elements"]) + 2*k*(k+1)/(resids_AICc["num_elements"]-k-1)
resids_AICc.to_csv(make_dir(*filepath))
return resids_AICc
def boxplot_AICc_wrap(AICcs, pred_const_data, descriptions, brain_locs, models,
undersampled):
AICcs_lst = copy.deepcopy(AICcs)
dpi = 96
res = (1920,1080)
uniq_cols = pred_const_data.index.get_level_values("time_window").unique()
#reorder so I get index with patient category labels with the same order as AICs
midx = (pred_const_data.copy()
.unstack("patient_catg")
.reindex(AICcs[0].index)
.stack("patient_catg")
).index
for i in range(len(AICcs_lst)):
AICcs_lst[i].index = midx
data = pd.concat(AICcs_lst, keys = tuple(models)).unstack("time_window")
print(data)
temp = list(data.index.names)
temp[0] = "model"
data.index.names = temp
data = data.reorder_levels(["patient_catg","brain_loc","model","ID"])
s = "_undersampled" if undersampled else ""
dirpath = f"{os.getcwd()}/PvUgraphs&tables{s}/boxplots_AICc/"
models_in_order = pd.unique(data.index.get_level_values('model'))
# Convert only that level to a categorical with a frozen order
data.index = data.index.set_levels(
pd.CategoricalIndex(models_in_order, categories=models_in_order, ordered=True),
level='model'
)
for desc in descriptions:
for brain_loc in brain_locs:
fig, axs = plt.subplots(1, uniq_cols.size,
num = f'{desc} {brain_loc}',
sharey = "row",
figsize = (res[0]/dpi, res[1]/dpi),
dpi = dpi)
#the groupby inside boxplot sorts the data which is annoying
axs = data.loc[(desc, brain_loc)].boxplot(by = "model",
ax = axs,
rot = 45,
grid = False)
for j, ax in enumerate(axs):
plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")
ax.set_title(f"Window: {uniq_cols[j]}s")
ax.set(xlabel = "")
ax.grid(axis = "y")
for w, model in enumerate(models):
scat_data = data.loc[(desc, brain_loc, model), uniq_cols[j]]
x = np.random.normal(w+1, 0.05, len(scat_data))
ax.scatter(x = x,
y = scat_data,
alpha = 0.4)
axs[0].set_ylabel("AICc")
fig.suptitle(f"Akaike information criterion for small datasets for fits at different time windows in {brain_loc} for {desc}")
dirpath_d = f"{dirpath}{desc}/"
filepath = make_dir(dirpath_d) + f"AICc_model_comparison_{brain_loc}.png"
#only closes the figure when saving
#if not os.path.isfile(filepath):
print(f"Saving {filepath}")
plt.savefig(filepath)
plt.close(fig)
@filecheck_csv
def get_BIC(filepath, resids, model_type):
"""
Calculate the Akaike information criterion using ΔAIC = 2*k + n*ln(RSS)
where RSS = residual sum of squares, n is the number of samples and k
number of parameters of the model.
Parameters
----------
filepath : str
resids : pd.DataFrame
Array of residuals.
model_type : str
Returns
-------
resids_BIC : pd.DataFrame
"""
#model_type : k
num_model_params = {"Intravascular": 1,
"Patlak": 2,
"2CUptake": 3,
"ExTofts": 3,
"2CExchange": 4,
"2CUptakeDisp": 4,
"2CUptakePlug": 5}
sum_list = []
for cutoff in resids.index.get_level_values("time_window").unique():
#sum up the residuals corresponding to a window -> Series
window_data_slice = resids.loc[cutoff].pow(2).sum(axis = 1)
df_slice = pd.DataFrame(window_data_slice, columns = ["sum"])
#add the number of elements with index value greater than cutoff
df_slice["num_elements"] = (resids.columns.values.astype(float) >= int(cutoff)).sum()
sum_list.append(df_slice)
#this WILL store AIC information, doesn't have it yet
resids_BIC = pd.concat(sum_list)
#return time_window index information
resids_BIC.index = resids.index
k = num_model_params[model_type]
#2 * number of parameters + number of data points
# https://www.sciencedirect.com/science/article/pii/S2468042719300508
resids_BIC["BIC"] = resids_BIC["num_elements"]*np.log(resids_BIC["sum"]/resids_BIC["num_elements"]) + k*np.log(resids_BIC["num_elements"])
resids_BIC.to_csv(make_dir(*filepath))
return resids_BIC
def boxplot_BIC_wrap(BICs, pred_const_data, descriptions, brain_locs, models,
undersampled):
BICs_lst = copy.deepcopy(BICs)
dpi = 96
res = (1920,1080)
uniq_cols = pred_const_data.index.get_level_values("time_window").unique()
#reorder so I get index with patient category labels with the same order as AICs
midx = (pred_const_data.copy()
.unstack("patient_catg")
.reindex(BICs[0].index)
.stack("patient_catg")
).index
for i in range(len(BICs_lst)):
BICs_lst[i].index = midx
data = pd.concat(BICs_lst, keys = tuple(models)).unstack("time_window")
print(data)
temp = list(data.index.names)
temp[0] = "model"
data.index.names = temp
data = data.reorder_levels(["patient_catg","brain_loc","model","ID"])
s = "_undersampled" if undersampled else ""
dirpath = f"{os.getcwd()}/PvUgraphs&tables{s}/boxplots_BIC/"
models_in_order = pd.unique(data.index.get_level_values('model'))
# Convert only that level to a categorical with a frozen order
data.index = data.index.set_levels(
pd.CategoricalIndex(models_in_order, categories=models_in_order, ordered=True),
level='model'
)
for desc in descriptions:
for brain_loc in brain_locs:
fig, axs = plt.subplots(1, uniq_cols.size,
num = f'{desc} {brain_loc}',
sharey = "row",
figsize = (res[0]/dpi, res[1]/dpi),
dpi = dpi)
#the groupby inside boxplot sorts the data which is annoying
axs = data.loc[(desc, brain_loc)].boxplot(by = "model",
ax = axs,
rot = 45,
grid = False)
for j, ax in enumerate(axs):
plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")
ax.set_title(f"Window: {uniq_cols[j]}s")
ax.set(xlabel = "")
ax.grid(axis = "y")
for w, model in enumerate(models):
scat_data = data.loc[(desc, brain_loc, model), uniq_cols[j]]
x = np.random.normal(w+1, 0.05, len(scat_data))
ax.scatter(x = x,
y = scat_data,
alpha = 0.4)
axs[0].set_ylabel("BIC")
fig.suptitle(f"Bayesian information criterion for fits at different time windows in {brain_loc} for {desc}")
dirpath_d = f"{dirpath}{desc}/"
filepath = make_dir(dirpath_d) + f"BIC_model_comparison_{brain_loc}.png"
#only closes the figure when saving
#if not os.path.isfile(filepath):
print(f"Saving {filepath}")
plt.savefig(filepath)
plt.close(fig)

functions.py at commit d7989f3, no license · at the source

Overview

Authors: Martin Kozár^1,2, Ioana‐Emilia Mosneag^1,2, Sarah Al‐Bachari³, Igor Chernyavsky^4,5, Geoff J M Parker^6,7, Hervé Boutin^2,8,9, Laura M Parkes^2,10, Ingo Schiessl^2,9, Ben R Dickie^1,2,11

ORCID iDs: Martin Kozár, Ioana‐Emilia Mosneag, Geoff J M Parker, Hervé Boutin, Laura M Parkes, Ingo Schiessl, Ben R Dickie

Division of Imaging, Informatics and Data Sciences, School of Health Sciences, Faculty of Biology, Medicine and Health, University of Manchester, Manchester, UK
Geoffrey Jefferson Brain Research Centre, Faculty of Biology, Medicine and Health, University of Manchester, Manchester Academic Health Science Centre, Manchester, UK
Department of Clinical and Movement Neurosciences, University College London, London, UK
Department of Mathematics, School of Natural Sciences, Faculty of Science and Engineering, University of Manchester, Manchester, UK
Maternal and Fetal Health Research Centre, School of Medical Sciences, University of Manchester, Manchester, UK
Quantitative Imaging Group, UCL Hawkes Institute, Department of Medical Physics & Biomedical Engineering, University College London, London, UK
Bioxydyn Limited, Manchester, UK
Imaging Brain and Neuropsychiatry iBraiN 1253, Université de Tours, Inserm, Bat Planiol, UFR de Médecine, Tours, France
Division of Neuroscience, School of Biological Sciences, Faculty of Biology, Medicine and Health, Manchester Academic Health Science Centre, The University of Manchester, Manchester, UK
Division of Psychology, Communication and Human Neuroscience, School of Health Sciences, Faculty of Biology, Medicine and Health, University of Manchester, Manchester, UK
Division of Pharmacy and Optometry, School of Health Sciences, Faculty of Biology, Medicine and Health, University of Manchester, Manchester, UK

Institutions: Manchester Academic Health Science Centre (United Kingdom); University of Manchester (United Kingdom); University College London (United Kingdom); Université de Tours (France); Inserm (France); Imaging, Brain, and Neuropsychiatry (France)

Journal: Magnetic resonance in medicine, volume 96, issue 5, pages 2309-2320

Dates: received 23 July 2025; accepted 2 July 2026; published online 18 July 2026; in print November 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1002/mrm.70514 · PMID 42470230 · PMCID PMC13527286 · OpenAlex W7169658861

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: structural MRI / diffusion (modality), human (organism), mouse (organism), stroke (population), Parkinson's (population)

Methods: Connectivity, Statistics, fMRI & imaging

MeSH: Blood-Brain Barrier*, Brain*, Contrast Media*, Dynamic Contrast Enhanced Magnetic Resonance Imaging*, Magnetic Resonance Imaging*, Parkinson Disease*, Animals, Feasibility Studies, Female, Humans, Male, Metabolic Clearance Rate, Mice, Reproducibility of Results, Stroke (* major topic)

Topic: MRI in cancer diagnosis (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex

Funding: Medical Research Council DTP Studentship: MR/N013751/1; Sydney Driscoll Neuroscience Foundation and the Engineering and Physical Sciences Research Council (EP/M005909/1); Medical Research Council

Citations: cited by 1 paper (Europe PMC); 61 references in the paper

Abstract

Purpose: To determine the feasibility of measuring brain clearance of gadolinium contrast agent using standard intravenous DCE‐MRI acquisitions and evaluate the contribution of blood–brain barrier (BBB) and non‐BBB clearance routes.

Methods: Uptake and extended Tofts models were fit to DCE‐MRI data from people with Parkinson's disease, post‐stroke and controls. Models were compared using the Akaike information criterion. Key parameters (K trans, v p, v e, and k) were extracted from gray and white matter ROIs. In mice, two‐photon microscopy of intravenously injected Sulforhodamine 101 was used to confirm small tracer clearance kinetics without the confounds of partial volume effects or water exchange.

Results: The extended Tofts model provided a superior fit compared to all uptake‐only models. Extended Tofts estimates of the extravascular extracellular volume fraction v e were underestimated by a factor of 10–20 compared to known literature values (1.5% vs. 20%–30%). We hypothesized this bias may be due to an additional competing non‐BBB clearance mechanism not accounted for by the extended Tofts model. In both MRI and two‐photon data, we found that v e estimates trended toward literature values as BBB permeability (and thus BBB clearance) increased, supporting our hypothesis.

Conclusion: Modeling clearance of extravasated contrast agent from the brain improves fit quality compared to models that neglect clearance. Estimates of interstitial volume fraction were underestimated at low K trans but tended toward literature values as BBB permeability (and clearance) increased. This study indicates that estimating brain clearance of contrast agent using DCE‐MRI is feasible, and that estimates reflect a combination of BBB and non‐BBB clearance pathways.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

nimo-group/DCE-MRI-brain-clearance-2025

License: none: the authors keep all their rights

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: d7989f3228dd041cd6f599d1fbeba779535e513c, 15 December 2025

Languages: Python (5)

Size: 293 files, 5 scripts

Software Heritage: not archived

Found in: “Data Availability Statement”

Holds: README, environment (pyproject.toml, uv.lock)

Not found: license file, CITATION.cff, tests, continuous integration, documentation

Tools: Matplotlib (4 files), NumPy (4 files), pandas (4 files), SciPy (2 files), seaborn (2 files)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

6 files

functions.py, Python, 1,562 lines, 1 match
main.py, Python, 206 lines
multiprocess_function.py, Python, 31 lines
rebuttal_functions.py, Python, 1,325 lines, 1 match
rebuttal_main.py, Python, 215 lines
README.md, Text, 21 lines

At the source: github.com/nimo-group/DCE-MRI-brain-clearance-2025

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
5 scripts, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data Availability Statement

Raw data used within this study are available upon reasonable request from the corresponding author. Derived data and analysis scripts can be found on https://github.com/nimo‐group/DCE‐MRI‐brain‐clearance‐2025 (https://github.com/nimo-group/DCE-MRI-brain-clearance-2025).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

Publisher: n/a → Wiley

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 15 MeSH terms, 3 funders, 61 references.

Cite

This paper

Kozár, M., Mosneag, I., Al‐Bachari, S., Chernyavsky, I., Parker, G. J. M., Boutin, H., Parkes, L. M., Schiessl, I., & Dickie, B. R. (2026). Brain Clearance of Contrast Agent in Intravenous DCE-MRI Is Measurable and Cannot Be Explained by Clearance Across the BBB Alone. Magnetic resonance in medicine, 96(5), 2309-2320. https://doi.org/10.1002/mrm.70514

BibTeX

@article{kozar2026brain, author = {Kozár, Martin and Mosneag, Ioana‐Emilia and Al‐Bachari, Sarah and Chernyavsky, Igor and Parker, Geoff J M and Boutin, Hervé and Parkes, Laura M and Schiessl, Ingo and Dickie, Ben R}, title = {{Brain Clearance of Contrast Agent in Intravenous DCE-MRI Is Measurable and Cannot Be Explained by Clearance Across the BBB Alone}}, journal = {Magnetic resonance in medicine}, year = {2026}, month = jul, volume = {96}, number = {5}, pages = {2309--2320}, publisher = {Wiley}, issn = {0740-3194}, doi = {10.1002/mrm.70514}, url = {https://doi.org/10.1002/mrm.70514}, pmid = {42470230}, pmcid = {PMC13527286} }

RIS

TY - JOUR AU - Kozár, Martin AU - Mosneag, Ioana‐Emilia AU - Al‐Bachari, Sarah AU - Chernyavsky, Igor AU - Parker, Geoff J M AU - Boutin, Hervé AU - Parkes, Laura M AU - Schiessl, Ingo AU - Dickie, Ben R TI - Brain Clearance of Contrast Agent in Intravenous DCE-MRI Is Measurable and Cannot Be Explained by Clearance Across the BBB Alone T2 - Magnetic resonance in medicine J2 - Magn Reson Med PY - 2026 DA - 2026/07/18 VL - 96 IS - 5 SP - 2309 EP - 2320 SN - 0740-3194 PB - Wiley DO - 10.1002/mrm.70514 UR - https://doi.org/10.1002/mrm.70514 LA - en ER -

CSL-JSON

{ "id": "10.1002/mrm.70514", "type": "article-journal", "title": "Brain Clearance of Contrast Agent in Intravenous DCE-MRI Is Measurable and Cannot Be Explained by Clearance Across the BBB Alone", "container-title": "Magnetic resonance in medicine", "author": [ { "family": "Kozár", "given": "Martin" }, { "family": "Mosneag", "given": "Ioana‐Emilia" }, { "family": "Al‐Bachari", "given": "Sarah" }, { "family": "Chernyavsky", "given": "Igor" }, { "family": "Parker", "given": "Geoff J M" }, { "family": "Boutin", "given": "Hervé" }, { "family": "Parkes", "given": "Laura M" }, { "family": "Schiessl", "given": "Ingo" }, { "family": "Dickie", "given": "Ben R" } ], "container-title-short": "Magn Reson Med", "volume": "96", "issue": "5", "page": "2309-2320", "DOI": "10.1002/mrm.70514", "PMID": "42470230", "PMCID": "PMC13527286", "ISSN": "0740-3194", "publisher": "Wiley", "URL": "https://doi.org/10.1002/mrm.70514", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 18 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/nbm.70277 [code]: Hierarchical Bayesian Modelling Improves Microstructural Parameter Mapping in Diffusion and Exchange MRI Data.
Journal: NMR in biomedicine
In common: seaborn, SciPy, Matplotlib, 1 other tool, structural MRI / diffusion, 2 references, 2 authors
[2] doi:10.1073/pnas.2516601123 [code]: Unveiling the glymphatic system's role in brain aging: A comprehensive biomarker and modifiable intervention target.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: seaborn, pandas, SciPy, 2 other tools, structural MRI / diffusion, 5 references
[3] doi:10.1126/sciadv.aeb0404 [code]: MR-AIV reveals in vivo brain-wide fluid flow with physics-informed AI.
Journal: Science advances
In common: pandas, SciPy, Matplotlib, 1 other tool, structural MRI / diffusion, 4 references
[4] doi:10.1073/pnas.2526239123 [code]: Quantitative assessment of flow between cerebrospinal and interstitial fluid compartments in humans.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: structural MRI / diffusion, 5 references
[5] doi:10.1002/alz.71745 [code]: Regional astrocyte dysregulation and altered glymphatic-related markers in Alzheimer's disease frontal cortex.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: seaborn, pandas, SciPy, 2 other tools, 3 references
[6] doi:10.3389/fpsyt.2026.1848053: Neuroanatomical substrates of perivascular space index: a morphometric and tractography study.
Journal: Frontiers in psychiatry
In common: structural MRI / diffusion, 4 references
[7] doi:10.1038/s41467-026-76306-9 [code]: Non-invasive characterization of perivascular subarachnoid spaces.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, structural MRI / diffusion, 3 references
[8] doi:10.1177/0271678x261455452 [code]: Meningeal CSF transport varies across parasagittal dura subregions with age in humans.
Journal: Journal of cerebral blood flow and metabolism : official journal of the International Society of Cerebral Blood Flow and Metabolism
In common: structural MRI / diffusion, 4 references
[9] doi:10.1038/s41593-026-02358-1 [code]: Cerebral venous blood flow regulates intracerebral pressure and brain clearance via meningeal lymphatic vessels.
Journal: Nature neuroscience
In common: SciPy, Matplotlib, NumPy, stroke, structural MRI / diffusion, mouse, 2 references
[10] doi:10.1093/brain/awaf432: Global network and local vulnerabilities underlie brain atrophy across Parkinson's disease stages.
Journal: Brain : a journal of neurology
In common: Parkinson's, structural MRI / diffusion, author Laura M Parkes

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.