Under the assumption that \(_\) and \(_^\) are constant within the acquisition time window (e.g., not influenced by exchange processes), the MR signals of ME-GE and ME-SE sequences can be expressed as integrals [3]

$$S_ \left( t \right) = \mathop \int \limits_^ \,\Omega_ \left( ^ } \right)\,\exp \left[ ^ } \right]\;dT_^$$

(1)

$$S_ \left( t \right) = \mathop \int \limits_^ \,\Omega_ \left( } \right)\,S_ \left( ,T_ ,T_ } \right)\;}T_$$

(2)

where \(_\) and \(_\) denote the \(_\) and \(_^\) relaxation spectra (or distributions) of an imaging voxel. The first equation assumes that the effective transverse relaxation is governed by a simple multi-exponential decay. Thus, the composite ME-GE signal is the forward Laplace transform of the \(_^\) distribution. In case of perfect 180° refocusing, the ME-SE signal would follow the same signal model with \(_\) instead of \(_^\). However, it is typically not possible to obtain perfect refocusing pulses in real experiments due to subject-induced transmit field inhomogeneities. Therefore, it is necessary to account for signal deviations originating from stimulated echoes generated in the CPMG echo train [12]. These signal components can be conveniently modeled with the extended phase graph (EPG) algorithm [13, 14]. Accordingly, \(_\left(t;_,_,_\right)\) in Eq. (2) denotes the ME-SE signal at the echo time \(t\) computed with the EPG formalism. It depends on the refocusing flip angle, \(_\), and the relaxation times \(_\) and \(_\). However, the \(_\) dependence is weak in case of brain tissue with \(_\gg _\), and it is sufficient to approximate \(_\approx 1\) s with minimal effect on the signal amplitude (and, therefore, subsequent MWF estimation) [15].

Discretization of the integrals in Eqs. (1) and (2) for a finite number of echo times and a grid of relaxation times, \(_\) and \(_^\), leads to linear equations

$$S_ \left( } \right) = A_ \left( ,T_^ } \right) \cdot \Omega_ \left( ^ } \right)$$

(3)

$$S_ \left( } \right) = A_ \left( ,T_ } \right) \cdot \Omega_ \left( } \right)$$

(4)

where \(_\) and \(_\) are operator matrices solving the respective forward problems, \(_\) and \(_\) are signal vectors at respective echo times \(_\) and \(_\), and \(_\) and \(_\) are the relaxation spectra vectors on a suitable grid of relaxation times, \(_\) and \(_^\), respectively. The forward operators, \(_\) and \(_\), can be precomputed for efficiently solving the nonlinear inverse problem.

The ME-GE exponential signal model of Eq. (1) does not account for compartmental frequency shifts within the imaging voxel. It was previously shown that susceptibility-induced frequency shifts give rise to a non-exponential ME-GE signal [16], resulting in biased MWF estimates from ME-GE magnitude data. Therefore, we extend the ME-GE signal model to account for frequency shifts, \(\Delta \omega\), in the MW and AEW compartments. Following the approach from Nam et al. [7], we extend Eq. (3) to a complex signal model by introducing a complex forward operator

$$_\left(t,_^,\Delta \omega ,\Delta _\right)=\text\left[-t/_^\right]\text\left[-}\left(\Delta \omega t+\Delta _\right)\right]\hspace$$

(5)

where \(\Delta _\) is a spatially dependent constant phase offset. In the following parameterized forward model, we assume constant frequency shifts, \(\Delta _\) and \(\Delta _\), for the MW and AEW compartments, respectively [7, 16].

We first present the method for the single inversion, i.e., the method to recover the MWF from \(_\) (respectively \(_^\)) data only, before describing the joint inversion. To model the \(_\) (respectively \(_^\)) distribution, we resort to a parametric model, assuming the presence of two compartments with different \(_\) values, one for the MW pool (fast compartment) and one for axonal/extra-cellular water pool (slow compartment). Similar to [6], we propose a bi-Gaussian parametric model, which strongly improves the optimization stability compared to a bi-impulse (delta function) model.

Within this assumption, the \(_\) distribution is modeled by six parameters: the \(_\) means (\(_\), \(_\)) and standard deviations (\(_\), \(_\)) of both compartments, the integral of the axonal/extra-cellular compartment (\(_\)), and finally the \(MWF\), which is the ratio of the integral of the myelin compartment over the sum of the two compartment integrals, \(_/\left(_+_\right)\). Most references in the literature, using either \(_\) or \(_^\) parametric models, first invert for all compartment integrals (or amplitudes if delta functions are used instead of Gaussian functions) and based on that compute MWF [2, 7, 17]. Instead, we used the MWF directly as one of the model parameters for the inversion. We found that this parametrization leads to a slightly better resolution and spatial continuity on the final 2D MWF map than inverting for both compartment integrals separately.

