High-Dimensional MR Spatiospectral Imaging by Integrating Physics-Based Modeling and Data-Driven Machine Learning: Current progress and future directions

A. Chemical Shift and MRSI

In a typical MRI experiment, the object of interest is commonly represented as a spatial function $\rho \left(r\right)$, whose values at a spatial location $r$ depend on the abundance of water molecules (or often referred to as spin density) and a few other biophysical parameters of tissue water, e.g., $T_1$, $T_2$, and $\chi$ (susceptibility). A detailed discussion on this can be found in [15]. However, there can be many different molecules other than water present in the object. Furthermore, one molecule can have the same nucleus (e.g., ¹H’s) in different functional groups, which are surrounded by varying numbers of orbiting electrons. These orbiting electrons “locally” perturb the magnetic field experienced by the nuclei. This effect, called electron shielding, illustrated in Fig. 1(a), makes different nuclei in the same or different molecules resonate at a range of frequencies specified by

where $B_0$ is the strength of the main magnetic field that the object experiences, $\sout{\gamma }$ is the gyromagnetic ratio (only nucleus dependent) and $\sigma$ is the electron shielding constant (dependent on nuclei, molecules as well as the position of a nucleus in a molecule). This frequency dispersion gives rise to the chemical shift phenomenon fundamental to MR spectroscopy. As an example, different protons in N-acetylaspartate (NAA, one of the most abundant metabolites in the brain) exhibit several frequencies, producing a unique resonance structure (Fig. 1(b)). Different molecules have their unique resonance structures.

Considering the chemical shift distribution in each imaging voxel, the image function of interest becomes a spatiospectral function $\rho \left(r,f\right)$. Mapping this function allows one to obtain spatially-resolved spectra, from which we can measure various molecules and quantify their relative abundance, providing useful insight into the physiological conditions of the tissues/organs of interest.

In practical MRSI acquisition with pulse excitations, we use a radiofrequency (RF) coil to pick up time-domain signals called free induction decays (FIDs). Therefore, we often use the Fourier counterpart of $\rho \left(r,f\right)$, i.e., $\tilde{\rho }\left(r,t\right)$, to represent the spatiotemporal function of interest (where $t$ denotes the FID dimension, sometimes referred to as the spin clock). Thus, the imaging data collected is in $\left(k,t\right)$-space, and related to $\tilde{\rho }\left(r,t\right)$ by

where $V$ denotes the imaging volume of interest, $\Delta f\left(r\right)$ denotes the macroscopic $B_0$ field inhomogeneity distribution (due to system imperfection and object induced magnetic field perturbation), $k_p=\left(p_x\Delta k_x,p_y\Delta k_y,p_z\Delta k_z\right)$ and $t_q=q\Delta t$ denote the sampling location in the high-dimensional $\left(k,t\right)$-space respectively ($\Delta k_x$, $\Delta k_y$, $\Delta k_z$ and $\Delta t$ are the corresponding sampling intervals), and $ϵ(\cdot )$ represent the measurement noise (modeled as complex white Gaussian). Based on Eq. (2), the MR spectroscopic imaging problem is to recover $\tilde{\rho }\left(r,t\right)$ from a set of $\left(k,t\right)$-space measurement ${\left\{s\left(k_p,t_q\right)\right\}}_{p=1,q=1}^{P,Q}$, also known as the spatiospectral encodings. Coil sensitivity can be included in Eq. (2) such that $s(\cdot )$ becomes $s_c(\cdot )$. Without loss of generality, we ignore the coil index in the following discussion. An illustration of the spatiospectral imaging problem is shown in Fig. 2.

B. Challenges for MRSI

As shown in Eq. (2) and Fig. 2, MRSI is a high-dimensional imaging problem due to the need to encode and decode both spatial and spectral information. Because $\Delta t$ needs to be small enough to satisfy the Nyquist sampling requirement along the spectral dimension (i.e., 1 to 2 kHz bandwidth for 3T and even higher for ultrahigh-field systems), it is not realistic to cover extended $k$-space during each readout. Therefore, a standard approach is to acquire all the time points needed for a single or a few $k$-space locations after each excitation and repeat the process many times. Furthermore, many FID $\left(t\right)$ points (on the order of hundreds) need to be acquired to achieve sufficient spectral resolutions for accurately differentiating signals from different molecules. As a result, the imaging time for an MRSI experiment is long and grows exponentially as the desired spatiospectral resolution increases in the conventional Fourier imaging paradigm.

