Skip to content

refactor: clean up DeltaM_model.py remove side effects, add docstrings and type hints #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Jitmisra wants to merge 1 commit into
base: main
Choose a base branch
from
Open

refactor: clean up DeltaM_model.py remove side effects, add docstrings and type hints #3

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
341 changes: 229 additions & 112 deletions src/bbb_exchange/DeltaM_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)