Introducing Two Compartment exchange Model

.pre-commit-config.yaml

@@ -1,9 +1,9 @@
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.4.8  # Use the latest version
+    rev: v0.4.8 # Use the latest version
     hooks:
       - id: ruff
-        args: [--fix]  # Optional: to enable lint fixes
+        args: [--fix] # Optional: to enable lint fixes
       - id: ruff-format
   - repo: https://github.com/pre-commit/pre-commit-hooks
     rev: v4.6.0

docs/examples/tissue/plot_two_cxm.py

@@ -0,0 +1,51 @@
+"""
+=============
+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, 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.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()
+plt.show()

docs/references/models/tissue_models/2cxm.md

@@ -0,0 +1,17 @@
+# osipi.Two Compartment Exchange Model
+
+::: osipi.two_compartment_exchange_model
+
+
+## Example using `osipi.two_compartment_exchange_model`
+
+<figure class="mkd-glr-thumbcontainer" tooltip="Simulating tissue concentrations from two compartment exchange model with different settings.">
+    <img alt="The Two Compartment Exchange Model" src="..\..\..\..\generated\gallery\tissue\images\thumb\mkd_glr_plot_two_cxm_thumb.png" />
+    <figcaption class="caption">
+        <span class="caption-text">
+            <a class="reference internal" href="..\..\..\..\generated\gallery\tissue\plot_two_cxm">
+                <span class="std std-ref">The Two Compartment Exchange Model</span>
+            </a>
+        </span>
+    </figcaption>
+</figure>

docs/references/models/tissue_models/index.md

@@ -1,3 +1,4 @@
 
 - [tofts](tofts.md)
 - [extended_tofts](extended_tofts.md)
+- [two_cxm](2cxm.md)

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
 

src/osipi/__init__.py

@@ -7,7 +7,8 @@
 
 from ._tissue import (
     tofts,
-    extended_tofts
+    extended_tofts,
+    two_compartment_exchange_model
 )
 
 from ._signal import (

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
 
@@ -332,3 +333,146 @@ def extended_tofts(
                 ct = ct_func(t)
 
     return ct
+
+
+def two_compartment_exchange_model(
+    t: np.ndarray,
+    ca: np.ndarray,
+    Fp: float,
+    PS: float,
+    ve: float,
+    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). [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 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."),
+            stacklevel=2,
+        )
+
+    # Convert units
+    fp_per_s = Fp / (60.0 * 100.0)
+    ps_per_s = PS / (60.0 * 100.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
+
+    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 + 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
+        * 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)
+        )
+        / (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)) / upsample_factor
+
+    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)]
+
+    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

tests/test_tissue.py

@@ -124,8 +124,39 @@ def test_tissue_extended_tofts():
     assert np.allclose(ct_conv, ca * 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
+    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__":
     test_tissue_tofts()
     test_tissue_extended_tofts()
+    test_tissue_two_compartment_exchange_model()
 
     print("All tissue concentration model tests passed!!")