From 1639ce57698101d69d592b84e58827e2b7566067 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Thu, 29 Aug 2024 10:25:29 +0300 Subject: [PATCH 01/13] Initial implementation for 2CXM --- .pre-commit-config.yaml | 2 +- src/osipi/__init__.py | 3 +- src/osipi/_tissue.py | 118 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 121 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index c2429e4..3d66742 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -12,7 +12,7 @@ repos: args: - --exclude=venv,.git,__pycache__,.eggs,.mypy_cache,.pytest_cache - --max-line-length=100 - - --ignore=E203 + - --ignore=E203, W503 - --per-file-ignores=__init__.py:F401 - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.6.0 diff --git a/src/osipi/__init__.py b/src/osipi/__init__.py index 454bb7e..06fdc38 100644 --- a/src/osipi/__init__.py +++ b/src/osipi/__init__.py @@ -7,7 +7,8 @@ from ._tissue import ( tofts, - extended_tofts + extended_tofts, + two_cxm ) from ._signal import ( diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 3e9e13d..7e43bfa 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -334,3 +334,121 @@ def extended_tofts( ct = ct_func(t) return ct + + +def two_cxm( + t: np.ndarray, + ca: np.ndarray, + ps: float, + fp: float, + ve: float, + vp: float, + Ta: float = 30.0, + discretization_method: str = "conv", +) -> np.ndarray: + """The 2CX model allows bi-directional exchange of indicator between vascular and + extra vascular extracellular compartments. Indicator is assumed to be + well mixed within each compartment. + + Args: + t (np.ndarray): array of time points in units of sec.[OSIPI code Q.GE1.004] + ca (np.ndarray): Arterial concentrations in mM for each + time point in t.[OSIPI code Q.IC1.001] + ps (float): Permeability surface area product in units of 1/min. [OSIPI code Q.PH1.009] + fp (float): Plasma flow in units of mL/min/100g. [OSIPI code Q.PH1.010] + ve (float): Relative volume fraction of the + extracellular extravascular compartment (e). [OSIPI code Q.PH1.001.[e]] + vp (float): Relative volyme fraction of the plasma + compartment (p). [OSIPI code Q.PH1.001.[p]] + Ta (float, optional): Arterial delay time, i.e., difference in onset time between + tissue curve and AIF in units of sec. Defaults to 30 seconds. [OSIPI code Q.PH1.007] + discretization_method (str, optional): Defines the discretization method. Options include + – 'conv': Numerical convolution (default) [OSIPI code G.DI1.001] + – 'exp': Exponential convolution [OSIPI code G.DI1.006] + """ + + if not np.allclose(np.diff(t), np.diff(t)[0]): + warnings.warn( + ("Non-uniform time spacing detected. Time array may be resampled."), stacklevel=2 + ) + + K_plus = ( + 1 + / 2 + * ( + (fp + ps) / vp + + ps / vp + + np.sqrt((fp + ps) / vp + ps / vp) ** 2 + - 4 * fp * ps / (vp * ve) + ) + ) + K_minus = ( + 1 + / 2 + * ( + (fp + ps) / vp + + ps / vp + - np.sqrt((fp + ps) / vp + ps / vp) ** 2 + - 4 * fp * ps / (vp * ve) + ) + ) + E_ = (K_plus + fp / vp) / (K_plus + K_minus) + + exp_k_plus = np.exp(-K_plus * t) + exp_k_minus = np.exp(-K_minus * t) + + imp = fp * exp_k_plus + E_ * (exp_k_plus - exp_k_minus) + + # Shift the AIF by the arterial delay time (if not zero) + if Ta != 0: + f = interp1d( + t, + ca, + kind="linear", + bounds_error=False, + fill_value=0, + ) + ca = (t > Ta) * f(t - Ta) + + if np.allclose(np.diff(t), np.diff(t)[0]): + # Convolve impulse response with AIF + convolution = np.convolve(ca, imp) * t[1] + + ct = fp * convolution + + else: + # Resample at the smallest spacing + dt = np.min(np.diff(t)) + t_resampled = np.linspace(t[0], t[-1], int((t[-1] - t[0]) / dt)) + ca_func = interp1d( + t, + ca, + kind="quadratic", + bounds_error=False, + fill_value=0, + ) + imp_func = interp1d( + t, + imp, + kind="quadratic", + bounds_error=False, + fill_value=0, + ) + ca_resampled = ca_func(t_resampled) + imp_resampled = imp_func(t_resampled) + # Convolve impulse response with AIF + convolution = np.convolve(ca_resampled, imp_resampled) * t_resampled[1] + + # Discard unwanted points and make sure time spacing is correct + ct_resampled = fp * convolution[0 : len(t_resampled)] + # Restore time grid spacing + ct_func = interp1d( + t_resampled, + ct_resampled, + kind="quadratic", + bounds_error=False, + fill_value=0, + ) + ct = ct_func(t) + + return ct From 5216dadf6169bed81565095278a165ffb67f5827 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sun, 1 Sep 2024 14:25:26 +0300 Subject: [PATCH 02/13] Add test for 2cxm --- src/osipi/_tissue.py | 51 +++++++++++++++++--------------------------- tests/test_tissue.py | 7 ++++++ 2 files changed, 26 insertions(+), 32 deletions(-) diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 7e43bfa..c6c133f 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -339,10 +339,10 @@ def extended_tofts( def two_cxm( t: np.ndarray, ca: np.ndarray, - ps: float, - fp: float, - ve: float, - vp: float, + E: float, + Fp: float, + Ve: float, + Vp: float, Ta: float = 30.0, discretization_method: str = "conv", ) -> np.ndarray: @@ -354,11 +354,11 @@ def two_cxm( t (np.ndarray): array of time points in units of sec.[OSIPI code Q.GE1.004] ca (np.ndarray): Arterial concentrations in mM for each time point in t.[OSIPI code Q.IC1.001] - ps (float): Permeability surface area product in units of 1/min. [OSIPI code Q.PH1.009] - fp (float): Plasma flow in units of mL/min/100g. [OSIPI code Q.PH1.010] - ve (float): Relative volume fraction of the + E (float): Extraction fraction in units of mL/min/100g. + Fp (float): Plasma flow in units of mL/min/100g. [OSIPI code Q.PH1.003] + Ve (float): Relative volume fraction of the extracellular extravascular compartment (e). [OSIPI code Q.PH1.001.[e]] - vp (float): Relative volyme fraction of the plasma + Vp (float): Relative volyme fraction of the plasma compartment (p). [OSIPI code Q.PH1.001.[p]] Ta (float, optional): Arterial delay time, i.e., difference in onset time between tissue curve and AIF in units of sec. Defaults to 30 seconds. [OSIPI code Q.PH1.007] @@ -372,32 +372,19 @@ def two_cxm( ("Non-uniform time spacing detected. Time array may be resampled."), stacklevel=2 ) - K_plus = ( - 1 - / 2 - * ( - (fp + ps) / vp - + ps / vp - + np.sqrt((fp + ps) / vp + ps / vp) ** 2 - - 4 * fp * ps / (vp * ve) - ) - ) - K_minus = ( - 1 - / 2 - * ( - (fp + ps) / vp - + ps / vp - - np.sqrt((fp + ps) / vp + ps / vp) ** 2 - - 4 * fp * ps / (vp * ve) - ) - ) - E_ = (K_plus + fp / vp) / (K_plus + K_minus) + Tp = Vp / Fp * (1 - E) + Te = Ve * (1 - E) / (Fp * E) + Tb = Vp / Fp + + K_plus = 0.5 * ((1 / Tp) + (1 / Te) + np.sqrt(((1 / Tp) + (1 / Te)) ** 2 - (4 / (Te * Tb)))) + K_minus = 0.5 * ((1 / Tp) + (1 / Te) - np.sqrt(((1 / Tp) + (1 / Te)) ** 2 - (4 / (Te * Tb)))) + + A = (K_plus - (1 / Tb)) / (K_plus - K_minus) exp_k_plus = np.exp(-K_plus * t) exp_k_minus = np.exp(-K_minus * t) - imp = fp * exp_k_plus + E_ * (exp_k_plus - exp_k_minus) + imp = exp_k_plus + A * (exp_k_minus - exp_k_plus) # Shift the AIF by the arterial delay time (if not zero) if Ta != 0: @@ -414,7 +401,7 @@ def two_cxm( # Convolve impulse response with AIF convolution = np.convolve(ca, imp) * t[1] - ct = fp * convolution + ct = Fp * convolution else: # Resample at the smallest spacing @@ -440,7 +427,7 @@ def two_cxm( convolution = np.convolve(ca_resampled, imp_resampled) * t_resampled[1] # Discard unwanted points and make sure time spacing is correct - ct_resampled = fp * convolution[0 : len(t_resampled)] + ct_resampled = Fp * convolution[0 : len(t_resampled)] # Restore time grid spacing ct_func = interp1d( t_resampled, diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 304cba9..2dd669a 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -124,6 +124,13 @@ def test_tissue_extended_tofts(): assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) +def test_tissue_2compartment_model(): + t = np.arange(0, 6 * 60, 1, dtype=float) + ca = osipi.aif_parker(t) + ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) + assert np.round(np.max(ct)) < np.round(np.max(ca)) + + if __name__ == "__main__": test_tissue_tofts() test_tissue_extended_tofts() From 1a9107022c3f4997d544eb76393d2b7a3b548039 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sun, 1 Sep 2024 15:43:09 +0300 Subject: [PATCH 03/13] Add some tests --- src/osipi/_tissue.py | 2 +- tests/test_tissue.py | 20 ++++++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index c6c133f..6f76c87 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -401,7 +401,7 @@ def two_cxm( # Convolve impulse response with AIF convolution = np.convolve(ca, imp) * t[1] - ct = Fp * convolution + ct = Fp * convolution[0 : len(t)] else: # Resample at the smallest spacing diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 2dd669a..ceb4d66 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -125,11 +125,31 @@ def test_tissue_extended_tofts(): def test_tissue_2compartment_model(): + # 1. Basic operation of the function - test that the peak tissue t = np.arange(0, 6 * 60, 1, dtype=float) ca = osipi.aif_parker(t) ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) assert np.round(np.max(ct)) < np.round(np.max(ca)) + # 2. Basic operation of the function - test with non-uniform spacing of + t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 + ca = osipi.aif_parker(t) + ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) + assert np.round(np.max(ct)) < np.round(np.max(ca)) + + # 3. The offset option - test that the tissue concentration is shifted + t = np.arange(0, 6 * 60, 1) + ca = osipi.aif_parker(t) + ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3, Ta=60.0) + assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 + + # 4. Test that the area under the ct curve is approximately the sum of + # the extracellular volume and the plasma volume + t = np.arange(0, 6 * 60, 1) + ca = osipi.aif_parker(t) + ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) + assert math.isclose(np.trapz(ct, t)/np.trapz(ca, t), 0.2 + 0.3, abs_tol=1e-1) + if __name__ == "__main__": test_tissue_tofts() From 6dd676326d0abf429ee01459a21725a9dd543c96 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sun, 1 Sep 2024 17:31:57 +0300 Subject: [PATCH 04/13] Added 2CXM documentation and example --- docs/examples/tissue/plot_two_cxm.py | 38 ++++++++++++++++++ docs/references/models/tissue_models/2cxm.md | 16 ++++++++ docs/references/models/tissue_models/index.md | 1 + mkdocs.yml | 1 + src/osipi/_tissue.py | 39 +++++++++++++++++++ tests/test_tissue.py | 2 +- 6 files changed, 96 insertions(+), 1 deletion(-) create mode 100644 docs/examples/tissue/plot_two_cxm.py create mode 100644 docs/references/models/tissue_models/2cxm.md diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py new file mode 100644 index 0000000..e6b7d2a --- /dev/null +++ b/docs/examples/tissue/plot_two_cxm.py @@ -0,0 +1,38 @@ +""" +============= +The Two Compartment Exchange Model +============= + +Simulating tissue concentrations from two compartment models with different settings. +""" + +import matplotlib.pyplot as plt + +# %% +# Import necessary packages +import numpy as np +import osipi + +# %% +# Generate Parker AIF with default settings. + +# Define time points in units of seconds - in this case we use a time +# resolution of 1 sec and a total duration of 6 minutes. +t = np.arange(0, 6 * 60, 1) + +# Create an AIF with default settings +ca = osipi.aif_parker(t) + +# %% +# Plot the tissue concentrations for an extracellular volume fraction +# of 0.2, plasma volume fraction of 0.3, extraction fraction of 0.15 +# and flow rate of 0.2 ml/min +E = 0.15 # Extraction fraction +Fp = 0.2 # Flow rate in ml/min +Ve = 0.2 # Extracellular volume fraction +Vp = 0.3 # Plasma volume fraction +ct = osipi.two_cxm(t, ca, E=E, Fp=Fp, Ve=Ve, Vp=Vp) +plt.plot(t, ct, "b-", label=f"E = {E}, Fp = {Fp}, Ve = {Ve}, Vp = {Vp}") +plt.xlabel("Time (sec)") +plt.ylabel("Tissue concentration (mM)") +plt.show() diff --git a/docs/references/models/tissue_models/2cxm.md b/docs/references/models/tissue_models/2cxm.md new file mode 100644 index 0000000..38bc95f --- /dev/null +++ b/docs/references/models/tissue_models/2cxm.md @@ -0,0 +1,16 @@ +# osipi.Two Compartment Exchange Model + +::: osipi.two_cxm + +## Example using `osipi.two_cxm` + +
+ The Two Compartment Exchange Model +
+ + + The Two Compartment Exchange Model + + +
+
diff --git a/docs/references/models/tissue_models/index.md b/docs/references/models/tissue_models/index.md index 298b575..5584f43 100644 --- a/docs/references/models/tissue_models/index.md +++ b/docs/references/models/tissue_models/index.md @@ -1,3 +1,4 @@ - [tofts](tofts.md) - [extended_tofts](extended_tofts.md) +- [two_cxm](2cxm.md) diff --git a/mkdocs.yml b/mkdocs.yml index a742a1c..183e016 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -75,6 +75,7 @@ nav: - references/models/tissue_models/index.md - osipi.tofts: references/models/tissue_models/tofts.md - osipi.extended_tofts: references/models/tissue_models/extended_tofts.md + - osipi.two_compartment_exchange: references/models/tissue_models/2cxm.md - Examples: generated/gallery diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 41561b5..220fac4 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -363,6 +363,45 @@ def two_cxm( discretization_method (str, optional): Defines the discretization method. Options include – 'conv': Numerical convolution (default) [OSIPI code G.DI1.001] – 'exp': Exponential convolution [OSIPI code G.DI1.006] + + Returns: + np.ndarray: Tissue concentrations in mM for each time point in t. + + See Also: + `tofts`, `extended_tofts` + + References: + - Lexicon url: + https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/#2CXM + - Lexicon code: M.IC1.009 + - OSIPI name: Two Compartment Model + - Adapted from contributions by: LEK_UoEdinburgh_UK, ST_USyd_AUS, MJT_UoEdinburgh_UK + + Example: + + Create an array of time points covering 6 min in steps of 1 sec, + calculate the Parker AIF at these time points, calculate tissue concentrations + using the 2CX model and plot the results. + + Import packages: + + >>> import matplotlib.pyplot as plt + >>> import osipi + + Calculate AIF: + + >>> t = np.arange(0, 6 * 60, 1) + >>> ca = osipi.aif_parker(t) + + Calculate tissue concentrations and plot: + + >>> E = 0.15 # in units of mL/min/100g + >>> Fp = 0.2 # in units of mL/min/100g + >>> Ve = 0.2 # takes values from 0 to 1 + >>> Vp = 0.3 # takes values from 0 to 1 + >>> ct = osipi.two_cxm(t, ca, E, Fp, Ve, Vp) + >>> plt.plot(t, ca, "r", t, ct, "b") + """ if not np.allclose(np.diff(t), np.diff(t)[0]): diff --git a/tests/test_tissue.py b/tests/test_tissue.py index ceb4d66..09f8a98 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -148,7 +148,7 @@ def test_tissue_2compartment_model(): t = np.arange(0, 6 * 60, 1) ca = osipi.aif_parker(t) ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) - assert math.isclose(np.trapz(ct, t)/np.trapz(ca, t), 0.2 + 0.3, abs_tol=1e-1) + assert math.isclose(np.trapz(ct, t) / np.trapz(ca, t), 0.2 + 0.3, abs_tol=1e-1) if __name__ == "__main__": From 2188bf59f76189469096d19da3db2b0d14eab42a Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sun, 15 Sep 2024 19:11:22 +0300 Subject: [PATCH 05/13] Added use case and corrected units --- docs/examples/tissue/plot_two_cxm.py | 14 ++++++---- docs/references/models/tissue_models/2cxm.md | 3 ++- src/osipi/__init__.py | 2 +- src/osipi/_tissue.py | 27 +++++++++++++------- tests/test_tissue.py | 26 +++++++++++-------- 5 files changed, 46 insertions(+), 26 deletions(-) diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py index e6b7d2a..d5572d2 100644 --- a/docs/examples/tissue/plot_two_cxm.py +++ b/docs/examples/tissue/plot_two_cxm.py @@ -18,7 +18,7 @@ # Define time points in units of seconds - in this case we use a time # resolution of 1 sec and a total duration of 6 minutes. -t = np.arange(0, 6 * 60, 1) +t = np.arange(0, 6 * 60, 1, dtype=float) # Create an AIF with default settings ca = osipi.aif_parker(t) @@ -27,12 +27,16 @@ # Plot the tissue concentrations for an extracellular volume fraction # of 0.2, plasma volume fraction of 0.3, extraction fraction of 0.15 # and flow rate of 0.2 ml/min -E = 0.15 # Extraction fraction -Fp = 0.2 # Flow rate in ml/min +E = 0.3 # Extraction fraction +Fp = 0.1 # Flow rate in ml/min Ve = 0.2 # Extracellular volume fraction -Vp = 0.3 # Plasma volume fraction -ct = osipi.two_cxm(t, ca, E=E, Fp=Fp, Ve=Ve, Vp=Vp) +Vp = 0.1 # Plasma volume fraction +ct = osipi.two_compartment_exchange_model(t, ca, E=E, Fp=Fp, Ve=Ve, Vp=Vp) plt.plot(t, ct, "b-", label=f"E = {E}, Fp = {Fp}, Ve = {Ve}, Vp = {Vp}") +t = np.arange(0, 6 * 60, 1, dtype=float) +ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.2) +plt.plot(t, ct, "r-", label="E = 0.15, Fp = 0.2, Ve = 0.2, Vp = 2") plt.xlabel("Time (sec)") plt.ylabel("Tissue concentration (mM)") +plt.legend() plt.show() diff --git a/docs/references/models/tissue_models/2cxm.md b/docs/references/models/tissue_models/2cxm.md index 38bc95f..0b3d636 100644 --- a/docs/references/models/tissue_models/2cxm.md +++ b/docs/references/models/tissue_models/2cxm.md @@ -1,6 +1,7 @@ # osipi.Two Compartment Exchange Model -::: osipi.two_cxm +::: osipi.two_compartment_exchange_model + ## Example using `osipi.two_cxm` diff --git a/src/osipi/__init__.py b/src/osipi/__init__.py index 06fdc38..ffe955b 100644 --- a/src/osipi/__init__.py +++ b/src/osipi/__init__.py @@ -8,7 +8,7 @@ from ._tissue import ( tofts, extended_tofts, - two_cxm + two_compartment_exchange_model ) from ._signal import ( diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 220fac4..8441e3e 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -334,7 +334,7 @@ def extended_tofts( return ct -def two_cxm( +def two_compartment_exchange_model( t: np.ndarray, ca: np.ndarray, E: float, @@ -342,9 +342,8 @@ def two_cxm( Ve: float, Vp: float, Ta: float = 30.0, - discretization_method: str = "conv", ) -> np.ndarray: - """The 2CX model allows bi-directional exchange of indicator between vascular and + """The 2 compartment exchange model allows bi-directional exchange of indicator between vascular and extra vascular extracellular compartments. Indicator is assumed to be well mixed within each compartment. @@ -352,8 +351,8 @@ def two_cxm( t (np.ndarray): array of time points in units of sec.[OSIPI code Q.GE1.004] ca (np.ndarray): Arterial concentrations in mM for each time point in t.[OSIPI code Q.IC1.001] - E (float): Extraction fraction in units of mL/min/100g. - Fp (float): Plasma flow in units of mL/min/100g. [OSIPI code Q.PH1.003] + E (float): Extraction fraction in units of mL/min/100mL. + Fp (float): Plasma flow in units of mL/min/100mL. [OSIPI code Q.PH1.003] Ve (float): Relative volume fraction of the extracellular extravascular compartment (e). [OSIPI code Q.PH1.001.[e]] Vp (float): Relative volyme fraction of the plasma @@ -395,15 +394,25 @@ def two_cxm( Calculate tissue concentrations and plot: - >>> E = 0.15 # in units of mL/min/100g - >>> Fp = 0.2 # in units of mL/min/100g - >>> Ve = 0.2 # takes values from 0 to 1 + >>> E = 0.15 # in units of mL/min/100mL + >>> Fp = 5 # in units of mL/min/100mL + >>> Ve = 0.1 # takes values from 0 to 1 >>> Vp = 0.3 # takes values from 0 to 1 - >>> ct = osipi.two_cxm(t, ca, E, Fp, Ve, Vp) + >>> ct = osipi.two_compartment_exchange_model(t, ca, E, Fp, Ve, Vp) >>> plt.plot(t, ca, "r", t, ct, "b") """ + if Vp == 0: + Ktrans = E * Fp + + ct = tofts(t, ca, Ktrans, Ve) + return ct + + # Convert units + t /= 60.0 + Ta /= 60.0 + if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( ("Non-uniform time spacing detected. Time array may be resampled."), stacklevel=2 diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 09f8a98..e7f9537 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -124,35 +124,41 @@ def test_tissue_extended_tofts(): assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) -def test_tissue_2compartment_model(): +def test_tissue_two_compartment_exchange_model(): # 1. Basic operation of the function - test that the peak tissue + # concentration is less than the peak AIF t = np.arange(0, 6 * 60, 1, dtype=float) ca = osipi.aif_parker(t) - ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) + ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) assert np.round(np.max(ct)) < np.round(np.max(ca)) # 2. Basic operation of the function - test with non-uniform spacing of + # time array t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 ca = osipi.aif_parker(t) - ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) + ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) assert np.round(np.max(ct)) < np.round(np.max(ca)) # 3. The offset option - test that the tissue concentration is shifted - t = np.arange(0, 6 * 60, 1) + # from the AIF by the specified offset time + t = np.arange(0, 6 * 60, 1, dtype=float) ca = osipi.aif_parker(t) - ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3, Ta=60.0) + ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.1, Ta=60.0) assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 - # 4. Test that the area under the ct curve is approximately the sum of - # the extracellular volume and the plasma volume - t = np.arange(0, 6 * 60, 1) + # 4. Test specific use cases + + # Test when Vp is 0 the output matches tofts model + t = np.arange(0, 6 * 60, 1, dtype=float) ca = osipi.aif_parker(t) - ct = osipi.two_cxm(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) - assert math.isclose(np.trapz(ct, t) / np.trapz(ca, t), 0.2 + 0.3, abs_tol=1e-1) + ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0) + ct_tofts = osipi.tofts(t, ca, Ktrans=0.03, ve=0.2) + assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) if __name__ == "__main__": test_tissue_tofts() test_tissue_extended_tofts() + test_tissue_two_compartment_exchange_model() print("All tissue concentration model tests passed!!") From 2ae41e2b0d37ee7bed041a0428a2195a5345f248 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Mon, 7 Oct 2024 08:20:24 +0300 Subject: [PATCH 06/13] Added initial implementation for two compartment exchange model --- docs/examples/tissue/plot_two_cxm.py | 15 +-- src/osipi/_tissue.py | 174 +++++---------------------- 2 files changed, 37 insertions(+), 152 deletions(-) diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py index d5572d2..5724755 100644 --- a/docs/examples/tissue/plot_two_cxm.py +++ b/docs/examples/tissue/plot_two_cxm.py @@ -27,15 +27,12 @@ # Plot the tissue concentrations for an extracellular volume fraction # of 0.2, plasma volume fraction of 0.3, extraction fraction of 0.15 # and flow rate of 0.2 ml/min -E = 0.3 # Extraction fraction -Fp = 0.1 # Flow rate in ml/min -Ve = 0.2 # Extracellular volume fraction -Vp = 0.1 # Plasma volume fraction -ct = osipi.two_compartment_exchange_model(t, ca, E=E, Fp=Fp, Ve=Ve, Vp=Vp) -plt.plot(t, ct, "b-", label=f"E = {E}, Fp = {Fp}, Ve = {Ve}, Vp = {Vp}") -t = np.arange(0, 6 * 60, 1, dtype=float) -ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.2) -plt.plot(t, ct, "r-", label="E = 0.15, Fp = 0.2, Ve = 0.2, Vp = 2") +Ps = 0.15 # Extraction fraction +Fp = 20 # Flow rate in ml/min +Ve = 0.1 # Extracellular volume fraction +Vp = 0.02 # Plasma volume fraction +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp, Ps=Ps, Ve=Ve, Vp=Vp) +plt.plot(t, ct, "b-", label=f" Fp = {Fp},Ps = {Ps}, Ve = {Ve}, Vp = {Vp}") plt.xlabel("Time (sec)") plt.ylabel("Tissue concentration (mM)") plt.legend() diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 8441e3e..430cd1b 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -2,6 +2,7 @@ import numpy as np from scipy.interpolate import interp1d +from scipy.signal import convolve from ._convolution import exp_conv @@ -337,151 +338,38 @@ def extended_tofts( def two_compartment_exchange_model( t: np.ndarray, ca: np.ndarray, - E: float, Fp: float, + Ps: float, Ve: float, Vp: float, Ta: float = 30.0, ) -> np.ndarray: - """The 2 compartment exchange model allows bi-directional exchange of indicator between vascular and - extra vascular extracellular compartments. Indicator is assumed to be - well mixed within each compartment. - - Args: - t (np.ndarray): array of time points in units of sec.[OSIPI code Q.GE1.004] - ca (np.ndarray): Arterial concentrations in mM for each - time point in t.[OSIPI code Q.IC1.001] - E (float): Extraction fraction in units of mL/min/100mL. - Fp (float): Plasma flow in units of mL/min/100mL. [OSIPI code Q.PH1.003] - Ve (float): Relative volume fraction of the - extracellular extravascular compartment (e). [OSIPI code Q.PH1.001.[e]] - Vp (float): Relative volyme fraction of the plasma - compartment (p). [OSIPI code Q.PH1.001.[p]] - Ta (float, optional): Arterial delay time, i.e., difference in onset time between - tissue curve and AIF in units of sec. Defaults to 30 seconds. [OSIPI code Q.PH1.007] - discretization_method (str, optional): Defines the discretization method. Options include - – 'conv': Numerical convolution (default) [OSIPI code G.DI1.001] - – 'exp': Exponential convolution [OSIPI code G.DI1.006] - - Returns: - np.ndarray: Tissue concentrations in mM for each time point in t. - - See Also: - `tofts`, `extended_tofts` - - References: - - Lexicon url: - https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/#2CXM - - Lexicon code: M.IC1.009 - - OSIPI name: Two Compartment Model - - Adapted from contributions by: LEK_UoEdinburgh_UK, ST_USyd_AUS, MJT_UoEdinburgh_UK - - Example: - - Create an array of time points covering 6 min in steps of 1 sec, - calculate the Parker AIF at these time points, calculate tissue concentrations - using the 2CX model and plot the results. - - Import packages: - - >>> import matplotlib.pyplot as plt - >>> import osipi - - Calculate AIF: - - >>> t = np.arange(0, 6 * 60, 1) - >>> ca = osipi.aif_parker(t) - - Calculate tissue concentrations and plot: - - >>> E = 0.15 # in units of mL/min/100mL - >>> Fp = 5 # in units of mL/min/100mL - >>> Ve = 0.1 # takes values from 0 to 1 - >>> Vp = 0.3 # takes values from 0 to 1 - >>> ct = osipi.two_compartment_exchange_model(t, ca, E, Fp, Ve, Vp) - >>> plt.plot(t, ca, "r", t, ct, "b") - - """ - - if Vp == 0: - Ktrans = E * Fp - - ct = tofts(t, ca, Ktrans, Ve) - return ct - - # Convert units - t /= 60.0 - Ta /= 60.0 - - if not np.allclose(np.diff(t), np.diff(t)[0]): - warnings.warn( - ("Non-uniform time spacing detected. Time array may be resampled."), stacklevel=2 - ) - - Tp = Vp / Fp * (1 - E) - Te = Ve * (1 - E) / (Fp * E) - Tb = Vp / Fp - - K_plus = 0.5 * ((1 / Tp) + (1 / Te) + np.sqrt(((1 / Tp) + (1 / Te)) ** 2 - (4 / (Te * Tb)))) - K_minus = 0.5 * ((1 / Tp) + (1 / Te) - np.sqrt(((1 / Tp) + (1 / Te)) ** 2 - (4 / (Te * Tb)))) - - A = (K_plus - (1 / Tb)) / (K_plus - K_minus) - - exp_k_plus = np.exp(-K_plus * t) - exp_k_minus = np.exp(-K_minus * t) - - imp = exp_k_plus + A * (exp_k_minus - exp_k_plus) - - # Shift the AIF by the arterial delay time (if not zero) - if Ta != 0: - f = interp1d( - t, - ca, - kind="linear", - bounds_error=False, - fill_value=0, - ) - ca = (t > Ta) * f(t - Ta) - - if np.allclose(np.diff(t), np.diff(t)[0]): - # Convolve impulse response with AIF - convolution = np.convolve(ca, imp) * t[1] - - ct = Fp * convolution[0 : len(t)] - - else: - # Resample at the smallest spacing - dt = np.min(np.diff(t)) - t_resampled = np.linspace(t[0], t[-1], int((t[-1] - t[0]) / dt)) - ca_func = interp1d( - t, - ca, - kind="quadratic", - bounds_error=False, - fill_value=0, - ) - imp_func = interp1d( - t, - imp, - kind="quadratic", - bounds_error=False, - fill_value=0, - ) - ca_resampled = ca_func(t_resampled) - imp_resampled = imp_func(t_resampled) - # Convolve impulse response with AIF - convolution = np.convolve(ca_resampled, imp_resampled) * t_resampled[1] - - # Discard unwanted points and make sure time spacing is correct - ct_resampled = Fp * convolution[0 : len(t_resampled)] - # Restore time grid spacing - ct_func = interp1d( - t_resampled, - ct_resampled, - kind="quadratic", - bounds_error=False, - fill_value=0, - ) - ct = ct_func(t) - - return ct + fp_per_s = Fp / (60.0 * 100.0) + ps_per_s = Ps / 60.0 + v = Ve + Vp + T = v / fp_per_s + tc = Vp / fp_per_s + te = Ve / ps_per_s + + sig_p = ((T + te) + np.sqrt((T + te) ** 2 - 4 * tc * te)) / (2 * tc * te) + sig_n = ((T + tc) - np.sqrt((T + tc) ** 2 - 4 * tc * te)) / (2 * tc * te) + + irf_cp = ( + Vp + * sig_p + * sig_n + * ((1 - te * sig_n) * np.exp(-t * sig_n) + (te * sig_p - 1.0) * np.exp(-t * sig_p)) + / (sig_p - sig_n) + ) + + irf_ce = Ve * sig_p * sig_n * (np.exp(-t * sig_n) - np.exp(-t * sig_p)) / (sig_p - sig_n) + irf_cp[[0]] /= 2 + irf_ce[[0]] /= 2 + + dt = np.min(np.diff(t)) + + Cp = dt * convolve(ca, irf_cp, mode="full", method="auto")[: len(t)] + Ce = dt * convolve(ca, irf_ce, mode="full", method="auto")[: len(t)] + + Ct = Cp + Ce + return Ct From f40ae1acb96e3f6ed6d504e65b541722a836b00d Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Mon, 7 Oct 2024 14:41:16 +0300 Subject: [PATCH 07/13] Added TCXM adapted from MJT implementation --- docs/examples/tissue/plot_two_cxm.py | 8 ++-- docs/references/models/tissue_models/2cxm.md | 2 +- src/osipi/_tissue.py | 41 ++++++++++++++++++++ tests/test_tissue.py | 30 +------------- 4 files changed, 47 insertions(+), 34 deletions(-) diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py index 5724755..fb47451 100644 --- a/docs/examples/tissue/plot_two_cxm.py +++ b/docs/examples/tissue/plot_two_cxm.py @@ -27,10 +27,10 @@ # Plot the tissue concentrations for an extracellular volume fraction # of 0.2, plasma volume fraction of 0.3, extraction fraction of 0.15 # and flow rate of 0.2 ml/min -Ps = 0.15 # Extraction fraction -Fp = 20 # Flow rate in ml/min -Ve = 0.1 # Extracellular volume fraction -Vp = 0.02 # Plasma volume fraction +Ps = 0.05 # Permeability surface area product in ml/min +Fp = 10 # Flow rate in ml/min +Ve = 0.2 # Extracellular volume fraction +Vp = 0.01 # Plasma volume fraction ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp, Ps=Ps, Ve=Ve, Vp=Vp) plt.plot(t, ct, "b-", label=f" Fp = {Fp},Ps = {Ps}, Ve = {Ve}, Vp = {Vp}") plt.xlabel("Time (sec)") diff --git a/docs/references/models/tissue_models/2cxm.md b/docs/references/models/tissue_models/2cxm.md index 0b3d636..68f2f7f 100644 --- a/docs/references/models/tissue_models/2cxm.md +++ b/docs/references/models/tissue_models/2cxm.md @@ -3,7 +3,7 @@ ::: osipi.two_compartment_exchange_model -## Example using `osipi.two_cxm` +## Example using `osipi.two_compartment_exchange_model`
The Two Compartment Exchange Model diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 430cd1b..c134c41 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -344,9 +344,44 @@ def two_compartment_exchange_model( Vp: float, Ta: float = 30.0, ) -> np.ndarray: + """Two compartment exchange model + + Tracer flows from the AIF to the blood plasma compartment; two-way leakage + between the plasma and extracellular compartments(EES) is permitted. + + Args: + t: 1D array of times(s). + ca: Arterial concentrations in mM for each time point in t. + Fp: Blood plasma flow rate into a unit tissue volume in ml/min. + Ps: Permeability surface area product in ml/min. + Ve: Extracellular volume fraction. + Vp: Plasma volume fraction. + Ta: Arterial delay time, i.e., difference in onset time between tissue curve and AIF in units of sec. + + Returns: + Ct: Tissue concentrations in mM + + References: + - Lexicon url: https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/#indicator-kinetic-models + - Lexicon code: M.IC1.009 + - OSIPI name: Two Compartment Exchange Model + - Adapted from contributions by: MJT_UoEdinburgh_UK + + """ + if not np.allclose(np.diff(t), np.diff(t)[0]): + warnings.warn( + ("Non-uniform time spacing detected. Time array may be" " resampled."), + stacklevel=2, + ) + + # Convert units fp_per_s = Fp / (60.0 * 100.0) ps_per_s = Ps / 60.0 + + # Calculate the impulse response function v = Ve + Vp + + # Mean transit time T = v / fp_per_s tc = Vp / fp_per_s te = Ve / ps_per_s @@ -354,6 +389,8 @@ def two_compartment_exchange_model( sig_p = ((T + te) + np.sqrt((T + te) ** 2 - 4 * tc * te)) / (2 * tc * te) sig_n = ((T + tc) - np.sqrt((T + tc) ** 2 - 4 * tc * te)) / (2 * tc * te) + # Calculate the impulse response function for the plasma compartment and EES + irf_cp = ( Vp * sig_p @@ -363,13 +400,17 @@ def two_compartment_exchange_model( ) irf_ce = Ve * sig_p * sig_n * (np.exp(-t * sig_n) - np.exp(-t * sig_p)) / (sig_p - sig_n) + irf_cp[[0]] /= 2 irf_ce[[0]] /= 2 dt = np.min(np.diff(t)) + # get concentration in plasma and EES + Cp = dt * convolve(ca, irf_cp, mode="full", method="auto")[: len(t)] Ce = dt * convolve(ca, irf_ce, mode="full", method="auto")[: len(t)] + # get tissue concentration Ct = Cp + Ce return Ct diff --git a/tests/test_tissue.py b/tests/test_tissue.py index e7f9537..8fbf2b5 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -125,35 +125,7 @@ def test_tissue_extended_tofts(): def test_tissue_two_compartment_exchange_model(): - # 1. Basic operation of the function - test that the peak tissue - # concentration is less than the peak AIF - t = np.arange(0, 6 * 60, 1, dtype=float) - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - - # 2. Basic operation of the function - test with non-uniform spacing of - # time array - t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.3) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - - # 3. The offset option - test that the tissue concentration is shifted - # from the AIF by the specified offset time - t = np.arange(0, 6 * 60, 1, dtype=float) - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0.1, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 - - # 4. Test specific use cases - - # Test when Vp is 0 the output matches tofts model - t = np.arange(0, 6 * 60, 1, dtype=float) - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, E=0.15, Fp=0.2, Ve=0.2, Vp=0) - ct_tofts = osipi.tofts(t, ca, Ktrans=0.03, ve=0.2) - assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) + pass if __name__ == "__main__": From 50cd6656122ffb6422c5c6acfb0d8013571cf047 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Fri, 15 Nov 2024 23:34:19 +0200 Subject: [PATCH 08/13] Corrected some errors Fixed some errors in the implementation, some units and naming --- docs/examples/tissue/plot_two_cxm.py | 32 ++++++++---- src/osipi/_tissue.py | 77 +++++++++++++++++++++------- 2 files changed, 80 insertions(+), 29 deletions(-) diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py index fb47451..bdf4549 100644 --- a/docs/examples/tissue/plot_two_cxm.py +++ b/docs/examples/tissue/plot_two_cxm.py @@ -18,21 +18,33 @@ # Define time points in units of seconds - in this case we use a time # resolution of 1 sec and a total duration of 6 minutes. -t = np.arange(0, 6 * 60, 1, dtype=float) +t = np.arange(0, 5 * 60, 1, dtype=float) # Create an AIF with default settings ca = osipi.aif_parker(t) # %% -# Plot the tissue concentrations for an extracellular volume fraction -# of 0.2, plasma volume fraction of 0.3, extraction fraction of 0.15 -# and flow rate of 0.2 ml/min -Ps = 0.05 # Permeability surface area product in ml/min -Fp = 10 # Flow rate in ml/min -Ve = 0.2 # Extracellular volume fraction -Vp = 0.01 # Plasma volume fraction -ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp, Ps=Ps, Ve=Ve, Vp=Vp) -plt.plot(t, ct, "b-", label=f" Fp = {Fp},Ps = {Ps}, Ve = {Ve}, Vp = {Vp}") +#Plot the tissue concentrations for an extracellular volume fraction +#of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min +#and flow rate of 10 ml/min +PS = 15 # Permeability surface area product in ml/min +Fp = [10, 25] # Flow rate in ml/min +ve = 0.1 # Extracellular volume fraction +vp = [0.1,0.02] # Plasma volume fraction + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[0], PS=PS, ve=ve, vp=vp[0]) +plt.plot(t, ct, "b-", label=f" Fp = {Fp[0]},PS = {PS}, ve = {ve}, vp = {vp[0]}") + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[1], PS=PS, ve=ve, vp=vp[0]) +plt.plot(t, ct, "r-", label=f" Fp = {Fp[1]},PS = {PS}, ve = {ve}, vp = {vp[0]}") + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[0], PS=PS, ve=ve, vp=vp[1]) +plt.plot(t, ct, "g-", label=f" Fp = {Fp[0]},PS = {PS}, ve = {ve}, vp = {vp[1]}") + +ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[1], PS=PS, ve=ve, vp=vp[1]) +plt.plot(t, ct, "y-", label=f" Fp = {Fp[1]},PS = {PS}, ve = {ve}, vp = {vp[1]}") + + plt.xlabel("Time (sec)") plt.ylabel("Tissue concentration (mM)") plt.legend() diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index c134c41..379e612 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -339,9 +339,9 @@ def two_compartment_exchange_model( t: np.ndarray, ca: np.ndarray, Fp: float, - Ps: float, - Ve: float, - Vp: float, + PS: float, + ve: float, + vp: float, Ta: float = 30.0, ) -> np.ndarray: """Two compartment exchange model @@ -350,23 +350,51 @@ def two_compartment_exchange_model( between the plasma and extracellular compartments(EES) is permitted. Args: - t: 1D array of times(s). - ca: Arterial concentrations in mM for each time point in t. - Fp: Blood plasma flow rate into a unit tissue volume in ml/min. - Ps: Permeability surface area product in ml/min. - Ve: Extracellular volume fraction. - Vp: Plasma volume fraction. - Ta: Arterial delay time, i.e., difference in onset time between tissue curve and AIF in units of sec. + t: 1D array of times(s). [OSIPI code is Q.GE1.004] + ca: Arterial concentrations in mM for each time point in t. [OSIPI code is Q.IC1.001.[a]] + Fp: Blood plasma flow rate into a unit tissue volume in ml/min. [OSIPI code is Q.PH1.002] + PS: Permeability surface area product in ml/min. [OSIPI code is Q.PH1.004] + ve: Extracellular volume fraction. [OSIPI code Q.PH1.001.[e]] + vp: Plasma volume fraction. [OSIPI code Q.PH1.001.[p]] + Ta: Arterial delay time, i.e., + difference in onset time between tissue curve and AIF in units of sec. [OSIPI code Q.PH1.007] Returns: Ct: Tissue concentrations in mM + See Also: + `extended_tofts` + References: - Lexicon url: https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/#indicator-kinetic-models - Lexicon code: M.IC1.009 - OSIPI name: Two Compartment Exchange Model - Adapted from contributions by: MJT_UoEdinburgh_UK + Example: + Create an array of time points covering 6 min in steps of 1 sec, + calculate the Parker AIF at these time points, calculate tissue concentrations + using the Two Compartment Exchange model and plot the results. + + >>> import matplotlib.pyplot as plt + >>> import osipi + + Calculate AIF + + >>> t = np.arange(0, 6 * 60, 0.1) + >>> ca = osipi.aif_parker(t) + + Plot the tissue concentrations for an extracellular volume fraction + of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min + and flow rate of 10 ml/min + + >>> PS = 5 # Permeability surface area product in ml/min + >>> Fp = 10 # Flow rate in ml/min + >>> ve = 0.2 # Extracellular volume fraction + >>> vp = 0.1 # Plasma volume fraction + >>> ct = osipi.two_compartment_exchange_model(t, ca, Fp, PS, ve, vp) + >>> plt.plot(t, ca, 'r', t, ct, 'g') + """ if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( @@ -376,41 +404,52 @@ def two_compartment_exchange_model( # Convert units fp_per_s = Fp / (60.0 * 100.0) - ps_per_s = Ps / 60.0 + ps_per_s = PS / (60.0 * 100.0) # Calculate the impulse response function - v = Ve + Vp + v = ve + vp # Mean transit time T = v / fp_per_s - tc = Vp / fp_per_s - te = Ve / ps_per_s + tc = vp / fp_per_s + te = ve / ps_per_s + + upsample_factor = 1 + n = t.size + n_upsample = (n - 1) * upsample_factor + 1 + t_upsample = np.linspace(t[0], t[-1], n_upsample) + tau_upsample = t_upsample - t[0] sig_p = ((T + te) + np.sqrt((T + te) ** 2 - 4 * tc * te)) / (2 * tc * te) - sig_n = ((T + tc) - np.sqrt((T + tc) ** 2 - 4 * tc * te)) / (2 * tc * te) + sig_n = ((T + te) - np.sqrt((T + te) ** 2 - 4 * tc * te)) / (2 * tc * te) # Calculate the impulse response function for the plasma compartment and EES irf_cp = ( - Vp + vp * sig_p * sig_n - * ((1 - te * sig_n) * np.exp(-t * sig_n) + (te * sig_p - 1.0) * np.exp(-t * sig_p)) + * ((1 - te * sig_n) * np.exp(-tau_upsample * sig_n) + (te * sig_p - 1.0) * np.exp(-tau_upsample * sig_p)) / (sig_p - sig_n) ) - irf_ce = Ve * sig_p * sig_n * (np.exp(-t * sig_n) - np.exp(-t * sig_p)) / (sig_p - sig_n) + irf_ce = ve * sig_p * sig_n * (np.exp(-tau_upsample * sig_n) - np.exp(-tau_upsample * sig_p)) / (sig_p - sig_n) irf_cp[[0]] /= 2 irf_ce[[0]] /= 2 - dt = np.min(np.diff(t)) + dt = np.min(np.diff(t)) / upsample_factor # get concentration in plasma and EES Cp = dt * convolve(ca, irf_cp, mode="full", method="auto")[: len(t)] Ce = dt * convolve(ca, irf_ce, mode="full", method="auto")[: len(t)] + t_upsample = np.linspace(t[0], t[-1], n_upsample) + + Cp = np.interp(t, t_upsample, Cp) + Ce = np.interp(t,t_upsample, Ce) + # get tissue concentration Ct = Cp + Ce return Ct From 3535067b1312904b14e491aa1644a0aa069a4c3b Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Fri, 15 Nov 2024 23:35:49 +0200 Subject: [PATCH 09/13] some fixes --- docs/examples/tissue/plot_two_cxm.py | 8 ++++---- src/osipi/_tissue.py | 17 +++++++++++++---- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py index bdf4549..8e608a3 100644 --- a/docs/examples/tissue/plot_two_cxm.py +++ b/docs/examples/tissue/plot_two_cxm.py @@ -24,13 +24,13 @@ ca = osipi.aif_parker(t) # %% -#Plot the tissue concentrations for an extracellular volume fraction -#of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min -#and flow rate of 10 ml/min +# Plot the tissue concentrations for an extracellular volume fraction +# of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min +# and flow rate of 10 ml/min PS = 15 # Permeability surface area product in ml/min Fp = [10, 25] # Flow rate in ml/min ve = 0.1 # Extracellular volume fraction -vp = [0.1,0.02] # Plasma volume fraction +vp = [0.1, 0.02] # Plasma volume fraction ct = osipi.two_compartment_exchange_model(t, ca, Fp=Fp[0], PS=PS, ve=ve, vp=vp[0]) plt.plot(t, ct, "b-", label=f" Fp = {Fp[0]},PS = {PS}, ve = {ve}, vp = {vp[0]}") diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 379e612..55af8fd 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -393,7 +393,7 @@ def two_compartment_exchange_model( >>> ve = 0.2 # Extracellular volume fraction >>> vp = 0.1 # Plasma volume fraction >>> ct = osipi.two_compartment_exchange_model(t, ca, Fp, PS, ve, vp) - >>> plt.plot(t, ca, 'r', t, ct, 'g') + >>> plt.plot(t, ca, "r", t, ct, "g") """ if not np.allclose(np.diff(t), np.diff(t)[0]): @@ -429,11 +429,20 @@ def two_compartment_exchange_model( vp * sig_p * sig_n - * ((1 - te * sig_n) * np.exp(-tau_upsample * sig_n) + (te * sig_p - 1.0) * np.exp(-tau_upsample * sig_p)) + * ( + (1 - te * sig_n) * np.exp(-tau_upsample * sig_n) + + (te * sig_p - 1.0) * np.exp(-tau_upsample * sig_p) + ) / (sig_p - sig_n) ) - irf_ce = ve * sig_p * sig_n * (np.exp(-tau_upsample * sig_n) - np.exp(-tau_upsample * sig_p)) / (sig_p - sig_n) + irf_ce = ( + ve + * sig_p + * sig_n + * (np.exp(-tau_upsample * sig_n) - np.exp(-tau_upsample * sig_p)) + / (sig_p - sig_n) + ) irf_cp[[0]] /= 2 irf_ce[[0]] /= 2 @@ -448,7 +457,7 @@ def two_compartment_exchange_model( t_upsample = np.linspace(t[0], t[-1], n_upsample) Cp = np.interp(t, t_upsample, Cp) - Ce = np.interp(t,t_upsample, Ce) + Ce = np.interp(t, t_upsample, Ce) # get tissue concentration Ct = Cp + Ce From c703ed69682edb3d388c9d21a432d89afc043cd0 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Wed, 27 Nov 2024 20:32:13 +0200 Subject: [PATCH 10/13] Add some tests --- src/osipi/_tissue.py | 16 +++++++++++++++- tests/test_tissue.py | 28 +++++++++++++++++++++++++++- 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 55af8fd..dd85f08 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -396,6 +396,11 @@ def two_compartment_exchange_model( >>> plt.plot(t, ca, "r", t, ct, "g") """ + if vp == 0: + E = 1 - np.exp(-PS / Fp) + Ktrans = E * Fp + return tofts(t, ca, Ktrans, ve, Ta, discretization_method="conv") + if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( ("Non-uniform time spacing detected. Time array may be" " resampled."), @@ -449,8 +454,17 @@ def two_compartment_exchange_model( dt = np.min(np.diff(t)) / upsample_factor - # get concentration in plasma and EES + if Ta != 0: + f = interp1d( + t, + ca, + kind="linear", + bounds_error=False, + fill_value=0, + ) + ca = (t > Ta) * f(t - Ta) + # get concentration in plasma and EES Cp = dt * convolve(ca, irf_cp, mode="full", method="auto")[: len(t)] Ce = dt * convolve(ca, irf_ce, mode="full", method="auto")[: len(t)] diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 8fbf2b5..143d0be 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -125,7 +125,33 @@ def test_tissue_extended_tofts(): def test_tissue_two_compartment_exchange_model(): - pass + # 1. Basic operation of the function - test that the peak tissue + # concentration is less than the peak AIF + t = np.linspace(0, 6 * 60, 360) + ca = osipi.aif_parker(t) + ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3) + assert np.round(np.max(ct)) < np.round(np.max(ca)) + + # 2. Basic operation of the function - test with non-uniform spacing of + # time array + t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 + ca = osipi.aif_parker(t) + ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3) + assert np.round(np.max(ct)) < np.round(np.max(ca)) + + # 3. The offset option - test that the tissue concentration is shifted + # from the AIF by the specified offset time + t = np.arange(0, 6 * 60, 1) + ca = osipi.aif_parker(t) + ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3, Ta=60.0) + assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 + + # 4. Handle case where VP is 0, it should behave like the tofts model + t = np.arange(0, 6 * 60, 1) + ca = osipi.aif_parker(t) + ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0) + ct_tofts = osipi.tofts(t, ca, Ktrans=3.93, ve=0.2) + assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) if __name__ == "__main__": From 9ca630a2e576922f9965fffbe5c7af4ff80b8ac2 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Wed, 19 Feb 2025 06:33:44 +0200 Subject: [PATCH 11/13] refactor(tests): implement pytest fixtures and parameterization --- tests/test_tissue.py | 152 ++++++++++++++++++++++--------------------- 1 file changed, 78 insertions(+), 74 deletions(-) diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 143d0be..36764ac 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -2,20 +2,31 @@ import numpy as np import osipi +import pytest -def test_tissue_tofts(): - """1. +@pytest.fixture +def time_array(): + return np.arange(0, 6 * 60, 1) - Basic operation of the function - test that the peak tissue concentration is less than the peak + +@pytest.fixture +def aif(time_array): + return osipi.aif_parker(time_array) + + +@pytest.fixture +def base_params(): + return {"Ktrans": 0.6, "ve": 0.2} + + +def test_tissue_tofts(time_array, aif, base_params): + """1. Basic operation of the function - test that the peak tissue concentration is less than the peak AIF """ - - t = np.linspace(0, 6 * 60, 360) - ca = osipi.aif_parker(t) - ct = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - assert np.round(np.max(ct)) < np.round(np.max(ca)) + ct = osipi.tofts(time_array, aif, discretization_method="exp", **base_params) + assert np.round(np.max(ct)) < np.round(np.max(aif)) # 2. Basic operation of the function - test with non-uniform spacing of # time array @@ -26,45 +37,42 @@ def test_tissue_tofts(): # 3. The offset option - test that the tissue concentration is shifted # from the AIF by the specified offset time - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 + ct = osipi.tofts(time_array, aif, Ta=60.0, **base_params) + assert (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 == 60.0 # 4. Test that the discretization options give almost the same result - # time step must be very small t = np.arange(0, 6 * 60, 0.01) ca = osipi.aif_parker(t) - ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, discretization_method="exp") + ct_conv = osipi.tofts(t, ca, **base_params) + ct_exp = osipi.tofts(t, ca, discretization_method="exp", **base_params) assert np.allclose(ct_conv, ct_exp, rtol=1e-4, atol=1e-3) # 5. Test that the ratio of the area under the ct and ca curves is # approximately the extracellular volume - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, discretization_method="exp") - assert math.isclose(np.trapz(ct_conv, t) / np.trapz(ca, t), 0.2, abs_tol=1e-1) - assert math.isclose(np.trapz(ct_exp, t) / np.trapz(ca, t), 0.2, abs_tol=1e-1) + + ct_conv = osipi.tofts(time_array, aif, **base_params) + ct_exp = osipi.tofts(time_array, aif, discretization_method="exp", **base_params) + assert math.isclose( + np.trapz(ct_conv, time_array) / np.trapz(aif, time_array), 0.2, abs_tol=1e-1 + ) + assert math.isclose(np.trapz(ct_exp, time_array) / np.trapz(aif, time_array), 0.2, abs_tol=1e-1) # 6. Test specific use cases - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.tofts(t, ca, Ktrans=0, ve=0.2) + ct_conv = osipi.tofts(time_array, aif, Ktrans=0, ve=0.2) assert np.count_nonzero(ct_conv) == 0 - ct_exp = osipi.tofts(t, ca, Ktrans=0, ve=0.2, discretization_method="exp") + ct_exp = osipi.tofts(time_array, aif, Ktrans=0, ve=0.2, discretization_method="exp") assert np.count_nonzero(ct_exp) == 0 - ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0) + ct_conv = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0) assert np.count_nonzero(ct_conv) == 0 - ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0, discretization_method="exp") + ct_exp = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0, discretization_method="exp") assert np.count_nonzero(ct_exp) == 0 -def test_tissue_extended_tofts(): +def test_tissue_extended_tofts(time_array, aif): # 1. Basic operation of the function - test that the peak tissue # concentration is less than the peak AIF t = np.linspace(0, 6 * 60, 360) @@ -81,10 +89,8 @@ def test_tissue_extended_tofts(): # 3. The offset option - test that the tissue concentration is shifted # from the AIF by the specified offset time - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 + ct = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, Ta=60.0) + assert (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 == 60.0 # 4. Test that the discretization options give almost the same result - # time step must be very small @@ -96,67 +102,65 @@ def test_tissue_extended_tofts(): # 5. Test that the ratio of the area under the ct and ca curves is # approximately the extracellular volume plus the plasma volume - - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) - ct_exp = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method="exp") + ct_conv = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3) + ct_exp = osipi.extended_tofts( + time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method="exp" + ) assert math.isclose( - np.trapz(ct_conv, t) / np.trapz(ca, t), + np.trapz(ct_conv, time_array) / np.trapz(aif, time_array), 0.2 + 0.3, abs_tol=1e-1, ) - assert math.isclose(np.trapz(ct_exp, t) / np.trapz(ca, t), 0.2 + 0.3, abs_tol=1e-1) + assert math.isclose( + np.trapz(ct_exp, time_array) / np.trapz(aif, time_array), 0.2 + 0.3, abs_tol=1e-1 + ) # 6. Test specific use cases - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct_conv = osipi.extended_tofts(t, ca, Ktrans=0, ve=0.2, vp=0.3) - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) - ct_exp = osipi.extended_tofts(t, ca, Ktrans=0, ve=0.2, vp=0.3, discretization_method="exp") - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) + ct_conv = osipi.extended_tofts(time_array, aif, Ktrans=0, ve=0.2, vp=0.3) + assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) - ct_conv = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0, vp=0.3) - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) + ct_exp = osipi.extended_tofts( + time_array, aif, Ktrans=0, ve=0.2, vp=0.3, discretization_method="exp" + ) + assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) - ct_exp = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0, vp=0.3, discretization_method="exp") - assert np.allclose(ct_conv, ca * 0.3, rtol=1e-4, atol=1e-3) + ct_conv = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0, vp=0.3) + assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) + ct_exp = osipi.extended_tofts( + time_array, aif, Ktrans=0.6, ve=0, vp=0.3, discretization_method="exp" + ) + assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) -def test_tissue_two_compartment_exchange_model(): - # 1. Basic operation of the function - test that the peak tissue - # concentration is less than the peak AIF - t = np.linspace(0, 6 * 60, 360) - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 2. Basic operation of the function - test with non-uniform spacing of - # time array - t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 +@pytest.mark.parametrize( + "time_points", + [ + (np.linspace(0, 6 * 60, 360), "uniform"), + (np.geomspace(1, 6 * 60 + 1, num=360) - 1, "non-uniform"), + ], +) +def test_2cxm_basic_operation(time_points): + """Test that the peak tissue concentration is less than the peak AIF for different time arrays""" + t, _ = time_points ca = osipi.aif_parker(t) ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3) assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 3. The offset option - test that the tissue concentration is shifted - # from the AIF by the specified offset time - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0.3, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(ca > 0.0)) - 1) * 1 == 60.0 - # 4. Handle case where VP is 0, it should behave like the tofts model - t = np.arange(0, 6 * 60, 1) - ca = osipi.aif_parker(t) - ct = osipi.two_compartment_exchange_model(t, ca, Fp=10, PS=5, ve=0.2, vp=0) - ct_tofts = osipi.tofts(t, ca, Ktrans=3.93, ve=0.2) +def test_2cxm_time_offset(time_array, aif): + """Test that the tissue concentration is shifted from the AIF by the specified offset time""" + ct = osipi.two_compartment_exchange_model(time_array, aif, Fp=10, PS=5, ve=0.2, vp=0.3, Ta=60.0) + assert (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 == 60.0 + + +def test_2cxm_zero_vp_matches_tofts(time_array, aif): + """Test that 2CXM with vp=0 behaves like the Tofts model""" + ct = osipi.two_compartment_exchange_model(time_array, aif, Fp=10, PS=5, ve=0.2, vp=0) + ct_tofts = osipi.tofts(time_array, aif, Ktrans=3.93, ve=0.2) assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) if __name__ == "__main__": - test_tissue_tofts() - test_tissue_extended_tofts() - test_tissue_two_compartment_exchange_model() - - print("All tissue concentration model tests passed!!") + pytest.main([__file__]) From e6e974b51a3a64ed476e824053b7d9d54b1b7405 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Mon, 3 Mar 2025 22:16:02 +0200 Subject: [PATCH 12/13] =?UTF-8?q?=F0=9F=94=A7=20=20Handle=20invalid=20inpu?= =?UTF-8?q?ts?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- docs/examples/tissue/plot_two_cxm.py | 2 +- src/osipi/_tissue.py | 19 +++++++++++++++ tests/test_tissue.py | 35 ++++++++++++++++++++++++++++ 3 files changed, 55 insertions(+), 1 deletion(-) diff --git a/docs/examples/tissue/plot_two_cxm.py b/docs/examples/tissue/plot_two_cxm.py index 8e608a3..3e0639b 100644 --- a/docs/examples/tissue/plot_two_cxm.py +++ b/docs/examples/tissue/plot_two_cxm.py @@ -28,7 +28,7 @@ # of 0.2, plasma volume fraction of 0.1, permeability serface area of 5 ml/min # and flow rate of 10 ml/min PS = 15 # Permeability surface area product in ml/min -Fp = [10, 25] # Flow rate in ml/min +Fp = [5, 25] # Flow rate in ml/min ve = 0.1 # Extracellular volume fraction vp = [0.1, 0.02] # Plasma volume fraction diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index dd85f08..15bdf5f 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -396,6 +396,25 @@ def two_compartment_exchange_model( >>> plt.plot(t, ca, "r", t, ct, "g") """ + # Input validation + if t.size != ca.size: + raise ValueError(f"t and ca must have the same size, got {t.size} and {ca.size}") + + for param, name in [(Fp, "Fp"), (PS, "PS"), (ve, "ve"), (vp, "vp")]: + if not (np.isscalar(param) and np.isreal(param)): + raise TypeError(f"{name} must be a numeric scalar value, got {type(param).__name__}") + + if Fp <= 0: + raise ValueError(f"Fp must be positive, got {Fp}") + if PS < 0: + raise ValueError(f"PS must be non-negative, got {PS}") + if not (0 <= ve <= 1): + raise ValueError(f"ve must be in range [0, 1], got {ve}") + if not (0 <= vp <= 1): + raise ValueError(f"vp must be in range [0, 1], got {vp}") + if (ve + vp) > 1: + raise ValueError(f"Sum of ve ({ve}) and vp ({vp}) exceeds 1") + if vp == 0: E = 1 - np.exp(-PS / Fp) Ktrans = E * Fp diff --git a/tests/test_tissue.py b/tests/test_tissue.py index 36764ac..e2f413d 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -162,5 +162,40 @@ def test_2cxm_zero_vp_matches_tofts(time_array, aif): assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) +def test_2cxm_invalid_inputs(time_array, aif): + """Test invalid inputs raise appropriate errors""" + # Test mismatched array sizes + with pytest.raises(ValueError, match="t and ca must have the same size"): + osipi.two_compartment_exchange_model(time_array[:10], aif, Fp=10, PS=5, ve=0.2, vp=0.1) + + # Test invalid parameter types and values + invalid_params = [ + ({"Fp": [10]}, TypeError, "Fp must be a numeric scalar value"), # Non-scalar Fp + ({"PS": "five"}, TypeError, "PS must be a numeric scalar value"), # Non-numeric PS + ({"ve": np.array([0.2])}, TypeError, "ve must be a numeric scalar value"), # Array ve + ({"vp": None}, TypeError, "vp must be a numeric scalar value"), # None value + ({"Fp": -5}, ValueError, "Fp must be positive"), # Negative Fp + ({"PS": -1}, ValueError, "PS must be non-negative"), # Negative PS + ({"ve": -0.1}, ValueError, "ve must be in range \[0, 1\]"), # ve < 0 + ({"vp": 1.1}, ValueError, "vp must be in range \[0, 1\]"), # vp > 1 + ( + {"ve": 0.8, "vp": 0.3}, + ValueError, + "Sum of ve \(0.8\) and vp \(0.3\) exceeds 1", + ), # ve + vp > 1 + ] + + for params, err_type, err_msg in invalid_params: + with pytest.raises(err_type, match=err_msg): + osipi.two_compartment_exchange_model( + time_array, + aif, + Fp=10 if "Fp" not in params else params["Fp"], + PS=5 if "PS" not in params else params["PS"], + ve=0.2 if "ve" not in params else params["ve"], + vp=0.1 if "vp" not in params else params["vp"], + ) + + if __name__ == "__main__": pytest.main([__file__]) From 3c679eb86824a261577924cef4702be1e479499b Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Sat, 22 Mar 2025 17:35:37 +0200 Subject: [PATCH 13/13] Add negative tests for each model --- src/osipi/_tissue.py | 36 +++++- tests/test_tissue.py | 268 +++++++++++++++++++++++++------------------ 2 files changed, 189 insertions(+), 115 deletions(-) diff --git a/src/osipi/_tissue.py b/src/osipi/_tissue.py index 15bdf5f..dbc3e9f 100755 --- a/src/osipi/_tissue.py +++ b/src/osipi/_tissue.py @@ -74,6 +74,22 @@ def tofts( >>> plt.plot(t, ca, "r", t, ct, "b") """ + + # Input validation + if t.size != ca.size: + raise ValueError(f"t and ca must have the same size, got {t.size} and {ca.size}") + + for param, name in [(Ktrans, "Ktrans"), (ve, "ve"), (Ta, "Ta")]: + if not (np.isscalar(param) and np.isreal(param)): + raise TypeError(f"{name} must be a numeric scalar value, got {type(param).__name__}") + + if Ktrans < 0: + raise ValueError(f"Ktrans must be non-negative, got {Ktrans}") + if not (0 <= ve <= 1): + raise ValueError(f"ve must be in range [0, 1], got {ve}") + if Ta < 0: + raise ValueError(f"Ta must be non-negative, got {Ta}") + if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( ("Non-uniform time spacing detected. Time array may be" " resampled."), @@ -238,6 +254,24 @@ def extended_tofts( >>> plt.plot(t, ca, "r", t, ct, "b") """ + # Input validation + if t.size != ca.size: + raise ValueError(f"t and ca must have the same size, got {t.size} and {ca.size}") + + for param, name in [(Ktrans, "Ktrans"), (ve, "ve"), (vp, "vp")]: + if not (np.isscalar(param) and np.isreal(param)): + raise TypeError(f"{name} must be a numeric scalar value, got {type(param).__name__}") + + if not (0 <= ve <= 1): + raise ValueError(f"ve must be in range [0, 1], got {ve}") + if not (0 <= vp <= 1): + raise ValueError(f"vp must be in range [0, 1], got {vp}") + if not (0 <= ve + vp <= 1): + raise ValueError(f"Sum of ve ({ve}) and vp ({vp}) exceeds 1") + if Ta < 0: + raise ValueError(f"Ta must be non-negative, got {Ta}") + if Ktrans < 0: + raise ValueError(f"Ktrans must be non-negative, got {Ktrans}") if not np.allclose(np.diff(t), np.diff(t)[0]): warnings.warn( @@ -412,7 +446,7 @@ def two_compartment_exchange_model( raise ValueError(f"ve must be in range [0, 1], got {ve}") if not (0 <= vp <= 1): raise ValueError(f"vp must be in range [0, 1], got {vp}") - if (ve + vp) > 1: + if not (0 <= ve + vp <= 1): raise ValueError(f"Sum of ve ({ve}) and vp ({vp}) exceeds 1") if vp == 0: diff --git a/tests/test_tissue.py b/tests/test_tissue.py index e2f413d..3e5fb5f 100755 --- a/tests/test_tissue.py +++ b/tests/test_tissue.py @@ -15,123 +15,159 @@ def aif(time_array): return osipi.aif_parker(time_array) -@pytest.fixture -def base_params(): - return {"Ktrans": 0.6, "ve": 0.2} - +# Tests for the Tofts model -def test_tissue_tofts(time_array, aif, base_params): - """1. Basic operation of the function - test that the peak tissue concentration is less than the peak - AIF - """ - ct = osipi.tofts(time_array, aif, discretization_method="exp", **base_params) +@pytest.mark.parametrize( + "t", + [np.arange(0, 6 * 60, 1), np.geomspace(1, 6 * 60 + 1, num=360) - 1], + ids=["uniform", "non-uniform"], +) +def test_tofts_peak_concentration(t, aif): + """Test peak tissue concentration is less than peak AIF""" + ct = osipi.tofts(t, aif, Ktrans=0.6, ve=0.2) assert np.round(np.max(ct)) < np.round(np.max(aif)) - # 2. Basic operation of the function - test with non-uniform spacing of - # time array - t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 - ca = osipi.aif_parker(t) - ct = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 3. The offset option - test that the tissue concentration is shifted - # from the AIF by the specified offset time - ct = osipi.tofts(time_array, aif, Ta=60.0, **base_params) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 == 60.0 +def test_tofts_time_offset(time_array, aif): + """Test tissue concentration is shifted by specified offset time""" + ct = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0.2, Ta=60.0) + offset_diff = (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 + assert offset_diff == 60.0 + - # 4. Test that the discretization options give almost the same result - - # time step must be very small +def test_tofts_discretization_agreement(): + """Test different discretization methods give similar results with small time steps""" t = np.arange(0, 6 * 60, 0.01) ca = osipi.aif_parker(t) - ct_conv = osipi.tofts(t, ca, **base_params) - ct_exp = osipi.tofts(t, ca, discretization_method="exp", **base_params) + ct_conv = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2) + ct_exp = osipi.tofts(t, ca, Ktrans=0.6, ve=0.2, discretization_method="exp") assert np.allclose(ct_conv, ct_exp, rtol=1e-4, atol=1e-3) - # 5. Test that the ratio of the area under the ct and ca curves is - # approximately the extracellular volume - ct_conv = osipi.tofts(time_array, aif, **base_params) - ct_exp = osipi.tofts(time_array, aif, discretization_method="exp", **base_params) - assert math.isclose( - np.trapz(ct_conv, time_array) / np.trapz(aif, time_array), 0.2, abs_tol=1e-1 - ) - assert math.isclose(np.trapz(ct_exp, time_array) / np.trapz(aif, time_array), 0.2, abs_tol=1e-1) +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_tofts_volume_ratio(time_array, aif, method): + """Test ratio of AUC matches extracellular volume""" + ct = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0.2, discretization_method=method) + ratio = np.trapz(ct, time_array) / np.trapz(aif, time_array) + assert math.isclose(ratio, 0.2, abs_tol=1e-1) + - # 6. Test specific use cases - ct_conv = osipi.tofts(time_array, aif, Ktrans=0, ve=0.2) - assert np.count_nonzero(ct_conv) == 0 +@pytest.mark.parametrize("Ktrans, ve", [(0, 0.2), (0.6, 0)], ids=["Ktrans=0", "ve=0"]) +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_tofts_specific_cases(time_array, aif, Ktrans, ve, method): + """Test edge cases with parameter combinations""" + ct = osipi.tofts(time_array, aif, Ktrans=Ktrans, ve=ve, discretization_method=method) + assert np.count_nonzero(ct) == 0 - ct_exp = osipi.tofts(time_array, aif, Ktrans=0, ve=0.2, discretization_method="exp") - assert np.count_nonzero(ct_exp) == 0 - ct_conv = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0) - assert np.count_nonzero(ct_conv) == 0 +@pytest.mark.parametrize( + "param, value, match", + [ + ("Ktrans", -0.1, "non-negative"), + ("Ktrans", "invalid", "numeric scalar"), + ("ve", -0.2, "in range \[0, 1\]"), + ("ve", [0.2], "numeric scalar"), + ("Ta", -5.0, "non-negative"), + ], + ids=["Ktrans_negative", "Ktrans_type", "ve_negative", "ve_type", "Ta_negative"], +) +def test_tofts_invalid_parameters(time_array, aif, param, value, match): + """Test invalid parameter values/type for all model parameters""" + with pytest.raises((ValueError, TypeError), match=match): + # Construct kwargs with one invalid parameter at a time + kwargs = {"Ktrans": 0.6, "ve": 0.2, "Ta": 30.0} + kwargs[param] = value + osipi.tofts(time_array, aif, **kwargs) - ct_exp = osipi.tofts(time_array, aif, Ktrans=0.6, ve=0, discretization_method="exp") - assert np.count_nonzero(ct_exp) == 0 +# Tests for the Extended Tofts model -def test_tissue_extended_tofts(time_array, aif): - # 1. Basic operation of the function - test that the peak tissue - # concentration is less than the peak AIF - t = np.linspace(0, 6 * 60, 360) - ca = osipi.aif_parker(t) - ct = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) - assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 2. Basic operation of the function - test with non-uniform spacing of - # time array - t = np.geomspace(1, 6 * 60 + 1, num=360) - 1 +@pytest.mark.parametrize( + "t", + [np.linspace(0, 6 * 60, 360), np.geomspace(1, 6 * 60 + 1, num=360) - 1], + ids=["uniform", "non-uniform"], +) +def test_extended_tofts_peak_concentration(t): + """Test peak tissue concentration is less than peak AIF for different time arrays""" ca = osipi.aif_parker(t) ct = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) assert np.round(np.max(ct)) < np.round(np.max(ca)) - # 3. The offset option - test that the tissue concentration is shifted - # from the AIF by the specified offset time + +def test_extended_tofts_time_offset(time_array, aif): + """Test tissue concentration is shifted by specified offset time""" ct = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, Ta=60.0) - assert (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 == 60.0 + offset_diff = (np.min(np.where(ct > 0.0)) - np.min(np.where(aif > 0.0)) - 1) * 1 + assert offset_diff == 60.0 + - # 4. Test that the discretization options give almost the same result - - # time step must be very small +def test_extended_tofts_discretization_agreement(): + """Test different discretization methods give similar results with small time steps""" t = np.arange(0, 6 * 60, 0.01) ca = osipi.aif_parker(t) ct_conv = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3) ct_exp = osipi.extended_tofts(t, ca, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method="exp") assert np.allclose(ct_conv, ct_exp, rtol=1e-4, atol=1e-3) - # 5. Test that the ratio of the area under the ct and ca curves is - # approximately the extracellular volume plus the plasma volume - ct_conv = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3) - ct_exp = osipi.extended_tofts( - time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method="exp" - ) - assert math.isclose( - np.trapz(ct_conv, time_array) / np.trapz(aif, time_array), - 0.2 + 0.3, - abs_tol=1e-1, - ) - assert math.isclose( - np.trapz(ct_exp, time_array) / np.trapz(aif, time_array), 0.2 + 0.3, abs_tol=1e-1 - ) - # 6. Test specific use cases +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_extended_tofts_volume_ratio(time_array, aif, method): + """Test ratio of AUC matches ve + vp for different discretization methods""" + ct = osipi.extended_tofts( + time_array, aif, Ktrans=0.6, ve=0.2, vp=0.3, discretization_method=method + ) + ratio = np.trapz(ct, time_array) / np.trapz(aif, time_array) + assert math.isclose(ratio, 0.2 + 0.3, abs_tol=1e-1) - ct_conv = osipi.extended_tofts(time_array, aif, Ktrans=0, ve=0.2, vp=0.3) - assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) - ct_exp = osipi.extended_tofts( - time_array, aif, Ktrans=0, ve=0.2, vp=0.3, discretization_method="exp" +@pytest.mark.parametrize("Ktrans, ve", [(0, 0.2), (0.6, 0)], ids=["Ktrans=0", "ve=0"]) +@pytest.mark.parametrize("method", [None, "exp"], ids=["conv", "exp"]) +def test_extended_tofts_specific_cases(time_array, aif, Ktrans, ve, method): + """Test edge cases with parameter combinations""" + ct = osipi.extended_tofts( + time_array, aif, Ktrans=Ktrans, ve=ve, vp=0.3, discretization_method=method ) - assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) + expected = aif * 0.3 + assert np.allclose(ct, expected, rtol=1e-4, atol=1e-3) - ct_conv = osipi.extended_tofts(time_array, aif, Ktrans=0.6, ve=0, vp=0.3) - assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) - ct_exp = osipi.extended_tofts( - time_array, aif, Ktrans=0.6, ve=0, vp=0.3, discretization_method="exp" - ) - assert np.allclose(ct_conv, aif * 0.3, rtol=1e-4, atol=1e-3) +@pytest.mark.parametrize( + "param, value, match", + [ + ("Ktrans", -0.1, "non-negative"), + ("Ktrans", "invalid", "numeric scalar"), + ("ve", -0.2, "in range \[0, 1\]"), + ("ve", [0.2], "numeric scalar"), + ("vp", -0.2, "in range \[0, 1\]"), + ("vp", [0.2], "numeric scalar"), + ("vp", 1.1, "in range \[0, 1\]"), + ("vp", [0.2], "numeric scalar"), + ("Ta", -5.0, "non-negative"), + ], + ids=[ + "Ktrans_negative", + "Ktrans_type", + "ve_negative", + "ve_type", + "vp_negative", + "vp_type", + "vp_range", + "vp_type", + "Ta_negative", + ], +) +def test_extended_tofts_invalid_parameters(time_array, aif, param, value, match): + """Test invalid parameter values/type for all model parameters""" + with pytest.raises((ValueError, TypeError), match=match): + # Construct kwargs with one invalid parameter at a time + kwargs = {"Ktrans": 0.6, "ve": 0.2, "vp": 0.3, "Ta": 30.0} + kwargs[param] = value + osipi.extended_tofts(time_array, aif, **kwargs) + + +# Tests for the 2CXM model @pytest.mark.parametrize( @@ -162,39 +198,43 @@ def test_2cxm_zero_vp_matches_tofts(time_array, aif): assert np.allclose(ct, ct_tofts, rtol=1e-4, atol=1e-3) -def test_2cxm_invalid_inputs(time_array, aif): - """Test invalid inputs raise appropriate errors""" - # Test mismatched array sizes - with pytest.raises(ValueError, match="t and ca must have the same size"): - osipi.two_compartment_exchange_model(time_array[:10], aif, Fp=10, PS=5, ve=0.2, vp=0.1) - - # Test invalid parameter types and values - invalid_params = [ - ({"Fp": [10]}, TypeError, "Fp must be a numeric scalar value"), # Non-scalar Fp - ({"PS": "five"}, TypeError, "PS must be a numeric scalar value"), # Non-numeric PS - ({"ve": np.array([0.2])}, TypeError, "ve must be a numeric scalar value"), # Array ve - ({"vp": None}, TypeError, "vp must be a numeric scalar value"), # None value - ({"Fp": -5}, ValueError, "Fp must be positive"), # Negative Fp - ({"PS": -1}, ValueError, "PS must be non-negative"), # Negative PS - ({"ve": -0.1}, ValueError, "ve must be in range \[0, 1\]"), # ve < 0 - ({"vp": 1.1}, ValueError, "vp must be in range \[0, 1\]"), # vp > 1 - ( - {"ve": 0.8, "vp": 0.3}, - ValueError, - "Sum of ve \(0.8\) and vp \(0.3\) exceeds 1", - ), # ve + vp > 1 - ] - - for params, err_type, err_msg in invalid_params: - with pytest.raises(err_type, match=err_msg): - osipi.two_compartment_exchange_model( - time_array, - aif, - Fp=10 if "Fp" not in params else params["Fp"], - PS=5 if "PS" not in params else params["PS"], - ve=0.2 if "ve" not in params else params["ve"], - vp=0.1 if "vp" not in params else params["vp"], - ) +@pytest.mark.parametrize( + "param, value, match", + [ + ("Fp", [10], "numeric scalar"), # Non-scalar Fp + ("PS", "five", "numeric scalar"), # Non-numeric PS + ("ve", np.array([0.2]), "numeric scalar"), # Array ve + ("vp", None, "numeric scalar"), # None value + ("Fp", -5, "positive"), # Negative Fp + ("PS", -1, "non-negative"), # Negative PS + ("ve", -0.1, "in range \[0, 1\]"), # ve < 0 + ("vp", 1.1, "in range \[0, 1\]"), # vp > 1 + (("ve", "vp"), (0.8, 0.3), "Sum of ve \(0.8\) and vp \(0.3\) exceeds 1"), + ], + ids=[ + "Fp_type", + "PS_type", + "ve_type", + "vp_type", + "Fp_negative", + "PS_negative", + "ve_range", + "vp_range", + "ve_vp_sum", + ], +) +def test_2cxm_invalid_inputs(time_array, aif, param, value, match): + """Test invalid parameter values/type for all model parameters""" + with pytest.raises((ValueError, TypeError), match=match): + # Construct kwargs with default valid values + kwargs = {"Fp": 10, "PS": 5, "ve": 0.2, "vp": 0.1} + # Handle tuple parameters for combined validation case + if isinstance(param, tuple): + for p, v in zip(param, value): + kwargs[p] = v + else: + kwargs[param] = value + osipi.two_compartment_exchange_model(time_array, aif, **kwargs) if __name__ == "__main__":