Let \(G\left(\tau |\mu ,\sigma ,I\right)\) be a Gaussian function over the variable \(\tau\), defined by its mean value \(\mu\), standard deviation \(\sigma\) and integral value \(I\). Denoting \(x=\left(_,_,MWF,_,_,_\right)\) as the model vector i.e., the six parameters describing the two Gaussian functions, the \(_\) distribution is given by

$$_\left(_|x\right)=G\left(_|_,_,_\right)+G\left(_|_,_,_\right)\hspace$$

(6)

where \(_\) and \(_\) are the integrals of the MW pool and the axonal/extra-cellular water pool (AEW), respectively. \(_\) is not treated as a model parameter, but expressed according to the MWF definition as \(_=_MWF/\left(1-MWF\right)\). Let \(_\left(t\right)\) be the observed ME-SE signal with \(t\) being the echo times vector. For a single inversion, the unregularized minimization problem is then given by

$$\mathop \limits_ \;\left\| \left( t \right) - \hat_ \left( \right)^ } \right\|$$

where \(}_\left(t|x\right)=_\left(t,_\right)\cdot _\left(_|x\right)\) is the parameterized forward model from the \(_\) relaxometry distribution to the data space (ME-SE decay curve).

Since the problem is highly ill-posed, regularization has to be added to stabilize the inversion. However, classical Tikhonov regularization aims to minimize the model parameters and, therefore, results in low and unrealistic values for the myelin and axonal/extra-cellular \(_\) values (\(_\) and \(_\)). To overcome this problem, we added a constraint on the model to stay as close as possible to the initial model, but weighting this constraint mainly on the mean of the parameters. The regularized minimization problem is given by

$$\mathop \limits_ \left\| \;S_ \left( t \right) - \hat_ \left( \right)\right\|^ + \left\| \lambda \cdot \left( } \right)\right\|^$$

(7)

where \(_}\) is the initial model vector (cf. Table 1) and \(\lambda\) a model-sized regularization vector. We resort to this regularization for the fast and slow compartments mean values (\(_\),\(_\)) since we do have prior information on these parameters from the literature, and want to use this knowledge as initial values to better constraint the inversion. Regarding the standard deviations, which are initialized with small values (see Tab. 1), this regularization helps to maintain the values as small as possible, since we would expect delta functions in an ideal noise-free world. Finally, since we do not have knowledge on the MWF and slow compartment integral (\(MWF\) and \(_\)), we set the corresponding elements of \(\lambda\) for these two parameters to negligibly small, non-zero values (see Tab. 1).

Accordingly, the above equations also apply for the inversion of magnitude ME-GE data, where the model vector \(x\) contains the respective six parameters of the \(_^\) distribution. Incorporating the complex signal model from Eq. (5), the previous equation naturally extend for the inversion of complex ME-GE data:

$$\mathop \limits_ \left\|S_ \left( t \right) - \hat_ \left( \right)\right\|^ + \left\|\lambda \cdot \left( } \right)\right\|^$$

(8)

with the parameterized complex signal forward model

$$}_\left(t|x\right)=\left[G\left(_^|_,_,\frac_\hspaceMWF}\right)}^_^+}\Delta _t}+G\left(_^|_,_,_\right)}^_^+}\Delta _t}\right]}^}_}$$

(9)

where the model vector \(x=\left(_,_,\Delta _,MWF,_,_,_,\Delta _,\Delta _\right)\) now contains nine values for the parameterization, including the MW and AEW frequency shifts, \(\Delta _\) and \(\Delta _\), and the constant phase term \(\Delta _\).

In contrast to the single inversions, the joint inversion aims to find simultaneously distributions for both \(_\) and \(_^\) which fit the respective decay curves, with the same \(MWF\) value in both distributions. The model vector \(x\) has now 14 parameters since it represents two signals, the SE signal modeled by six parameters and the GE signal modeled by nine parameters, but the \(MWF\) parameter is the same in both. Then, the minimization problem of the joint inversion can be stated as the sum of both residuals with a regularization term (analog to the single inversions) using 14 parameters:

$$\mathop \limits_ \,\,\,\alpha\left\| S_ \left( } \right) - \hat_ \left( \right)\right\|^ + \left\|S_ \left( } \right) - \hat_ \left( \right)\right\|^ + \left\|\lambda \cdot \left( } \right)\right\|^$$

(10)

where \(_\) denote the echo times of the ME-SE and ME-GE acquisitions, respectively, and \(\alpha\) is a scalar weighting factor to adjust the influence of the first term, the ME-SE data consistency. Since ME-GE acquisition has the double amount of data (magnitude and phase) than the ME-SE acquisitions (magnitude only), \(\alpha =2\) was chosen for equal weighting of both terms for all simulations and in vivo experiments.