Another fundamental challenge with MRSI is low SNR. The MR signal from each molecule is proportional to the bulk magnetization $M_{z,m}^0$ given by [15]

where $N_m$ is the number of polarized spins that can be used to generate detectable signals for the $m$th molecules, $\hslash$ is the Planck constant, $K$ the Boltzmann constant and $T_s$ the absolute temperature. As $N_m$ for the molecules of interest in MRSI, e.g., metabolites such as NAA, creatine (Cr) and choline (Cho) and neurotransmitters such as glutamate (Glu) and GABA, are typically three orders of magnitude smaller than that for water, much weaker signals from these molecules are expected. As a result, a common practice is to prescribe large voxel sizes (poor spatial resolution) in order to maintain sufficient SNRs within a clinically feasible scan time. A related important issue (which is often less discussed) is the large dynamic range for signals from different molecular components due to their concentration differences. An outstanding example is the strong nuisance water and subcutaneous lipid signals in ¹H-MRSI, which can significantly affect our ability to reliably reconstruct and quantify metabolite signals of interest. Therefore, strong prior information is needed to address the dimensionality, SNR and dynamic range challenges to enable fast, high-resolution, high-SNR MRSI for practical applications. The well-established chemical and physical knowledge for the spectroscopic signals and advanced data-driven machine learning tools to leverage this knowledge present new opportunities to tackle these challenges.

A. Subspace Model

Under physiological conditions, the parameters $c_m,T_{2,m}^{\text{*}}\left(r\right)$ and $\delta f_m\left(r\right)$ have values in a narrow range (e.g., $T_{2,m}^{\text{*}}~30$ to 100ms), and $M$ is usually small, i.e., 10–15 for in vivo ¹H-MRSI. Therefore, the physiologically meaningful high-dimensional FIDs/spectra we acquired should reside in or close to a low-dimensional space. One approach to exploiting this low-dimensionality is to approximate the high-dimensional $\tilde{\rho }\left(r,t\right)$ using the partial separability (PS) model [7], [8], [17]:

which represent the FIDs/spectra at individual voxels as linear combinations of a small set of basis functions $\left\{v_l\left(t\right)\right\}$, thus in a low-dimensional subspace. The PS model leads to low-rank structures for the discretized representation of $\tilde{\rho }\left(r,t\right)$, which has been used for denoising [18] or super-resolution reconstruction [19]. Such a low-rank filtering approach requires joint estimation of both the spatial coefficients $\left\{u_l\left(r\right)\right\}$ and subspace basis $\left\{v_l\left(t\right)\right\}$ (subspace pursuit). Subspace pursuit often requires relatively high SNR data which is well beyond that of in vivo MRSI data acquired using accelerated acquisition schemes. Alternatively, we can leverage physics-based priors to pre-learn $\left\{v_l\right\}$ from high-SNR training data and transform the MRSI problem into the recovery of $\left\{u_l\right\}$ which is significantly lower-dimensional than $\tilde{\rho }\left(r,t\right)$, thus enabling more flexible designs of data acquisition and physics-motivated low-rank reconstruction to address the speed, resolution and SNR (SRS) challenges. We will focus on this type of reconstruction methods in this review and more detailed discussions on data acquisition can be found in [20]–[22].

