Fixed bug and added tests around molality scaling.

CHANGELOG.rst

@@ -2,6 +2,7 @@ This document contains the FSL-MRS release history in reverse chronological orde
 
 1.1.4 (WIP)
 ------------------------------
+- Fixed bug in calculation of molality concentration. Tissue mole fractions had been swapped for tissue volume fractions. Molar concentrations unaffected.
 - Fixed bug in mrs_tools split
 - Fixed bug in alignment of multi-dimensional data.
 - fsl_mrsi now outputs fitting nuisance parameters: phases, and shifts & linewidths for each metabolite group.

fsl_mrs/tests/test_utils_quantify.py

@@ -60,6 +60,59 @@ def test_QuantificationInfo():
     assert qci.add_corr == 5.0
 
 
+def test_volumefraction_calc():
+    qci = quant.QuantificationInfo(0.010, 3, ['NAA'], 298)
+    qci.set_fractions({'GM': 0.45, 'WM': 0.40, 'CSF': 0.15})
+    assert qci.f_GM == 0.45
+    assert qci.f_WM == 0.40
+    assert qci.f_CSF == 0.15
+
+
+def test_molefraction_calc():
+    qci = quant.QuantificationInfo(0.010, 3, ['NAA'], 298)
+    qci.set_fractions({'GM': 0.45, 'WM': 0.40, 'CSF': 0.15})
+
+    # Densitites are 'GM': 0.78, 'WM': 0.65, 'CSF': 0.97
+    sum_frac = (0.45 * 0.78 + 0.40 * 0.65 + 0.15 * 0.97)
+    assert np.isclose(qci.f_GM_H2O, 0.45 * 0.78 / sum_frac)
+    assert np.isclose(qci.f_WM_H2O, 0.40 * 0.65 / sum_frac)
+    assert np.isclose(qci.f_CSF_H2O, 0.15 * 0.97 / sum_frac)
+
+
+def test_corrected_water_conc():
+    # No relaxation
+    qci = quant.QuantificationInfo(1E-10, 1E5, ['NAA'], 298)
+    qci.set_fractions({'GM': 1.00, 'WM': 0.0, 'CSF': 0.0})
+
+    print(qci.relax_corr_water_molal)
+    print(qci.relax_corr_water_molar)
+    # Molality should be close to pure water as density term cancels
+    assert np.isclose(qci.relax_corr_water_molal, 55510)
+    # Molarity should be scaled by density term as volume fixed
+    assert np.isclose(qci.relax_corr_water_molar, 55510 * 0.78)
+
+    qci.set_fractions({'GM': 0.50, 'WM': 0.5, 'CSF': 0.0})
+
+    print(qci.relax_corr_water_molal)
+    print(qci.relax_corr_water_molar)
+    # Molality should be close to pure water as density term cancels
+    assert np.isclose(qci.relax_corr_water_molal, 55510)
+    # Molarity should be scaled by density terms as volume fixed
+    assert np.isclose(qci.relax_corr_water_molar, 55510 * (0.78 + 0.65) / 2)
+
+    qci = quant.QuantificationInfo(1E-10, 1, ['NAA'], 298)
+    qci.set_fractions({'GM': 0.50, 'WM': 0.5, 'CSF': 0.0})
+
+    print(qci.relax_corr_water_molal)
+    print(qci.relax_corr_water_molar)
+    # Molality should scaled by relaxation terms in proportion of mole fraction.
+    mf_gm = 0.5 * 0.78 / (0.5 * 0.78 + 0.5 * 0.65)
+    mf_wm = 0.5 * 0.65 / (0.5 * 0.78 + 0.5 * 0.65)
+    assert np.isclose(qci.relax_corr_water_molal, 55510 * (qci.R_H2O_GM * mf_gm + qci.R_H2O_WM * mf_wm))
+    # Molarity should be scaled by density terms * relaxation terms
+    assert np.isclose(qci.relax_corr_water_molar, 55510 * (0.78 * qci.R_H2O_GM + 0.65 * qci.R_H2O_WM) / 2)
+
+
 def test_quantifyWater():
     basis = mrsio.read_basis(basisfile)
     data = mrsio.read_FID(metabfile)
@@ -107,3 +160,4 @@ def test_quantifyWater():
 
     assert np.allclose(res.getConc(scaling='internal'), 1.0)
     assert np.allclose(res.getConc(scaling='molarity'), 10.78, atol=1E-1)
+    assert np.allclose(res.getConc(scaling='molality'), 10.78 * 1 / (0.6 * 0.78 + 0.4 * 0.65), atol=1E-1)

fsl_mrs/utils/quantify.py

@@ -321,7 +321,7 @@ class QuantificationInfo(object):
     def d_CSF(self):
         return self._densitites['CSF']
 
-    # tissue fractions f_GM, f_WM, f_CSF
+    # tissue volume fractions f_GM, f_WM, f_CSF
     @property
     def f_GM(self):
         if self._fractions:
@@ -343,6 +343,54 @@ class QuantificationInfo(object):
         else:
             return None
 
