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__":