simulating MRI signal

added forward model for C to R1
added example in user guide for simulating an MRI signal

docs/user-guide/gen_aif.md

@@ -9,6 +9,8 @@ import osipi
 t = np.arange(0, 6*60, 1)
 ca = osipi.aif_parker(t)
 plt.plot(t, ca)
+plt.xlabel('Time (s)')
+plt.ylabel('Indicator concentration (mM)')
 plt.show()
 ```
 
@@ -24,13 +26,55 @@ ca = osipi.aif_parker(t)
 Ktrans = 0.6
 ve = 0.2
 ct = osipi.tofts(t, ca, Ktrans=Ktrans/60, ve=ve)
-plt.plot(t, ct)
+fig, ax = plt.subplots(1, 2)
+ax[0].plot(t, ca)
+ax[0].set_xlabel('Time (s)')
+ax[0].set_ylabel('Indicator concentration (mM)')
+ax[0].set_title('AIF')
+ax[1].plot(t, ct)
+ax[1].set_xlabel('Time (s)')
+ax[1].set_ylabel('Indicator concentration (mM)')
+ax[1].set_title('Tissue')
+
+fig.tight_layout()
 plt.show()
 ```
 
 ## Generating an MRI signal
-!!! note "Coming Soon"
-    This section is under development and will be available soon.
+
+``` py
+import numpy as np
+import matplotlib.pyplot as plt
+import osipi
+
+t = np.arange(0, 6*60, 1)
+ca = osipi.aif_parker(t)
+Ktrans = 0.6
+ve = 0.2
+ct = osipi.tofts(t, ca, Ktrans=Ktrans/60, ve=ve)
+
+R10 = 0.5
+r1 = 4.5
+R1t = osipi.C_to_R1_linear_relaxivity(ct, R10, r1)
+
+S0 = 1
+TR = 0.004
+a = 12
+St = osipi.signal_SPGR(R1t, S0, TR, a)
+
+fig, ax = plt.subplots(1, 2)
+ax[0].plot(t, ct)
+ax[0].set_xlabel('Time (s)')
+ax[0].set_ylabel('Indicator concentration (mM)')
+ax[0].set_title('Concentration')
+ax[1].plot(t, St)
+ax[1].set_xlabel('Time (s)')
+ax[1].set_ylabel('MRI signal (a.u)')
+ax[1].set_title('MRI signal')
+
+fig.tight_layout()
+plt.show()
+```
 
 ## Adding measurement error
 !!! note "Coming Soon"

src/osipi/__init__.py

@@ -17,6 +17,10 @@
 
 from ._signal_to_concentration import (
     S_to_C_via_R1_SPGR,
-    S_to_R1_SPGR,
-    R1_to_C_linear_relaxivity
+    S_to_R1_SPGR
+)
+
+from ._electromagnetic_property import (
+    R1_to_C_linear_relaxivity,
+    C_to_R1_linear_relaxivity
 )

src/osipi/_electromagnetic_property.py

@@ -38,3 +38,37 @@ def R1_to_C_linear_relaxivity(
     elif not (r1 >= 0):
         raise ValueError("r1 must be positive")
     return (R1 - R10) / r1  # C
+
+
+def C_to_R1_linear_relaxivity(
+    C: NDArray[np.float64], R10: np.float64, r1: np.float64
+) -> NDArray[np.float64]:
+    """
+    Electromagnetic property forward model:
+    - longitudinal relaxation rate, linear with relaxivity
+
+    Converts tissue concentration to R1
+
+    Args:
+        C (1D array of np.float64):
+            Vector of indicator concentrations in units of mM. [OSIPI code Q.IC1.001]
+        R10 (np.float64):
+            Native longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.002]
+        r1 (np.float64):
+            Longitudinal relaxivity in units of /s/mM. [OSIPI code Q.EL1.015]
+
+    Returns:
+        NDArray[np.float64]:
+            Vector of longitudinal relaxation rate in units of /s. [OSIPI code Q.EL1.001]
+
+    References:
+        - Lexicon URL: https://osipi.github.io/OSIPI_CAPLEX/perfusionModels/#
+        - Lexicon code: M.EL1.003
+        - Adapted from equation given in lexicon
+    """
+    # Check C is a 1D array of floats
+    if not (isinstance(C, np.ndarray) and C.ndim == 1 and C.dtype == np.float64):
+        raise TypeError("C must be a 1D NumPy array of np.float64")
+    elif not (r1 >= 0):
+        raise ValueError("r1 must be positive")
+    return R10 + r1 * C  # R1