+    # tissue mole fractions f_GM_H2O, f_WM_H2O, f_CSF_H2O
+    # Implementing equation 5 in https://doi.org/10.1002/nbm.4257
+    # f_X_H2O = f_X*d_X / (f_GM*d_GM + f_WM*d_WM f_CSF*d_CSF)
+    @property
+    def f_GM_H2O(self):
+        """Grey matter (GM) mole fraction
+
+        :return: GM mole fraction
+        :rtype: float
+        """
+        if self._fractions:
+            return self.f_GM * self.d_GM\
+                / (self.f_GM * self.d_GM
+                   + self.f_WM * self.d_WM
+                   + self.f_CSF * self.d_CSF)
+        else:
+            return None
+
+    @property
+    def f_WM_H2O(self):
+        """White matter (WM) mole fraction
+
+        :return: WM mole fraction
+        :rtype: float
+        """
+        if self._fractions:
+            return self.f_WM * self.d_WM\
+                / (self.f_GM * self.d_GM
+                   + self.f_WM * self.d_WM
+                   + self.f_CSF * self.d_CSF)
+        else:
+            return None
+
+    @property
+    def f_CSF_H2O(self):
+        """Cerebral spinal fluid (CSF) mole fraction
+
+        :return: CSF mole fraction
+        :rtype: float
+        """
+        if self._fractions:
+            return self.f_CSF * self.d_CSF\
+                / (self.f_GM * self.d_GM
+                   + self.f_WM * self.d_WM
+                   + self.f_CSF * self.d_CSF)
+        else:
+            return None
+
     # Relaxation parameters
     # R_H2O_GM, R_H2O_WM, R_H2O_CSF, R_H2O, R_M
     # Assumes TR >> TE and excitation FA = 90 deg
@@ -365,6 +413,7 @@ class QuantificationInfo(object):
         return self._calc_relax(self.t1['H2O_CSF'], self.t2['H2O_CSF'])
 
     # A 50/50 average between WM/GM
+    # Used in the absence of fractions
     @property
     def R_H2O(self):
         t1 = (self.t1['H2O_WM'] + self.t1['H2O_GM']) / 2.0
@@ -377,17 +426,33 @@ class QuantificationInfo(object):
 
     # Water concentrations
     # relax_corr_water_molal, relax_corr_water_molar
+    # relax_corr_water_molal uses mole fractions
+    # relax_corr_water_molar uses volume fractions * densitites
     @property
     def relax_corr_water_molal(self):
+        """Relaxation (T1, T2) corrected water molality (mmol/kg).
+        If volume fractions aren't availible then relaxation correction will be based on
+        a 50/50 split of GM/WM T1/T2s and pure water will be assumed.
+
+        :return: concentration
+        :rtype: float
+        """
         if self._fractions is None:
             return self.R_H2O * H2O_MOLALITY
         else:
-            return H2O_MOLALITY * (self.f_GM * self.R_H2O_GM
-                                   + self.f_WM * self.R_H2O_WM
-                                   + self.f_CSF * self.R_H2O_CSF)
+            return H2O_MOLALITY * (self.f_GM_H2O * self.R_H2O_GM
+                                   + self.f_WM_H2O * self.R_H2O_WM
+                                   + self.f_CSF_H2O * self.R_H2O_CSF)
 
     @property
     def relax_corr_water_molar(self):
+        """Relaxation (T1, T2) corrected water molariyt (mmol/dm^3 = mM).
+        If volume fractions aren't availible then relaxation correction will be based on
+        a 50/50 split of GM/WM T1/T2s and pure water will be assumed.
+
+        :return: concentration
+        :rtype: float
+        """
         if self._fractions is None:
             return self.R_H2O * H2O_MOLALITY
         else:
@@ -477,8 +542,8 @@ def quantifyWater(mrs, results, quant_info, verbose=False):
 
     # Calculate concentration scalings
     # EQ 4 and 6 in https://doi.org/10.1002/nbm.4257
-    # conc_molal =  (SMObs *(Q.f_GM*Q.R_H2O_GM + Q.f_WM*Q.R_H2O_WM + Q.f_CSF*Q.R_H2O_CSF)\
-    #                       / (SH2OObs*(1-Q.f_CSF)*Q.R_M)) \
+    # conc_molal =  (SMObs *(Q.f_GM_H20*Q.R_H2O_GM + Q.f_WM_H20*Q.R_H2O_WM + Q.f_CSF_H20*Q.R_H2O_CSF)\
+    #                       / (SH2OObs*(1-Q.f_CSF_H20)*Q.R_M)) \
     #                 * (H2O_PROTONS/refProtons)\
     #                 * H2O_MOLALITY
 
@@ -487,6 +552,9 @@ def quantifyWater(mrs, results, quant_info, verbose=False):
     #                 * (H2O_PROTONS/refProtons)\
     #                 * H2O_MOLALITY
 
+    # Note the difference between Q.f_X and Q.f_X_H2O. Equation 5 of reference. With thanks to Alex Craig-Craven
+    # for pointing this out.
+
     conc_molal = (SMObs / SH2OObs) * (H2O_PROTONS / quant_info.ref_protons) * \
         quant_info.relax_corr_water_molal * quant_info.csf_corr * quant_info.add_corr * quant_info.relax_corr_metab
     conc_molar = (SMObs / SH2OObs) * (H2O_PROTONS / quant_info.ref_protons) * \
@@ -497,7 +565,8 @@ def quantifyWater(mrs, results, quant_info, verbose=False):
         metabRelaxCorr = quant_info.relax_corr_metab
         print(f'Metabolite area = {SMObs:0.2e}')
         print(f'Water area = {SH2OObs:0.2e}')
-        print(f'Relaxation corrected water concentration = {rcorwaterconc:0.2e}')
+        print(f'Relaxation corrected water concentration (molal) = {quant_info.relax_corr_water_molal:0.2e} mmol/kg')
+        print(f'Relaxation corrected water concentration (molar) = {rcorwaterconc:0.2e} mM')
         print(f'metabolite relaxation correction  = {metabRelaxCorr:0.2e}')
         print(f'H2O to ref molality scaling = {conc_molal:0.2e}')
         print(f'H2O to ref molarity scaling = {conc_molar:0.2e}')