While subspace learning has been well studied in machine learning [23], the unique physical properties of each FID allows for unique training data generation and learning strategies [8]. More specifically, we can draw samples from pre-specified distributions of $c_m,T_{2,m}^{\text{*}}$ and $\delta f_m$ and synthesize a large number of FIDs from which the subspace structure can be determined, e.g., through PCA, ICA or other dimensionality reduction techniques. The parameter distributions can be specified using literature data with empirical distribution assumptions (e.g.,Gaussian distributions with lower and upper bounds) [24], [25]. In addition, spectral fitting using Eq. (5) or its variants can be applied to experimentally acquired high-SNR training data to extract empirical distributions of the molecular parameters [8], [12] or to supplement the synthetic data. The subspace learned this way incorporates strong physics prior, captures the inherent molecular spectral features of interest, and is less susceptible to experimental artifacts than a direct application of dimensionality reduction to noisy experimental data.

B. Nonlinear Manifold Model

While linear subspace models have been demonstrated highly effective and enables flexible sampling of $\left(k,t\right)$-space for the determination of spatial coefficients. Their effectiveness degrades when the range of spectral parameters increases, which motivates more general models. This is illustrated in Fig. 3. We generated spectra of glutamate (Glu) with varying ranges of $c$, $T_2^{\text{*}}$ and $\delta f$, and compared the representation accuracy between linear subspace and a nonlinear manifold model [26]. As can be seen, as the ranges of $T_2^{\text{*}}$ and $\delta f$ increase, the approximation errors of subspace increases, as shown by the error spectra in black, while the nonlinear model maintains high accuracy for the same order. Therefore, learning a nonlinear low-dimensional model can offer a more efficient representation of general in vivo MR spectra. As large quantities of training samples can be generated using the physics-based method described above, data-driven learning of nonlinear low-dimensional representations is possible, leveraging recent progresses in deep learning.

Nonlinear dimensionality reduction has been used for classifying single-voxel spectroscopy (SVS) data and demonstrated superior performance than linear dimensionality reduction [27], [28]. But these works did not focus on representation accuracy and did not consider the use of low-dimensional representations for image reconstruction. A deep autoencoder (DAE) based approach was proposed recently to learn an accurate low-dimensional representation for ensembles of spectroscopic signals (FIDs) [26]. Denoting each voxel FID as a vector $\rho$, an encoder $E\left(\cdot ;{\theta }_e\right)$ and a decoder $D\left(\cdot ;{\theta }_d\right)$ network can be learned such that a low-dimensional embedding encoded from $E(\cdot )$ can accurately reconstruct $\rho$, i.e., $\rho \approx$$D\left(E\left(\rho ;{\theta }_e\right);{\theta }_d\right)$. The idea has been extended to multi-echo-time (multi-TE) spectroscopy data [29], and the improved representation efficiency of learned nonlinear low-dimensional models over linear subspace model has been shown, for not only ¹H but also other nuclei data [25], [26], [29]. Figure 4 provides a conceptual illustration of physics-model-based training data generation incorporating QM priors and empirical distributions of spectral parameters, and data-driven learning of low-dimensional representations for high-dimensional spectroscopic signals.

These learned low-dimensional models enabled new ways to formulate the recovery of the high-dimensional image function (a.k.a. spatiospectral reconstruction) from noisy or often sparsely sampled spatiospectral encodings, which are reviewed in the next sections.

A. Subspace-Based Reconstruction

Using the discretized representation, i.e., representing $\tilde{\rho }\left(r,t\right)$ at a set of sampled locations ${\left\{r_n,t_{q^{\prime }}\right\}}_{n,q^{\prime }=1}^{N,Q^{\prime }}$, the subspace model in Eq. (6) can be rewritten in a matrix form $\rho =U\hat{V}$ [17], where $U\in C^{N\times L}$ and $\hat{V}\in C^{L\times Q^{\prime }}$ are matrix notations for the spatial coefficients and learned subspace. Accordingly, the reconstruction problem can be formulated as estimating the unknown spatial coefficients (with a significantly less number of degrees-of-freedom than) [17]

