Added initial implementation for two compartment exchange model

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()

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