diff --git a/src/bbb_exchange/DeltaM_model.py b/src/bbb_exchange/DeltaM_model.py index 75b17b3..ef3ad8e 100644 --- a/src/bbb_exchange/DeltaM_model.py +++ b/src/bbb_exchange/DeltaM_model.py @@ -1,116 +1,233 @@ +""" +Single-TE ASL Kinetic Models (Chappell et al.) + +Implements the tissue and arterial compartment models for Arterial Spin Labelling +(ASL) signal modelling based on: + + Chappell, M. A., et al. (2010). "Separation of macrovascular signal in + multi-inversion time arterial spin labelling MRI." + Magnetic Resonance in Medicine, 63(5), 1357-1365. + https://doi.org/10.1002/mrm.22320 + +Equation [1]: Tissue compartment signal (dm_tiss) +Equation [2]: Arterial compartment signal (dm_art) + +Functions: + dm_tiss -- Tissue compartment DeltaM signal + dm_art -- Arterial compartment DeltaM signal + DeltaM_model_ext -- Combined tissue + arterial model + plot_deltaM_curves -- Visualization of model components +""" + import numpy as np -import matplotlib.pyplot as plt - - -def dm_tiss(t, Dt, tau, f, M0a, a, T1, T1a, k=0.9): - """ - Chappell equation [1] (tissue component) and [2] (arterial component) from the paper 'Separation of macrovascular signal in multi‐inversion time arterial spin labelling MRI', https://doi.org/10.1002/mrm.22320 - Modelling of the ASL signal (DeltaM) over time. - Inputs: - t - Time points (seconds) - att - Arterial transit time (Dt in paper), time until the blood arrives in the voxel [s] - cbf - Cerebral blood flow (CBF) in [ml/min/100g] - m0a - Scaling factor: M0 value of the arterial blood signal (e.g. from M0. nii) - tau - Duration of the labelling phase (tau) [s] - T1a - Longitudinal relaxation time of the arterial blood [s] (default: 1.65s) - lambd - Blood-tissue-water partition coefficient (default: 0.9) - alpha - Labelling efficiency (e.g. 0.85 * 0.8) - Output: - deltaM - Expected signal curve DeltaM(t) for a single voxel +from typing import Union + + +def dm_tiss( + t: np.ndarray, + Dt: float, + tau: float, + f: float, + M0a: float, + a: float, + T1: float, + T1a: float, + k: float = 0.9, +) -> np.ndarray: + """ + Tissue compartment of the ASL signal (Chappell Eq. [1]). + + Computes the expected DeltaM signal arising from labelled water that has + exchanged into the tissue compartment. + + Parameters + ---------- + t : array_like + Time points in seconds (post-labelling delays). + Dt : float + Arterial transit time (ATT) — time for blood to arrive in the voxel [s]. + tau : float + Duration of the labelling phase [s]. + f : float + Perfusion rate (flow) [ml/g/s]. Note: CBF [ml/min/100g] = f * 6000. + M0a : float + Equilibrium magnetisation of arterial blood (scaling factor). + a : float + Labelling efficiency (e.g. 0.85 * 0.8 = 0.68). + T1 : float + Longitudinal relaxation time of tissue [s]. + T1a : float + Longitudinal relaxation time of arterial blood [s]. + k : float, optional + Blood-tissue water partition coefficient (lambda). Default is 0.9. + + Returns + ------- + DM : np.ndarray + Expected tissue DeltaM signal at each time point. """ - t = np.asarray(t, dtype=float) - T1app = 1.0 / T1 + f / k - R = 1.0 / T1app - 1.0 / T1a - - DM = np.zeros_like(t) - - # case 2: Dt <= t <= Dt + tau - mask2 = (t >= Dt) & (t <= Dt + tau) - if np.any(mask2): - tt = t[mask2] - term = (np.exp(R * tt) - np.exp(R * Dt)) - DM[mask2] = (2 * a * M0a * f * np.exp(-tt / T1app) / R) * term - - # case 3: t > Dt + tau - mask3 = t > Dt + tau - if np.any(mask3): - tt = t[mask3] - term = (np.exp(R * (Dt + tau)) - np.exp(R * Dt)) - DM[mask3] = (2 * a * M0a * f * np.exp(-tt / T1app) / R) * term - - return DM - - -def dm_art(t, Dta, ta, aBV, M0a, a, T1a): - """ - Chappell equation [2] (arterial component) from the paper 'Separation of macrovascular signal in multi‐inversion time arterial spin labelling MRI', https://doi.org/10.1002/mrm.22320 - Modelling of the ASL signal (DeltaM) over time. - Inputs: - t - Time points (seconds) - att - Arterial transit time (Dt in paper), time until the blood arrives in the voxel [s] - cbf - Cerebral blood flow (CBF) in [ml/min/100g] - m0a - Scaling factor: M0 value of the arterial blood signal (e.g. from M0. nii) - tau - Duration of the labelling phase (tau) [s] - T1a - Longitudinal relaxation time of the arterial blood [s] (default: 1.65s) - Dta - arterial arrival time (s) - ta - artieral bolus time (s) - abV - arterial blood volume - lambd - Blood-tissue-water partition coefficient (default: 0.9) - alpha - Labelling efficiency (e.g. 0.85 * 0.8) - Output: - deltaM - Expected signal curve DeltaM(t) for a single voxel + t = np.asarray(t, dtype=float) + T1app = 1.0 / T1 + f / k + R = 1.0 / T1app - 1.0 / T1a + + DM = np.zeros_like(t) + + # Case 2: Dt <= t <= Dt + tau (during labelled bolus arrival) + mask2 = (t >= Dt) & (t <= Dt + tau) + if np.any(mask2): + tt = t[mask2] + term = np.exp(R * tt) - np.exp(R * Dt) + DM[mask2] = (2 * a * M0a * f * np.exp(-tt / T1app) / R) * term + + # Case 3: t > Dt + tau (after bolus has fully arrived) + mask3 = t > Dt + tau + if np.any(mask3): + tt = t[mask3] + term = np.exp(R * (Dt + tau)) - np.exp(R * Dt) + DM[mask3] = (2 * a * M0a * f * np.exp(-tt / T1app) / R) * term + + return DM + + +def dm_art( + t: np.ndarray, + Dta: float, + ta: float, + aBV: float, + M0a: float, + a: float, + T1a: float, +) -> np.ndarray: """ - t = np.asarray(t, dtype=float) - DM = np.zeros_like(t) - # Dta <= t <= Dta + ta, only if signal exists - mask = (t >= Dta) & (t <= Dta + ta) - if np.any(mask): - tt = t[mask] - DM[mask] = 2 * a * M0a * aBV * np.exp(-tt / T1a) - return DM - - -def DeltaM_model_ext(t, params): - - return dm_tiss(t, params['Dt'], params['tau'], params['f'], params['M0a'], - params['a'], params['T1'], params['T1a'], params.get('k', 0.9)) + \ - dm_art(t, params['Dta'], params['ta'], params['aBV'], params['M0a'], - params['a'], params['T1a']) - - -t = np.linspace(0, 3.0, 301) - -params = { - 'f': 0.01, # ml/g/s - 'Dt': 0.7, # s - 'tau': 1.0, # s - 'M0a': 1.0, # scaling - 'a': 0.95, # efficiency factor - 'T1': 1.3, # s - 'T1a': 1.6, # s - 'k': 0.9, # lambda - 'aBV': 0.01, # arterial blood volume fraction - 'Dta': 0.5, # arterial arrival time - 'ta': 1.0 # arteral bolus time -} + Arterial compartment of the ASL signal (Chappell Eq. [2]). + + Computes the expected DeltaM signal arising from labelled water still + within the arterial vasculature (macrovascular component). + + Parameters + ---------- + t : array_like + Time points in seconds. + Dta : float + Arterial arrival time — transit time to the arterial compartment [s]. + ta : float + Arterial bolus duration [s]. + aBV : float + Arterial blood volume fraction (dimensionless). + M0a : float + Equilibrium magnetisation of arterial blood (scaling factor). + a : float + Labelling efficiency. + T1a : float + Longitudinal relaxation time of arterial blood [s]. + + Returns + ------- + DM : np.ndarray + Expected arterial DeltaM signal at each time point. + """ + t = np.asarray(t, dtype=float) + DM = np.zeros_like(t) -""" -# Visualising the model -Dm_tissue = dm_tiss(t, params['Dt'], params['tau'], params['f'], params['M0a'], - params['a'], params['T1'], params['T1a'], params['k']) -Dm_art = dm_art(t, params['Dta'], params['ta'], params['aBV'], params['M0a'], - params['a'], params['T1a']) -Dm_tot = Dm_tissue + Dm_art - - -plt.figure(figsize=(8, 4)) -plt.plot(t, Dm_tot, label='DM_total (tissue + arterial)') -plt.plot(t, Dm_tissue, '--', label='tissue component') -plt.plot(t, Dm_art, ':', label='arterial component') -plt.xlabel('time (s)') -plt.ylabel('DeltaM') -plt.title('Simulated ASL kinetic curves (Equation 1 + 2)') -plt.legend() -plt.tight_layout() -plt.show() -""" + # Arterial signal present only during Dta <= t <= Dta + ta + mask = (t >= Dta) & (t <= Dta + ta) + if np.any(mask): + tt = t[mask] + DM[mask] = 2 * a * M0a * aBV * np.exp(-tt / T1a) + + return DM + + +def DeltaM_model_ext(t: np.ndarray, params: dict) -> np.ndarray: + """ + Extended ASL model combining tissue and arterial compartments. + + This is a convenience wrapper that sums the tissue (Eq. [1]) and arterial + (Eq. [2]) components of the Chappell model. + + Parameters + ---------- + t : array_like + Time points in seconds. + params : dict + Dictionary of model parameters. Required keys: + - 'Dt' : Arterial transit time [s] + - 'tau' : Labelling duration [s] + - 'f' : Perfusion rate [ml/g/s] + - 'M0a' : Equilibrium magnetisation of arterial blood + - 'a' : Labelling efficiency + - 'T1' : T1 of tissue [s] + - 'T1a' : T1 of arterial blood [s] + - 'Dta' : Arterial arrival time [s] + - 'ta' : Arterial bolus duration [s] + - 'aBV' : Arterial blood volume fraction + Optional keys: + - 'k' : Partition coefficient (default: 0.9) + + Returns + ------- + DM : np.ndarray + Combined tissue + arterial DeltaM signal. + """ + return dm_tiss( + t, params['Dt'], params['tau'], params['f'], params['M0a'], + params['a'], params['T1'], params['T1a'], params.get('k', 0.9) + ) + dm_art( + t, params['Dta'], params['ta'], params['aBV'], params['M0a'], + params['a'], params['T1a'] + ) + + +def plot_deltaM_curves(t: np.ndarray, params: dict) -> None: + """ + Plot the tissue, arterial, and combined DeltaM signal curves. + + Parameters + ---------- + t : array_like + Time points in seconds. + params : dict + Model parameters (same format as DeltaM_model_ext). + """ + import matplotlib.pyplot as plt + Dm_tissue = dm_tiss( + t, params['Dt'], params['tau'], params['f'], params['M0a'], + params['a'], params['T1'], params['T1a'], params.get('k', 0.9) + ) + Dm_arterial = dm_art( + t, params['Dta'], params['ta'], params['aBV'], params['M0a'], + params['a'], params['T1a'] + ) + Dm_total = Dm_tissue + Dm_arterial + + plt.figure(figsize=(8, 4)) + plt.plot(t, Dm_total, label='DM_total (tissue + arterial)') + plt.plot(t, Dm_tissue, '--', label='Tissue component') + plt.plot(t, Dm_arterial, ':', label='Arterial component') + plt.xlabel('Time (s)') + plt.ylabel('DeltaM') + plt.title('Simulated ASL kinetic curves (Chappell Eq. 1 + 2)') + plt.legend() + plt.tight_layout() + plt.show() + + +if __name__ == "__main__": + # Example parameters for visualising the model + t = np.linspace(0, 3.0, 301) + + params = { + 'f': 0.01, # Perfusion rate [ml/g/s] + 'Dt': 0.7, # Arterial transit time [s] + 'tau': 1.0, # Labelling duration [s] + 'M0a': 1.0, # Equilibrium magnetisation (scaling) + 'a': 0.95, # Labelling efficiency + 'T1': 1.3, # T1 tissue [s] + 'T1a': 1.6, # T1 blood [s] + 'k': 0.9, # Partition coefficient (lambda) + 'aBV': 0.01, # Arterial blood volume fraction + 'Dta': 0.5, # Arterial arrival time [s] + 'ta': 1.0, # Arterial bolus duration [s] + } + + plot_deltaM_curves(t, params)