From 50cd6656122ffb6422c5c6acfb0d8013571cf047 Mon Sep 17 00:00:00 2001 From: Mohamed Nasser Date: Fri, 15 Nov 2024 23:34:19 +0200 Subject: [PATCH] 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