where $B$ captures the $B_0$ field inhomogeneity induced linear phases (i.e., $B_{n{q}^{\prime }}=e^{i2\pi \Delta f\left(r_n\right)t_q^{\text{'}}}$), $F$ is a spatial Fourier transform operator and $\text{Ω}$ denotes a $\left(k,t\right)$-space sampling operator. Note that $q^{\prime }$ can be different from $q$ as truncated FID sampling can be easily implemented. The vector $d$ contains all the noisy spatiospectral encodings. $R(\cdot )$ is a regularization that can be used to impose spatial constraint on $U$ or the images at individual FID time point, e.g., an edge preserving penalty [20] or sparsity constraint [30], and $\lambda$ is the regularization parameter. The final reconstruction $\hat{\rho }$ can be formed as $\hat{\rho }=\hat{U}\hat{V}$. This joint subspace and spatial constrained reconstruction has achieved impressive improvement in resolution, SNR and speed for both ¹H and X-nuclei MRSI in brain and muscle applications over traditional methods [9], [20], [31].

It is possible to extend Eq. (7) to jointly update both $U$ and $V$ [30], which is a special form of low-rank matrix recovery. This type of generic low-rank-model-based methods have been successful in several MRI applications, e.g., dynamic imaging [32], [33], MR relaxometry [34] and hyperpolarized ¹³C-MRSI [35], where SNR is not a major bottleneck. However, they do not work well when using MRSI data with very low SNRs for subspace pursuit, and are not flexible for incorporatingfuture development physics-motivated subspace learning. Substantial bias may be present in the reconstructed spectra due to errors of subspace estimation from very noisy data, especially for less dominant but important spatiospectral features [20]. A comparison of reconstructions from simulated MRSI data produced by direct low-rank filtering and using a learned subspace is shown in Fig. 5 to demonstrate this effect. Locally low-rank models offer better representation accuracy and capability than global low-rank models [33], [34], [36], [37], but are still limited by the SNR challenge and subspace estimation errors in MRSI applications.

B. Nonlinear Manifold Constrained Reconstruction

As discussed above, generalizing the linear subspace model to nonlinear low-dimensional manifold models can enable more accurate approximation of general spectroscopic signal variations, thus reducing potential modeling errors for diverse physiological and pathological conditions. A key issue is how to integrate such a learned nonlinear model into the reconstruction formulation. One approach is to incorporate a “network representation error” penalty term into a regularized reconstruction formalism:

with $𝒞(\cdot )$ representing the learned DAE $D\left(E\left(\cdot ;{\hat{\theta }}_e\right);{\hat{\theta }}_d\right)$. While the first term imposes data consistency same as previous formulations, the second term enforces the prior that the desired FIDs of interest at each voxel, ${\rho }_n$, should yield small representation errors for the learned low-dimensional model captured by $𝒞$. If $𝒞$ reduces to a linear network, this regularization can be viewed as penalizing the error for projection onto a learned subspace, a softened version of the explicit subspace model in Eqs. (7) – (12). $R(\cdot )$ is an additional regularization term that can be flexibly incorporated and used to impose spatial constraints, either hand-crafted or learned (see the next section for more discussion). Similar to the reconstruction strategies discussed above for the learned subspace, because we kept the forward encoding model and the learned $𝒞$ is imposing a prior on the underlying true FIDs, once pre-trained, $C$ can be used synergically with different acquisition designs. This is an important difference with another popular unrolling approach, in which an unrolled network is motivated from but does not exactly solve a regularized least-squares problem (the regularization term can not be explicitly defined) and the network typically needs to be retrained for different acquisition schemes. However, the disadvantage of actually solving the optimization problem in Eq. (8) is the need to invert the deep network (involving backpropagation) for each reconstruction, thus higher computational cost compared to the trained end-toend mapping network that only requires forward pass at the inference stage. Figure 6 provides an example to demonstrate the superior denoising reconstruction performance by using the learned nonlinear model (over linear subspace) for an in vivo ¹H-MRSI data.

To effectively take advantage of the learned nonlinear model without significantly compromising computational efficiency, projected gradient descent (PGD) based algorithms can be considered. For example, the recently proposed RAIISE method (LeaRning nonlineAr representatIon and projectIon for faSt constrained MRSI rEconstruction) [38] seeks to learn a projection network to extract low-dimensional embedding from high-dimensional, noisy FIDs and use it in an accelerated PGD algorithm. More specifically, the network-constrained reconstruction can be formulated as

where $\rho \in D\left(Z;{\theta }_d\right)$ enforces the prior that the underlying signals should yield a low-dimensional representation $Z$ (residing on a low-dimensional manifold) and $D\left(\cdot ;{\theta }_d\right)$ is the decoder from the representation network $𝒞$. The PGD update step can be realized via a projection operator $\text{Proj}(\cdot ):=D\left(P\left(\cdot ;{\theta }_p\right);{\theta }_d\right)$, where $P(\cdot )$ is a learned projector that recovers the latent representations from noisy FIDs and trained as follows:

where $x_j$ and ${\tilde{x}}_j$ are individual training FIDs and their noisy counterparts, respectively. A range of SNRs were considered to generate ${\tilde{x}}_j$ for training $P\left(\cdot ;{\theta }_p\right)$ parameterized by ${\theta }_p,{ϵ}_1$ and ${ϵ}_2$ assess the “projection” errors for the low-dimensional features as well as the full signals, respectively, and $\lambda$ balances the two losses. The encoder $E(\cdot )$ and decoder D(.)are from the representation network discussed above.The projector can be trained with different structures adapted to the data characteristics (single or multi-TE FIDs), and ${ϵ}_1$ and ${ϵ}_2$ can be chosen separately. It has been demonstrated that leveraging GPU acceleration RAIISE achieved similar or slightly better denoising performance for both single-TE and multi-TE MRSI with dramatically improved computation time compared to directly minimizing ${\Vert 𝒞\left({\rho }_n\right)-{\rho }_n\Vert }_2^2$ in the regularization formulation of Eq. (8) [38].

C. Integrating Learned Spatial and Spectral Priors

While the effectiveness of learned spectral priors incorporating both physics modeling and machine learning have been well documented, the investigations into learned spatial priors have been scarce. One of the major challenges is the difficulty in acquiring high-resolution, high-SNR MRS images as the “gold standard”, for either supervised training as commonly done in deep-learning-based MRI reconstruction or self-supervised representation learning. One solution is to simulate spatial distributions of spectra by combining anatomical images and metabolite concentrations and lineshapes from tissue/brain-regionspecific literature values and/or experimental data. For example, a dense U-Net trained completely by synthetic data was used to produce high-resolution metabolite maps from the low-resolution counterpart plus a T1w MRI [10]. The main issue, however, is that the synthetic spatial distributions of metabolites may not capture sufficient variations in real data, especially those related to different pathological conditions.

An alternative approach is to construct generic network-based image representations that do not require high-SNR, high-resolution data for supervised training of image mapping. Recently, efforts have been made to combine subspace constraint and generative network based image representation. The work in [39] used deep image prior (DIP) to model the spatial function $U$ as $U=f\left(Z;\theta \right)$, where $f\left(\cdot ;\theta \right)$ is a network with trainable parameters $\theta$, the reconstruction in Eq. (7) is then changed to

With a chosen network architecture (e.g., U-net in [39]), the parameters $\theta$ are estimated from noisy data with $Z$ being an anatomical image. This method assumes that the network can account for the contrast difference between the desired spatial coefficient maps and an MR image to leverage high-resolution anatomical prior information. However, DIP can overfit noise by allowing all the network parameters to be updated and anatomical image being the input, the network may introduce undesirable bias into the spatial reconstruction. Thus, choosing a proper early stopping criterion to balance the two factors may be challenging [40]. Another strategy is to leverage a generative network model pretrained using anatomical, multicontrast MRIs with sufficient representation power for images with different contrasts, and solve for the latent variable $Z$ instead of network parameters $\theta$ during reconstruction. As $Z$ is usually lower dimensional than the voxel-wise image representation, noise reduction is possible. For example, a GAN-regularized subspace reconstruction can be formulated as [41]

where $f\left(Z;\theta \right)$ is a pretrained GAN (e.g., using T1 and T2w images). This problem can be solved using alternating minimization and the step for updating $Z$ is a GAN inversion problem, which is an active research topic at the intersection of inverse problem and deep learning. While different GAN designs can be considered for $f\left(\cdot ;\theta \right)$, StyleGAN combined with intermediate layer optimization (ILO) has been shown to achieve the best inversion performance [42]. A more comprehensive review on this topic can be found in [43]. One important feature of these reconstruction methods is the integration of the physical spatiospectral encoding model and deep learning priors, in contrast to training an end-to-end network for direct mapping from data/noisy images to reconstructions. As a result, the pretrained model can work for any sampling pattern, sequence parameters, and different resolution and SNRs, without the need of retraining for each case. Meanwhile, the computation time required for the network parameter update or GAN inversion is one limitation/consideration.

A. Low-Rank Tensor Models with Learned Subspaces

To tackle the increased dimensionality problem, the PS model in Eq. (6) can be generalized to represent new image functions of interest $\tilde{\rho }\left(r,t,T\right)$, with $T$ denoting the additional time dimension $\left(T=1,2,\dots ,N_T\right)$: [9],[49]

where $N_T$ is the total number of time frame, ${\left\{g_s\left(T\right)\right\}}_{s=1}^S$ denote the basis functions that span the space of temporal variations and ${\left\{c_{m,l,s}\right\}}_{m,l,s=1}^{M,L,S}$ are the model coefficients. Eq. (13) implies a low-rank tensor structure (generalized from the low-rank matrices) that exploits the spatial-spectral-temporal correlations in dynamic MRSI data with $\left\{c_{m,l,s}\right\}$ referred to as the core tensor [49], [50]. $M$, $L$ and $S$ are the respective model orders (typically much smaller than $N$ and $Q$). Similarly as the spectral basis $\left\{v_l\left(t\right)\right\}$ can be learned by incorporating physics-based model, the temporal basis $\left\{g_s\left(T\right)\right\}$ can be predetermined from temporal training data [9], leaving only $\left\{u_m\left(r\right)\right\}$ and $\left\{c_{m,l,s}\right\}$ as the unknown in the imaging problem. Eq. (13) can be further adapted into a molecule-dependent tensor model [51]

where $M$ is the number of molecules of interest, $L_m$ and $S_m$ are the spectral and temporal model orders of the $m$th molecule, $c_{m,l,s}\left(r\right)$ combines the core tensor coefficients and molecule-dependent spatial basis. This model offers up to two orders-of-magnitude reduction in the degrees of freedom (DOF) [51] compared to the canonical truncated Fourier series representation, thus significantly improved resolution and SNR tradeoffs. Furthermore, the model in Eq. (14) provides the flexibility to incorporate molecule-specific physics- and biochemistry-based priors [44]. The spectral basis can be obtained as described in section III.A and the temporal basis can be derived by exploiting existing physiological and metabolic knowledge. More specifically, the temporal variations of the spatially-resolved spectra are often results of metabolic process in living organs, which can be approximated by a set of parametric curves specified by some time constants [51], e.g., the exponential recovery of PCr after ischemia or exercise. The variable $s_m\left(r\right)$ can be used to impose molecule-specific spatial support and distribution priors from high-resolution reference images available, further improving the performance (see [51] for more details). Accordingly, the dynamic MRSI reconstruction can be formulated as follows:

where $𝒜$ describes the forward model as discussed above (capturing field inhomogeneity effects and spatiospectral encoding process), $C_m$, $V_m$, $G_m$ and $s_m$ are tensor forms of $\left\{c_{m,l,s}\left(r\right)\right\}$, $\left\{v_{m,l}\left(t\right)\right\}$, $\left\{g_{m,s}\left(T\right)\right\}$ and $\left\{s_m\left(r\right)\right\}$ respectively, where $C_m\in C^{N\times L_m\times S_m}$, $V_m\in C^{L_m\times Q^{\prime }}$, $G_m\in C^{S_m\times N_T}$ and $s_m\in C^{N\times 1}$. $\odot$ and $\otimes$ denote column-wise and tensor products. The algorithm in [51] first performed a series of reconstruction of time-dependent spatial coefficients (using per-estimated $V_m$), i.e., $\left\{{\hat{c}}_{m,l}\left(T\right)\right\}$, and then fit the temporal variations to biochemical models for each molecule to determine $G_m$. Finally, Eq. (15) was solved using estimated spectral and temporal basis. $R\left(\left\{C_m\right\}\right)$ imposes regularization to encourage spatial smoothness or sparsity.

Considering the unique properties of the dynamic spectroscopic signals, Eq. (14) can be reformulated by imposing strict separability between the spatial-spectral and spatial-temporal variations [48], i.e.,

where $a_{m,l}$ and $b_{m,s}$ are the new spatial coefficients for the spectral and temporal bases, respectively. With this model, the number of unknowns at each voxel is reduced from ${\sum }_{m=1}^M{L}_m\times S_m$ to ${\sum }_{m=1}^M\left(L_m+S_m\right)$. This does not necessarily increase modeling error for dynamic MRSI specifically, because it is accurate to assume all FID/frequency points corresponding to the same molecule should share the same dynamics at each voxel. Using matrix notations for $a_{m,l}$ and $b_{m,s}$, a similar regularized least-squares problem to Eq. (15) can be formulated.

B. Combining Low-Rank Tensor and Manifold Modeling

As in the “static” MRSI case, it is straightforward to incorporate nonlinear manifold constraints on spectral variations into the reconstruction problems in Eqs. (15) or (16). For example, the learned network can be used to regularize individual FID’s at every time point, i.e.,

This formulation assumes that the desired FIDs at each voxel and each time point (i.e., ${\rho }_{n,T}$) reside on the same learned manifolds and can be solved using alternating minimization. Note that the ensemble temporal variations are still captured by the linear subspace model. To further exploit the prior information available along individual dimensions, the subspace model for temporal variations can be generalized to nonlinear manifolds as well. To this end, the separability representation in Eq. (16) offers additional flexibility to learn and incorporate low-dimensional manifold constraints on molecule-specific temporal dynamics. Specifically, a set of molecule-dependent DAEs, denoted as ${𝒢}_m(\cdot )$, were trained (using both synthetic and experimental dynamic data) and integrated into the following low-rank-tensor-based reconstruction problem [48]:

where $A_m$ and $B_m$ are matrix forms of $\left\{a_{m,l}\left(r\right)\right\}$ and $\left\{b_{m,s}\left(r\right)\right\}$ with $A_m\in C^{N\times L_m}$ and $B_m\in C^{N\times S_m}$. $\circ$ denotes the operation defined in Eq. (16), and ${\lambda }_1$ and ${\lambda }_2$ are regularization parameters controlling the trade-off between data consistency and manifold (spectral and temporal) representation errors. Additional regularization terms on the spatial distribution of $a_{m,l}$ and $b_{m,s}$ may also be included. The problem in Eq. (18) can be solved via alternating minimization of $A_m$ and $B_m$, which also involves backpropagation of the learned networks. Details of the exact algorithm can be found in [48]. While parallelizable, the reconstruction can be computationally demanding due to the need to invert multiple networks for signals at many time points. The formulation used in RAIISE described in Section IV.B may be helpful in addressing this issue. A set of denoising dynamic ²H-MRSI results is shown in Fig. 8 to illustrate the capability of the low-rank-tensor-model-based reconstruction using learned spectral and temporal subspaces. Comparisons of estimated metabolite maps, spatially-resolved spectra and temporal dynamics from standard Fourier reconstruction, direct low-rank tensor truncation and learned-subspace-based methods were made.

C. The Motion Issue

The formulations discussed above do not fully consider motions (involuntary or voluntary) during the time window when the dynamic metabolic activities are monitored. This may be an issue when the imaging time is long (e.g., more than just a few minutes) and/or the organ of interest is more susceptible to motion. Prospective and retrospective motion tracking and correction may be used to mitigate its effects [22], [52]. Another possibility is to incorporate motion models into the $T$ dimension, which might lead to an increased model order and degradation in the accuracy of capturing metabolite dynamics. A new dimension can be introduced to represent different motion phases for motion-resolved dynamic MRSI using a higher-order tensor model [50].