Joint multi-field T_1 quantification for fast field-cycling MRI
Markus Bödenler, Oliver Maier, Rudolf Stollberger, Lionel M. Broche, P. James Ross, Mary-Joan MacLeod, Hermann Scharfetter
FF U L L P A P E R
M a g n e t i c Re s o n a n c e i n M e d i c i n e
Joint multi-field T quantification for fast field-cycling MRI Markus Bödenler | Oliver Maier |Rudolf Stollberger | Lionel M. Broche |P. James Ross | Mary-Joan MacLeod |Hermann Scharfetter Institute of Medical Engineering, GrazUniversity of Technology, Graz, Austria Institute of eHealth, University of AppliedSciences FH JOANNEUM, Graz, Austria Aberdeen Biomedical Imaging Centre,University of Aberdeen, Foresterhill, AB252ZD, Aberdeen, UK Institute of Medical Sciences, University ofAberdeen, Foresterhill, AB25 2ZD,Aberdeen, UK BioTechMed-Graz, Mozartgasse 12/II,A-8010 Graz, Austria
Correspondence
Markus Bödenler, Institute of eHealth,University of Applied Sciences FHJOANNEUM, Graz, AustriaEmail: [email protected]
Funding information
ÖAW DOC Fellowship 24966
WORD COUNT:
Approximately 4900
Recent developments in hardware design enable the use ofFast Field-Cycling (FFC) techniques in MRI to exploit thedifferent relaxation rates at very low field strength, achiev-ing novel contrast. The method opens new avenues forin vivo characterisations of pathologies but at the expenseof longer acquisition times. To mitigate this we propose amodel-based reconstruction method that fully exploits thehigh information redundancy offered by FFC methods. Thisis based on joining spatial information from all fields basedon TGV regularization. The algorithm was tested on brainstroke images, both simulated and acquired from FFC pa-tients scans using an FFC spin echo sequences. The re-sults are compared to non-linear least squares combinedwith k-space filtering. The proposed method shows excel-lent abilities to remove noise while maintaining sharp imagefeatures with large SNR gains at low-field images, clearlyoutperforming the reference approach. Especially patientdata shows huge improvements in visual appearance overall fields. The proposed reconstruction technique largelyimproves FFC image quality, further pushing this new tech-nology towards clinical standards. * Equally contributing authors. a r X i v : . [ ee ss . SP ] F e b K E Y W O R D S fast field-cycling, dispersion, T quantification, model-basedreconstruction, low-field MRI | INTRODUCTION
The magnetic field dependency of the longitudinal and transverse relaxation times, also referred to as Nuclear Mag-netic Relaxation Dispersion (NMRD), provides insight into the underlying structural order and dynamics of a widerange of molecular systems [1, 2]. In recent years, the T dispersion of protons in particular has experienced increasedinterest for the investigation of biomarkers related to various pathological processes [3, 4, 5, 6]. The field depen-dent properties of such biomarkers are invisible to traditional MRI scanners, which operate only at one fixed mainmagnetic field strength and are restricted to the measurement of relaxation times corresponding to the B field em-ployed. However, new MRI-derived technologies are emerging that allow exploring different magnetic fields withina single system. One such technology is Fast Field-Cycling Magnetic Resonance Imaging, also known as FFC-MRIor FFC imaging, which enables a modulation of the main magnetic field during an imaging sequence giving accessto field-dependent relaxation properties as a novel contrast mechanism [7]. FFC imaging derives from MRI but usesradically different technologies to generate the main magnetic field, and both types of scanner offer different viewson biological tissues.Indeed, varying the main magnetic field within a defined range requires dedicated hardware and various ap-proaches exist to realize FFC imaging systems [8]. In the clinical field range, FFC imaging is implemented by means of a B insert coil together with the superconducting magnet provided by a commercial MRI system for 1.5 T [9, 10, 11] or3 T [12]. This approach, also referred to as delta relaxation enhanced MR (dreMR), has auspicious applications for thedetection and quantification of contrast agents with increased specificity and sensitivity [9, 13, 14, 15]. Several sys-tems were also developed to access the endogenous T dispersion of tissues in the low-field regime [16, 17, 18, 19, 20].Recently, a whole-body FFC scanner approved for clinical imaging studies was reported, capable of reaching any fieldfrom 50 µ T to 0.2 T [21]. Controlled variations of the magnetic field with this single resistive magnet design allowfor multi-field T quantification over a wide range of field strength while retaining image quality down to ultra-lowfields. Pilot studies show promising potential for innovations in the imaging of osteoarthritis [4], sarcoma [22] or brainstroke [21], with potentially important applications in medicine as an in vivo assessment method of multi-field T and T dispersion information.Compared to conventional MRI, implementation of fast field-cycling poses additional demands on power supplies,control electronics, magnet design, pulse sequences and image quality [8]. Signal to noise ratio (SNR) is therefore animportant issue for FFC scanners to satisfy the latter. High fields benefit from an inherently high SNR as they rely onstable and homogeneous acquisition fields provided by superconducting magnets. Although the SNR is not a limitingfactor for the individual images, contrast in dreMR is obtained by image subtraction and strongly depends on the T dispersion of the contrast agent in use [23]. The magnitude of the dreMR signal is rather small in comparison to theindividual images (e.g. about 2.5 % in [12]) and retaining sufficient high SNR may become an issue. Similarly, low-fieldFFC imaging systems operate with acquisition fields of 0.2 T or less, which limits the SNR compared to conventionalclinical fields due to its dependency to B . Moreover, image quality deteriorates because of poor magnetic fieldhomogeneity, field instabilities during operation and delays in the field ramps between different phases in the pulsesequence especially for ultra-low evolution fields.For all these reasons, both high- and low-field FFC scanners may strongly benefit from SNR-enhancing methods and a vast number of these have been developed for high field MRI images in the recent years [24, 25, 26, 27, 28,29]. These can be divided into denoising and reconstruction-based approaches. The former takes a series of noisyinput images and tries to find a denoised solution by making use of a priori knowledge in the form of regularization.Regularization can utilizes either spatial information, information from the acquired series, or a combination of both.Denoising approaches are generally simpler to implement and computation time is lower compared to reconstructionapproaches, but they can not recover structures that were missed by the k-space to image space transformation dueto poor SNR. To this end, recent approaches rely on a constrained reconstruction process, incorporating the a prioriknowledge in the image generation process [26, 27, 28, 29]. In the case of quantitative MRI, this approach can betaken one step further by including the non-linear MRI signal model into the reconstruction process, thus, directlyacting on the parameter maps of interest [30, 31, 32, 33]. This kind of fitting approach is known as model-basedreconstruction in the high field MRI regime. Most regularization strategies rely on some sort of sparsifying transformto separate image content from noise and artefacts. Commonly used transforms include finite differences basedapproaches [27, 34, 35] and applications of the wavelet transformation [36, 37, 38]. The regularization functionalused highly influences the appearance of the final reconstruction and should be chosen based on a priori knowledgeabout the given parameter map. It was shown that Total Generalized Variation (TGV) [35] priors are a superb choicesfor both image reconstruction [27] and quantitative MRI [33], leading to high quality reconstruction results withoutthe stair-casing artefacts of total variation (TV) [27].These methods apply well to the estimation of T maps. Standard T quantification with inversion recovery se-quences requires the acquisition of an image series with different inversion times, leading to high redundancy in theinformation collected that can be exploited by the regularisation algorithm. FFC imaging adds an extra dimension toMRI by varying the magnetic field during the relaxation phase of the pulse sequence, thus providing an additional fieldseries. This multi-field data offers new possibilities to exploit information redundancies to improve the quantificationprocess. These redundancies could be utilized by a model-based reconstruction approach, incorporating the data fromall measurements at different field strengths into one large optimization problem. Each field leads to an individual T map which shares common information with the other fields, e.g. most edges in the T map should coincide. This in-formation can be exploited by joining the individual regularization functionals in parametric dimension via a Frobeniusnorm. Such an approach was shown to further improve the quality of the resulting parameter maps in the context of T mapping from highly subsampled data [33].Herein, we formulate the multi-field FFC imaging parameter quantification as a non-linear model-based recon-struction problem with Frobenius type TGV regularization. With this formulation as a single optimization problemit is possible to exploit all the joint spatial information of the additional field dimension to stabilize the quantificationprocess and hence enhance the image quality. The proposed method is evaluated on simulated numerical FFC imagingdata as well as on in vivo datasets from two stroke patients and compared to standard non field-dependent methods.The results show improved stability of the parameter quantification with excellent noise suppression properties. Inparticular, the proposed method reveals remarkable contrast between the lesion and surrounding tissues in case ofultra-low fields. | THEORY2.1 | FFC imaging signal model
In the most general case of a FFC imaging pulse sequence, the main magnetic field is rapidly cycled between threedifferent levels: polarization field B P , evolution field B E and signal detection field B D . A designated pre-polarization of the sample magnetization is not necessarily required in the high SNR regime and the polarization field can be setto the detection field, i.e., B P = B D . For simplicity we will assume that this is the case and we will refer to these fieldsas B and the corresponding equilibrium magnetization as M , respectively, as this does not alter the validity of ourapproach in the case of low-field systems since the effect of polarisation at a different field can be compensated bya polarisation efficiency term that blends into the inversion efficiency parameter. A schematic of a typical inversionrecovery FFC imaging pulse sequence can be seen in Figure 1: following an inversion RF pulse, the main magnetic fieldis cycled to the desired strength B E and the spin system undergoes a relaxation associated with the applied evolutionfield during a given evolution time t evo . The longitudinal magnetization M z at the end of this evolution period is givenby M z ( t evo ) = [− αM − M E ] e − t evo T E + M E , (1)where M and M E are the equilibrium magnetizations for the detection and evolution field, respectively. The T relaxation time, corresponding to the evolution field applied, is given by T E and α corrects M for field dependenteffects from ramping the field combined with non-ideal inversion efficiency of the RF pulse [39]. The equilibriummagnetization is proportional to the strength of the applied magnetic field, i.e., M = C B and M E = C B E for thedetection and evolution field, respectively. This can also be written as M B = M E B E = C . (2)After the evolution period the signal is acquired at B to ensure that the Larmor frequency of the spins corresponds tothe tuning frequency of the receive RF coil. With equations (1) and (2), and the unknown parameters u = ( C , α , T E ) ,the acquired signal S B E , t evo ( u ) can be modeled by the non-linear signal equation S : U → D , given by S B E , t evo ( u ) = F (cid:40) C [− α B e − tevoT E + B E ( − e − tevoT E ) ] (cid:41) , (3)with F representing the Fourier transformation and sampling of k-space. | Multi-field parameter fitting
Acquiring several time points t evo for a specific evolution field B E allows to quantify C , α , and T E . Typically, each B E field yields a different T E value, thus the fitting process must be repeated for each evolution field, omitting jointinformation in the unknowns such as structural information. By combining the separate fitting steps into a single op-timization problem it is possible to utilize the shared information to stabilize the quantification process. Furthermore,joint information between different parameter maps can be exploited by means of a Frobenius type functional. Sucha fitting and regularization strategy has been successfully applied in other multi-channel fitting problems [33, 40, 41].Thus we propose to apply a similar approach to quantify C , α , and T E from multi-field FFC imaging data, utilizingshared information, especially between the different T E maps. | Fixing notation
Throughout the course of this work we fix the following notations. The image dimensions in 2D are denoted as N x and N y , defining the image space U = (cid:195) N x × N y with p = ( x , y ) defining a point at location ( x , y ) ∈ (cid:206) . u ∈ U N u expresses the space of unknowns u = { C , α B E = ( α B E , α B E . . . , α B NE ) , T E = ( T E , T E , . . . , T E NE ) } ∈ U N u -maps,with N u = 1+2 N E and N E the number of evolution fields E . The measured data space is denoted as D = (cid:195) N kx × N ky × N d and is constituted of N d = (cid:205) N E i =1 N E i t measurements of 2D k-spaces N k x × N k y . Each measurement is a combinationof an evolution field B E i = ( B E , B E , . . . , B E NE ) ∈ (cid:210) N E + with an associated number of measurement time points t E i n = ( t E i , t E i , . . . , t E i N Eit ) ∈ (cid:210) N Eit + , and N E i t the number of time points acquired at a specific field strength E i . Tosimplify notation we will drop the indices and refer to α B E as α . | Model-based reconstruction framework
Assuming Gaussian noise corrupts the measurement data d , it is possible to quantify the unknown parameter u via aregularized non-linear, minimum-least-squares problemmin u (cid:107) A ( u ) − d (cid:107) + γR ( u ) , (4)which origins from a maximum a posteriori approach using Bayes’ theory. For multi-field FFC data u is linked tothe measurement data d = ( d , , d , , . . . , d , N E t , d , N E t , . . . , d N E , N ENEt ) ∈ U via S B Ei , t Ein : u ↦→ d i , n We denote allmeasured data as d and the corresponding mapping from unknown space to data space as S . Thus, the optimizationproblem is defined as min u , v (cid:107) S ( u ) − d (cid:107) + γ ( β (cid:107) (cid:43) u − v (cid:107) , , F + β (cid:107) E v (cid:107) , , F ) . (5) R ( u ) is chosen as TGV regularization with a joint Frobenius norm on all unknowns u . This type of regularizationwas shown to have favorable properties for multi-parameter model-based reconstruction [33]. The (cid:107) · (cid:107) , , F termsresemble the Frobenius type TGV functionals, joining common spatial information of the unknown parameter mapsby combining gradient information of all maps via an L -norm in parameter direction.Using the TGV model parameters β and β it is possible to balance the approximated first and second derivatives,avoiding the stair-casing artifacts of TV while maintaining its favorable edge-preserving features. The ratio β / β = 1 / of TGV is fixed throughout this work, as it was shown to yield good results for MRI image reconstruction [27]. Thenumerical solution is done in analogy to [33] via a Gauss-Newton (GN) approach. Leading to an inner GN problem ofthe form min u , v (cid:107) DS u − ˜ d k (cid:107) + γ k ( β (cid:107) (cid:43) u − v (cid:107) , , F + β | (cid:107) E v (cid:107) , , F ) + δ k (cid:107) u − u k (cid:107) M k . (6)The linearization is done via a Taylor series expansion of S w.r.t. each unknown in u at position u k . Constant terms arefused into ˜ d k to keep the notation clean. DS amounts to the Jacobian of S , evaluated at u k . Introducing an additional weighted L -norm penalty on u improves convexity of the function. The weighting matrix M k can be used to resemblea Levenberg-Marquart update if M k = di ag ( DS T DS ) .Using Fenchel duality the problem of non-differentiability can be overcome and equation 6 can be cast into asaddle-point form min u max y (cid:104) K u , y (cid:105) + G ( u ) − F ∗ ( y ) . (7)Problems in such a form can be solved via a primal-dual algorithm [42] using a line-search to speed up convergence [43].Mathematical details for each step are given in Appendix A. The update scheme written as pseudo code is given inAppendix B. | METHODS3.1 | Numerical FFC imaging data
To evaluate the proposed model-based reconstruction approach numerical FFC imaging data were simulated usingparameters measured from FFC imaging scans of brain stroke patients at the University of Aberdeen as part of aseparate study (PUFFINS study, see details in section 3.2). The numerical phantom followed a schematic geometryand dispersive characteristics of an axial head scan with four regions (see Figure 3), representing the subcutaneousfat (ROI 1), the tissues surrounding the brain (ROI 2), the brain (ROI 3) and a stroke-like lesion (ROI 4). T values weresimulated by means of a power-law dispersion, / T = aB E b , coarsely in line with proton T -NMRD profiles of fat(ROI 1), white (ROI 2) and grey matter (ROI 3) and stroke lesions (ROI 4) measured in vivo from the PUFFINS patientscohort (data to be published). The values retained for the different evolution fields and times are summarized in Table 1and Table 2, respectively, and served as a ground truth for the validation of the T quantification. The numerical FFCimaging phantom was first generated as a vector graphic to be subsequently converted to matrix data to allow forany desired sampling resolution. In this case we used an image resolution with matrix size of × pixel, whichis typical for the original FFC imaging of stroke patients. Tissue reference values were assigned for each ROI and adata series was generated using the signal equation in 3. Additionally, a small constant phase offset was introducedfor each α .Simulated proton density values were normalized to 1, resulting in a theoretical maximum signal amplitude of1 for the simulated series. Zero-mean Gaussian noise was added on both real and imaginary parts of the imagesto simulate the noise arising from the patient tissues and acquisition system. The noise amplitude was selected asa percentage of the theoretical maximum signal in the overall ground truth image series ranging from 1 % to 4 %,reflecting a typical range from the FFC images acquired. SNR directly after inversion at B ranged from 33.3 to 8.3 inthe white matter ROI and 66.7 to 16.7 for the grey matter ROI, respectively. Note that low-field FFC images tend toexhibit markedly lower signal strength than higher-field ones because of losses when the magnetic field switches, sotheir SNR was proportionally more affected using this approach. Finally, the image series was transformed to k-spacevia a 2D Fourier transformation, as input for the proposed method. The reference method was applied to image series. | In vivo FFC imaging data
The performance of the proposed method was tested on in vivo FFC imaging patient data. Two data sets obtainedfrom patients scanned for a brain stroke were selected, as part of the PUFFINS study currently taking place at the
University of Aberdeen. This study has been approved by the North of Scotland Research Ethics Committee (studynumber 16/NS/0136) and all the participants agreed for the FFC imaging data to be used anonymously for researchpurposes. The scans selected both present a lesion in the ultra-low field regime that could not be easily observed at200 mT, as illustrated in Figure 2 for patient I. Both measurements were performed using a whole-body FFC scanner[21] using a FFC inversion-recovery spin echo sequence [44] with an echo time of 24 ms, 20 kHz bandwidth, 8.37MHz acquisition frequency, 10 mm slide thickness and single slice acquisition. The images had a field of view of 290mm and a resolution of 128 x 128 pixel in-plane with 80 phase encode acquisitions and partial Fourier acquisition (80lines out of 128). The sample was pre-polarised at 200 mT for 300 ms before each evolution periods with the timingsas shown in Table 1, for an acquisition time of 40 min. | Data processing and corrections
The original raw image was reconstructed using partial Fourier completion to recover the correct image ratio. Phase-encode artefacts were removed using a method previously published [45] but the images had not been filtered orfurther modified. While the noisy images were used as input for the standard pixel-based fitting data, the correspond-ing noisy k-space data was used as input for the proposed fitting process.As reference method lsqnonlin of Matlab (The MathWorks, Inc.) was used for fitting equation 4, where R ( u ) wasreplaced with Tikhonov regularization on the unknowns to stabilize fitting. The method was applied to each fieldseparately. Prior to fitting, images were smoothed in k-space using the following filter function f ( k ) = 12 + 1 π arctan β k c − | k | k c , (8)with k c = 30 denoting the cutoff radius, k the k-space location, and β = 100 as parameter for the slope of the filter.The analyses with the proposed method were done by implementing the FFC signal model in PyQMRI [46]. Allfittings were performed on a desktop PC equipped with an Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz with 64gigabyte of RAM and a NVIDIA GeForce GTX 1080 Ti GPU with 12 gigabyte of RAM. | Optimization
The regularization weights γ k and δ k were reduced after each linearization step, following the iterative regularizedGauss Newton scheme [47]. γ k = 10 − and δ k = 1 were used as initial values and were reduced by a factor of . and . , respectively. To account for the typical smooth appearance of α , corresponding regularization weights weremultiplied by a factor of 10. The reduction steps were repeated down to γ min = 4 × − and δ min = 10 − . In total, 12linearization steps were performed. The number of primal-dual iterations for each sub-problem was doubled startingat 10 iterations up to 2000 iterations, i.e. i t er k = mi n ( ∗ k , ) . If the relative decrease in the primal problemor the decrease of the primal-dual gap was less than − , the inner iteration was terminated. The step sizes of theemployed primal-dual algorithm were determined via a line-search, described in Algorithm 2 [43]. | Code availability statement
The Code used for this publication is integrated in PyQMRI [46] and is available at https://github.com/IMTtugraz/PyQMRI . | RESULTS4.1 | Numerical FFC imaging data
The simulated high noise level can be seen as residual noise in the reconstructed T maps of the pixel-wise fittingapproach (figure 3). Simultaneously, a difference to the simulated reference is visually noticeable in the pixel-wisefitting approach. The proposed model-based method is able to reduce outliers throughout all noise levels and isvisually closer to the simulated reference values. Plots of C and α in Supporting Information Figure S1-S4 showsimilar results. The pixel-wise fitting approach even fails to capture the correct phase of the simulated phantom.A pixel-wise relative absolute difference plot (figure 4) confirms this visual impression of reduced noise using theproposed approach. The proposed method shows an up to 18 fold lower mean error in the phantom, computedover all pixels, than standard pixel-wise fitting. The error increases with increased noise level, as can be expected.Difference plots also reveal a slight bias of the proposed method. The bias of both methods is further assessed in 2Djoint histogram plots (figure 5). For these plots, T values of all fields are combined to form a single plot. The proposedmethod shows slight underestimation of high T values, as reflected by points lying below the identity line. However,noise could be greatly reduced compared to the pixel-wise fitting and different T ranges are clearly separated andshow a similar distribution as the simulated values. Fitting with the standard method took approximately 100 seconds.The proposed method took roughly 120 seconds. | In vivo FFC imaging data
The improvements in T estimations held true when processing real FFC imaging data from stroke patients. The T maps of unfiltered FFC images obtained using standard fitting-based processing methods could not resolve anatomicalfeatures inside the brain region, as seen in figures 6 and 7. The proposed method offers clear distinguishable structuresin T maps at 200 mT and is even able to recover some structural details in lower fields. It also assessed sharp featuresaround the lesion area appearing at 37 mT and below in both patients.The quality of the T maps obtained allowed estimating the T dispersion curves for different Regions Of Interest(ROI), as shown on Figure 8 for subcutaneous fat selected under the scalp, the area of the lesion observed at thelowest field strength, and white and grey matter as seen at the highest field strength (the ROIs are shown in Sup-porting Information Figure S5). The dispersion profiles of fatty tissues show large standard deviations, which may beattributed to the presence of various types of tissues within these ROIs, due to the relatively low resolution of theimage. It can also be noted that the relaxation time of the lesions were comparable to those of normal brain tissuesat 200 mT but changed markedly at lower magnetic fields. Fitting took approximately 65 and 150 seconds with thestandard method for patient I and II, respectively. The proposed method took 100 and 240 seconds for each patient,respectively. | DISCUSSION
The approach used here has high potential to serve as a new standard procedure for fast post-processing of FFC MRIdata. As the phantom simulations showed, the noise in the reconstructed T maps could be reduced very efficientlywhile preserving important anatomical details to a high extent. The algorithm outperforms established methods basedon pixel-wise fitting of the relaxation profiles yielding lower deviations from the reference values and significantly lessvariance (figure 3 , 4, and 5). The improved stability results from the combination of information from all acquired fields and exploiting the existence of similar topological structures in the different unknowns. The improved stability is alsoreflected by increased accuracy of recovered pseudo proton density C and correction factor α values (Supporting In-formation Figure S1-S4). Also the T maps show significantly reduced variance though there remains some bias whichmay be, at least partially, due to residual errors in α . Another remarkable feature of the proposed method is its abilityto accurately recover the phase information in C and α , making phase correction prior to fitting obsolete. This in turncan improve T maps as no normalization with a noise phase estimate is necessary.The advantages of the improved fitting approach become immanent in the in vivo applications (figure 6 and 7). Thestandard approach fails to reconstruct image details in both patients. In current practice, k-space windowing filtersare applied to recover usable information but this dramatically reduces image resolution by filtering out the high-frequency components of the image, which are responsible for the sharp features. In contrast, the joint regularizationapproach can recover clearly distinguishable grey and white matter regions at 200 mT on the two patient datasets,previously hidden in noise. The values obtained for the different regions of interest agree well between the patients,given the estimation of the error provided by the variation of the T values within each ROI (figure 8). The T valueswere systematically higher in patient I than in patient II, which may be attributed to patient variability and differentRF receive coil sensitivity relative to the used ROIs.As in the phantom images the proposed method proved to preserve anatomical features with high spatial frequen-cies because of using the existence of sharp edges for regularization. This approach is increasingly accurate with thenumber of views that can be compared showing the same object, either as a repetition of a recording or as differentacquisition of the same field of view, as it is the case here. Hence FFC imaging can benefit from the high informationredundancy obtained from the typical acquisition method, which repeats the measurement of the field of view atdifferent evolution times and fields.As the proposed approach is also model-based, and can therefore provide T directly from the raw images, it couldbe used to reduce the number of steps required to process the image and limits data losses. However model-basedapproaches also limit the amount of information that is extracted from the image, and properties not covered by thesignal equation, may be missed. For instance, brain tissues are known to follow bi-exponential relaxation because ofthe presence of the myelin sheath around the axons. Hence in a subsequent step, the model will be adapted to thetype of scan, or following a test for potential multi-exponential behaviour [48].Using direct reconstruction from k-space opens up the possibility of undersampled image acquisition while main-taining high quality in the reconstructed T maps [33]. The proposed method allows for different kinds of undersam-pling and is not limited to Cartesian sampling or single slice acquisitions. While a single receive coil hast been usedfor the current study, the extension to a multi-coil setup is straight forward [33]. The combination of multiple receivecoils and the potential of undersampling k-space could be used to reduce acquisition time in FFC imaging which shallbe subject of a future study. The gained acquisition time might lead to a clinical acceptable scan time using the threeor four fields shown in this work or could be spent to investigate a multitude of different field strength. However,such extensions would require modifications to the phase correction algorithm which is based on images.A potential limitation of the proposed approach is the risk of cross-contamination of the information betweenimages due to the joint regularization [33, 41, 49]. It is assumed that features share the same edge position. If thisassumption is violated in one parameter map, artificial edges might be introduced. The likelihood strongly dependson the used norm for joining the information. As we use a relative weak coupling by means of a Frobenius norm, suchcross-contamination is unlikely. It was shown in previous work that Frobenius norm based joint regularization doesnot show cross-contamination in practice [33, 41, 49]. It might only occur if way too strong regularization weights areused, however, such cases would be discarded in practice as images would look unnatural [41].The proposed reconstruction and fitting approach is integrated into a recently published Python framework for quantitative MRI [46]. This framework allows for an easy adaption to different signal models and thus, a broad appli-cation of the proposed method. Adaptions to the signal model can be made by simply editing text files. In addition,3D regularization strategies are possible which were shown to further improve reconstruction quality [33, 49]. | CONCLUSION
We have successfully introduced joint TGV regularization to multi-field T quantification from FFC imaging. Thehighly significant improvements in T estimation makes it now possible to obtain clinically usable multi-field T maps,and to produce reliable and comparable results. This shows exciting potential for the exploration of low magneticfields and T dispersion effects as illustrated here on two stroke patients. Acknowledgment
This article is based upon work from COST Action CA15209, supported by COST (European Cooperation in Scienceand Technology). Oliver Maier is a Recipient of a DOC Fellowship (24966) of the Austrian Academy of Sciences atthe Institute of Medical Engineering at TU Graz. The authors would like to acknowledge the NVIDIA CorporationHardware grant support. references [1] Steele RM, Korb JP, Ferrante G, Bubici S. New applications and perspectives of fast field cycling NMR relaxometry.Magnetic Resonance in Chemistry 2016;54(6):502–509.[2] Korb JP. Multiscale nuclear magnetic relaxation dispersion of complex liquids in bulk and confinement. Progress inNuclear Magnetic Resonance Spectroscopy 2018 feb;104:12–55. .[3] Broche LM, Ismail SR, Booth NA, Lurie DJ. Measurement of fibrin concentration by fast field-cycling NMR. MagneticResonance in Medicine 2011 oct;67(5):1453–1457. https://doi.org/ . /mrm. .[4] Broche LM, Ashcroft GP, Lurie DJ. Detection of osteoarthritis in knee and hip joints by fast field-cycling NMR. MagneticResonance in Medicine 2012;68(2):358–362.[5] Ruggiero MR, Baroni S, Pezzana S, Ferrante G, Geninatti Crich S, Aime S. Evidence for the Role of Intracellular WaterLifetime as a Tumour Biomarker Obtained by In Vivo Field-Cycling Relaxometry. Angewandte Chemie InternationalEdition 2018 jun;57(25):7468–7472. https://doi.org/ . /anie. .[6] Di Gregorio E, Ferrauto G, Lanzardo S, Gianolio E, Aime S. Use of FCC-NMRD relaxometry for early detection andcharacterization of ex-vivo murine breast cancer. Scientific Reports 2019;9(1):4624. https://doi.org/ . /s - - - .[7] Lurie DJ, Aime S, Baroni S, Booth NA, Broche LM, Choi CH, et al. Fast Field-Cycling Magnetic Resonance Imaging.Comptes Rendus Physique 2010;11(2):136–148.[8] Bödenler M, de Rochefort L, Ross PJ, Chanet N, Guillot G, Davies GR, et al. Comparison of fast field-cycling magneticresonance imaging methods and future perspectives. Molecular Physics 2019 apr;117(7-8):832–848. https://doi.org/ . / . . .[9] Hoelscher UC, Lother S, Fidler F, Blaimer M, Jakob P. Quantification and localization of contrast agents using deltarelaxation enhanced magnetic resonance at 1.5 T. Magnetic Resonance Materials in Physics, Biology and Medicine2012;25(3):223–231. [10] Harris CT, Handler WB, Araya Y, Martínez-Santiesteban F, Alford JK, Dalrymple B, et al. Development and optimizationof hardware for delta relaxation enhanced MRI. Magnetic Resonance in Medicine 2014 oct;72(4):1182–1190. http://doi.wiley.com/ . /mrm. .[11] Chanet N, Guillot G, Willoquet G, Jourdain L, Dubuisson RM, Reganha G, et al. Design of a fast field-cycling magneticresonance imaging system, characterization and methods for relaxation dispersion measurements around 1.5 T. Reviewof Scientific Instruments 2020 feb;91(2):24102. https://doi.org/ . / . .[12] Bödenler M, Basini M, Casula MF, Umut E, Gösweiner C, Petrovic A, et al. R1 dispersion contrast at high field with fastfield-cycling MRI. Journal of Magnetic Resonance 2018;290:68–75.[13] Alford JK, Rutt BK, Scholl TJ, Handler WB, Chronik BA. Delta relaxation enhanced mr: Improving activation - Speeificityof molecular probes through R 1 dispersion imaging. Magnetic Resonance in Medicine 2009;61(4):796–802.[14] Araya YT, Martínez-Santiesteban F, Handler WB, Harris CT, Chronik BA, Scholl TJ. Nuclear magnetic relaxation dispersionof murine tissue for development of T 1 ( R 1 ) dispersion contrast imaging.NMR in Biomedicine 2017 dec;30(12):e3789. http://doi.wiley.com/ . /nbm. .[15] Bödenler M, Malikidogo KP, Morfin JF, Aigner CS, Tóth É, Bonnet CS, et al. High-Field Detection of Biomarkers withFast Field-Cycling MRI: The Example of Zinc Sensing. Chemistry - A European Journal 2019 jun;25(35):8236–8239. https://doi.org/ . /chem. .[16] Carlson JW, Goldhaber DM, Brito A, Kaufman L. MR relaxometry imaging. Work in progress. Radiology 1992sep;184(3):635–639. https://doi.org/ . /radiology. . . .[17] Lurie DJ, Foster MA, Yeung D, Hutchison JMS. Design, construction and use of a large-sample field-cycled PEDRI imager.Physics in Medicine & Biology 1998;43(7):1877. http://stacks.iop.org/ - / /i= /a= .[18] Ungersma SE, Matter NI, Hardy JW, Venook RD, Macovski A, Conolly SM, et al. Magnetic resonance imaging with T1dispersion contrast. Magnetic Resonance in Medicine 2006;55(6):1362–71. .[19] Pine KJ, Goldie F, Lurie DJ. In vivo field-cycling relaxometry using an insert coil for magnetic field offset. MagneticResonance in Medicine 2014;72(5):1492–1497.[20] Romero JA, Rodriguez GG, Anoardo E. A fast field-cycling MRI relaxometer for physical contrasts design and pre-clinicalstudies in small animals. Journal of Magnetic Resonance 2020;311:106682. https://doi.org/ . /j.jmr. . .[21] Broche LM, Ross PJ, Davies GR, MacLeod MJ, Lurie DJ. A whole-body Fast Field-Cycling scanner for clinical molecularimaging studies. Scientific Reports 2019;9(1):10402. https://doi.org/ . /s - - - .[22] Masiewicz E, Ashcroft GP, Boddie D, Dundas SR, Kruk D, Broche LM. Towards applying NMR relaxometry as a diagnostictool for bone and soft tissue sarcomas: a pilot study. Scientific Reports 2020;10(1).[23] Aime S, Botta M, Esteban-Gómez D, Platas-Iglesias C. Characterisation of magnetic resonance imaging (MRI) contrastagents using NMR relaxometry. Molecular Physics 2018;117:898–909. . / . . .[24] Coupe P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An Optimized Blockwise Nonlocal Means Denoising Filter for3-D Magnetic Resonance Images. IEEE Transactions on Medical Imaging 2008;27(4):425–441.[25] Anand CS, Sahambi JS. Wavelet domain non-linear filtering for MRI denoising. Magnetic Resonance Imaging2010;28(6):842–861. X .[26] Fessler JA. Model-Based Image Reconstruction for MRI. IEEE Signal Processing Magazine 2010;27(4):81–89. [27] Knoll F, Bredies K, Pock T, Stollberger R. Second order total generalized variation (TGV) for MRI. Magnetic Resonancein Medicine 2011;65(2):480–491. http:https://doi.org/ . /mrm. .[28] Lingala SG, Hu Y, DiBella E, Jacob M. Accelerated Dynamic MRI Exploiting Sparsity and Low-Rank Structure: k-t SLR.IEEE Transactions on Medical Imaging 2011;30(5):1042–1054.[29] Schloegl M, Holler M, Schwarzl A, Bredies K, Stollberger R. Infimal convolution of total generalized variation functionalsfor dynamic MRI. Magnetic Resonance in Medicine 2017;78(1):142–155. https://onlinelibrary.wiley.com/doi/abs/ . /mrm. .[30] Doneva M, Börnert P, Eggers H, Stehning C, Sénégas J, Mertins A. Compressed sensing reconstruction for magneticresonance parameter mapping. Magnetic Resonance in Medicine 2010;64(4):1114–1120. https://onlinelibrary.wiley.com/doi/abs/ . /mrm. .[31] Sumpf TJ, Uecker M, Boretius S, Frahm J. Model-based nonlinear inverse reconstruction for T2 mapping using highlyundersampled spin-echo MRI. Journal of Magnetic Resonance Imaging 2011;34(2):420–428. https://onlinelibrary.wiley.com/doi/abs/ . /jmri. .[32] Wang X, Roeloffs V, Klosowski J, Tan Z, Voit D, Uecker M, et al. Model-based T1 mapping with sparsity constraintsusing single-shot inversion-recovery radial FLASH. Magnetic Resonance in Medicine 2018;79(2):730–740. https://onlinelibrary.wiley.com/doi/abs/ . /mrm. .[33] Maier O, Schoormans J, Schloegl M, Strijkers GJ, Lesch A, Benkert T, et al. Rapid T1 quantification from high resolution3D data with model-based reconstruction. Magnetic Resonance in Medicine 2019 mar;81(3):2072–2089. https://doi.org/ . /mrm. .[34] Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D: Nonlinear Phenomena1992;60(1):259–268. F .[35] Bredies K, Kunisch K, Pock T. Total Generalized Variation. SIAM Journal on Imaging Sciences 2010 jan;3(3):492–526. https://doi.org/ . / .[36] Lustig M, Donoho D, Pauly JM. Sparse MRI: The application of compressed sensing for rapid MR imaging. MagneticResonance in Medicine 2007 dec;58(6):1182–1195. https://doi.org/ . /mrm. .[37] Wang X, Roeloffs V, Klosowski J, Tan Z, Voit D, Uecker M, et al. Model-based T1 mapping with sparsity constraintsusing single-shot inversion-recovery radial FLASH. Magnetic Resonance in Medicine 2018 feb;79(2):730–740. https://doi.org/ . /mrm. .[38] Lai Z, Zhang X, Guo D, Du X, Yang Y, Guo G, et al. Joint sparse reconstruction of multi-contrast MRI images with graphbased redundant wavelet transform. BMC Medical Imaging 2018;18(1):7. https://doi.org/ . /s - - -y .[39] Hógáin DÓ, Davies GR, Baroni S, Aime S, Lurie DJ. The use of contrast agents with fast field-cycling magnetic resonanceimaging. Physics in Medicine and Biology 2010 nov;56(1):105–115. https://doi.org/ . / - / / / .[40] Bredies K. Recovering Piecewise Smooth Multichannel Images by Minimization of Convex Functionals with Total Gen-eralized Variation Penalty BT - Efficient Algorithms for Global Optimization Methods in Computer Vision. Berlin, Hei-delberg: Springer Berlin Heidelberg; 2014. p. 44–77.[41] Knoll F, Holler M, Koesters T, Otazo R, Bredies K, Sodickson DK. Joint MR-PET Reconstruction Using a Multi-ChannelImage Regularizer. IEEE Transactions on Medical Imaging 2017;36(1):1–16.[42] Chambolle A, Pock T. A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging. Journalof Mathematical Imaging and Vision 2011;40(1):120–145. https://doi.org/ . /s - - - . [43] Malitsky Y, Pock T. A First-Order Primal-Dual Algorithm with Linesearch. SIAM Journal on Optimization 2018jan;28(1):411–432. https://doi.org/ . / M .[44] Ross PJ, Broche LM, Lurie DJ. Rapid field-cycling MRI using fast spin-echo. Magnetic Resonance in Medicine2015;73(3):1120–1124. https://onlinelibrary.wiley.com/doi/abs/ . /mrm. .[45] Broche LM, Ross PJ, Davies GR, Lurie DJ. Simple algorithm for the correction of MRI image artefacts due to randomphase fluctuations. Magnetic Resonance Imaging 2017;44:55–59.[46] Maier O, Spann SM, Bödenler M, Stollberger R. PyQMRI: An accelerated Python based Quantitative MRI toolbox. Journalof Open Source Software 2020;5(56):2727. https://doi.org/ . /joss. .[47] Jin Q, Zhong M. On the iteratively regularized Gauss–Newton method in Banach spaces with applications to parameteridentification problems. Numerische Mathematik 2013;124(4):647–683.[48] Petrov OV, Stapf S. Parameterization of NMR relaxation curves in terms of logarithmic moments of the relaxation timedistribution. Journal of Magnetic Resonance 2017;279:29–38.[49] Huber R, Haberfehlner G, Holler M, Kothleitner G, Bredies K. Total generalized variation regularization for multi-modalelectron tomography. Nanoscale 2019;11(12):5617–5632. https://doi.org/ . /c nr k . A | MATHEMATICAL DERIVATIONS
The main computational burden lies in the optimization of the convex inner problems of each GN iteration. Eachconsists of an iterative solution to the following optimization taskmin u , v (cid:107) DS | u = u k u − ˜ d k (cid:107) + γ k ( β (cid:107) (cid:43) u − v (cid:107) , , F + β | (cid:107) E v (cid:107) , , F ) + δ k (cid:107) u − u k (cid:107) M k . (9)The (cid:107) · (cid:107) , , F terms resemble the Frobenius type TGV functionals, joining common spatial information of theunknown parameter maps and are defined as (cid:107) v (cid:107) , , F = (cid:213) x , y (cid:118)(cid:117)(cid:116) N u (cid:213) l =1 | v , lx , y | + | v , lx , y | (10)with v = ( v , l , v , l ) N u l =1 ∈ U × N u constituting the approximation of 2D spatial derivatives using (cid:43) , and for the sym-metrized gradient E χ = ( χ , l , χ , l , χ , l ) N u l =1 ∈ U × N u (cid:107) χ (cid:107) , , F = (cid:213) x , y (cid:118)(cid:117)(cid:116) N u (cid:213) l =1 | χ , lx , y | + | χ , lx , y | + 2 | χ , lx , y | . (11) (cid:43) : U N u → U × N u and E : U × N u → U × N u are defined as (cid:43) u = (cid:16) δ x + u l , δ y + u l (cid:17) N u l =1 and E v = (cid:16) δ x − v , l , δ y − v , l , δ y − v , l + δ x − v , l (cid:17) N u l =1 . The gradient (cid:43) and symmetrized gradient E operations are computed by taking the finite differences, defined as for-ward δ x + , δ y + and backward δ x − , δ y − differences with respect to the spatial directions along ( x , y ) . At the boundaries,the image is symmetrically extended. To achieve the saddle point formulationmin x max y (cid:104) K x , y (cid:105) + G ( x ) − F ∗ ( y ) , (12)as required for the applied primal-dual algorithm, the linearized problem of Eq. 6 can be reformulated by applying theconvex conjugate: min x = ( u , v ) (cid:107) DS | u = u k u − ˜ d k (cid:107) + γ k ( β (cid:107) (cid:43) u − v (cid:107) , , F + β | (cid:107) E v (cid:107) , , F ) + δ k (cid:107) u − u k (cid:107) M k ⇔ min x max y = ( z , z , r ) (cid:26)(cid:10) DS | u = u k | u = u k u , r (cid:11) − (cid:68) ˜ d k , r (cid:69) − (cid:107) r (cid:107) (cid:27) + (cid:104) K x , z (cid:105) − I {(cid:107)·(cid:107) ∞≤ β γk } ( z ) − I {(cid:107)·(cid:107) ∞≤ β γk } ( z ) + δ k (cid:107) u − u k (cid:107) M k ⇔ min x max y (cid:104) K x , y (cid:105) + G ( x ) − F ∗ ( y ) . with K = (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) DS − i d E (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) , K = (cid:32) (cid:43) − i d E (cid:33) , z = ( z , z ) T . F ∗ ( y ) = (cid:26)(cid:10) DS | u = u k u , r (cid:11) − (cid:68) ˜ d k , r (cid:69) + 12 (cid:107) r (cid:107) (cid:27) + I {(cid:107)·(cid:107) ∞≤ β γk } ( z ) + I {(cid:107)·(cid:107) ∞≤ β γk } ( z ) , G ( x ) = δ k (cid:107) u − u k (cid:107) M k . The projection I {(cid:107)·(cid:107) ∞≤ αpγk } ( z p ) stems from to the convex conjugate of the L -norm which amounts to the indicatorfunction of the L ∞ -norm unit ball, scaled by the corresponding regularization parameter α p γ k I {(cid:107)·(cid:107) ∞≤ αpγk } ( z p ) = (cid:107) z p (cid:107) ∞ ≤ α p γ k ∞ el se DS | u = u k is the Jacobian matrix of S evaluated at u = u k of the non-linear FFC signal equation for all fields B and theircorresponding inversions times t E i : DS : u = ( u l ) N u l =1 ↦→ (cid:169)(cid:173)(cid:171) N u (cid:213) l =1 ∂S B Ei , t Ein ∂u l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u = u k u l (cid:170)(cid:174)(cid:172) N d n =1 = ( ξ n ) N d n =1 . (13)To update steps of the PD algorithm are defined as y n +1 = ( i d + σ∂F ∗ ) − ( y n + σK x n ) x n +1 = ( i d + τ∂G ) − ( x n − τK H y n +1 ) x n +1 = x n +1 + θ ( x n +1 − x n ) , (14)with i d amounting to the identity matrix and θ ∈ [ , ] . To compute the updates, additional operations need to bedefined, which will be covered in the next few paragraphs.First, the adjoint operations to the linear operator K of the forward problem, termed K H with H being the Hermi-tian transpose operation, is defined as K H = (cid:32) DS H − div − i d − div (cid:33) , (15)where the divergence operators div and div are the negative adjoints of (cid:43) and E , respectively. The adjoint of theJacobi matrix DS H amounts to a simple complex transpose operation in matrix notation, which can be written inoperator notation as DS H : ξ = ( ξ n ) N d n =1 ↦→ (cid:169)(cid:173)(cid:173)(cid:171) N d (cid:213) n =1 ∂S B Ei , t Ein ( u ) ∂u l (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) u = u k ξ n (cid:170)(cid:174)(cid:174)(cid:172) N u l =1 = ( u l ) N u l =1 = u . Finally, the proximal maps, projecting on the function F ∗ amount to simple point wise operations P β ( ξ ) = ξ max (cid:16) , | ξ | β γ (cid:17) , P β ( ξ ) = ξ max (cid:16) , | ξ | β γ (cid:17) , P σL ( ξ ) = ξ − σ ˜ d k σ and the proximal map of G is defined as P G ( ξ ) = ( i d + τδ k M k ) − ( τδ k M k u k + ξ ) . Making use of the fact that M k is a diagonal matrix, the inversion of ( i d + τδ k M k ) is trivial and thus P G ( ξ ) can becomputed in a point-wise fashion. B | PSEUDO CODE Initialize: ( u , v ) , ( u , v ) , ( z , z , r ) , τ > , κ = 1 , θ = 1 , µ = 0 . Iterate: Primal Update: u m +1 ← P τ m G (cid:16) u m − τ m (cid:16) − div z m + DS H r m (cid:17)(cid:17) v m +1 ← v − τ m (cid:16) − div z m − z m (cid:17) Update κ and τ : κ m +1 ← κ m ( δ k τ m ) τ m +1 ← τ m (cid:113) κ m κ m +1 ( θ m ) Start Linesearch: Update θ : θ m +1 ← τ m +1 τ m Extrapolation: ( u m +1 , v m +1 ) ← ( u m +1 , v m +1 ) + θ m +1 ( ( u m +1 , v m +1 ) − ( u m , v m )) Dual Update: z m +10 ← P β (cid:16) z m + κ m +1 τ m +1 ( (cid:43) u m +1 − v m +1 ) (cid:17) z m +11 ← P β (cid:16) z m + κ m +1 τ m +1 ( E v m +1 ) (cid:17) r m +1 ← P κ m +1 τ m +1 L (cid:16) r m + κτ m +1 ( DS u m +1 ) (cid:17) break Linesearch if: √ κ m +1 τ m +1 (cid:107) K H y m +1 − K H y m (cid:107) ≤ (cid:107) y m +1 − y m (cid:107) else: τ m +1 ← τ m +1 µ Update: ( u m , v m , τ m ) ← ( u m +1 , v m +1 , τ m +1 ) ( z m , z m , r m ) ← ( z m +10 , z m +11 , r m +1 ) Algorithm 1:
Pseudo code for the solution of the inner problems of each GN step based on the primal-dual algo-rithm with linesearch. As described in the original publication, linearity of certain operations can be used to speedup the computation of the linesearch procedure. List of Figures B P (dashed line) or by settingthe polarization field to the detection field B D , i.e., B P = B D (solid line). After an inversion pulse at t the longitudinal magnetization M z evolves at the desired evolution field B E for a given evolution time.The following magnetization M z ( t evo ) can be detected by any MRI acquisition module. Note that theMRI signal is both inverted and detected at B D . M and M E represent the equilibrium magnetizationfor B D and B E , respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212 FFC images of a stroke patient from the PUFFINS study (patient I). Image obtained at 200 mT evolutionfield show good signal after inversion (b) and after 455 ms evolution time (a) but low contrast in thelesion, while low-field images at 21 mT (c) show good lesion contrast but most other tissue shows littleto no signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 Multi-field T maps obtained from simulated FFC imaging inversion recovery data. The reference T maps for three different evolution fields (200 mT, 21.1 mT and 2.2 mT) are shown at the top. In the leftcolumn, the multi-field T maps were obtained by pixel-wise fitting of each evolution field separately(standard), and the right column results from joint model-based reconstruction of all three evolutionfields together (proposed). The noise level increases from top to bottom from 1 % to 4 %. . . . . . . . . 224 Pixel-wise relative absolute difference to the ground truth T values. Numbers next to the differenceimages show mean relative absolute error within the phantom. All values are given in percent. . . . . 235 2D histogram evaluation for T maps in Figure 3, which were obtained from synthetic FFC imagingdata by pixel-wise fitting (standard) and joint model-based reconstruction (proposed). The dashed linerepresents identity. Shown are reference values on the ordinate versus results obtained with NLLS orTGV on the abscissa. Points below the identity line correspond to under-estimation, points above toover-estimation, respectively. All values are given in ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . 246 In vivo multi-field T maps of a transverse slice of the brain of stroke patient I. T maps were obtainedat three different evolution fields B E = { , . , . } mT by pixel-wise fitting of the signal model foreach B E separately (top row) and by the proposed multi-field model-based reconstruction approachutilizing the joint information of all three evolution fields (bottom row). . . . . . . . . . . . . . . . . . . . 257 In vivo multi-field T maps of a transverse slice of the brain of stroke patient II. T maps were obtainedat four different evolution fields B E = { , , . , . } mT by pixel-wise fitting of the signal modelfor each B E separately (top row) and by the proposed multi-field model-based reconstruction approachutilizing the joint information of all four evolution fields (bottom row). . . . . . . . . . . . . . . . . . . . 268 T dispersion profiles obtained from the patients I (solid lines) and II (dashed lines) in the regions ofsubcutaneous fat (SC fat, orange) measured between the scalp and the brain, grey matter (GM, grey)measured over two centimeters of the cortical region, white matter (WM, blue) measured over theinner region of the lobes and the lesion (red). The error bars stand for twice the standard deviation ofthe T values measured across the ROIs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 List of Tables List of Supporting Figures
S1 Absolute value of C map obtained from simulated FFC imaging inversion recovery data. The reference C map is shown at the top. In the left column, the multi-field C maps were obtained by pixel-wisefitting of each evolution field separately (standard), and the right column results from joint model-basedreconstruction of all three evolution fields together (proposed). The noise level increases from top tobottom from 1 % to 4 %. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27S2 Phase value of C map obtained from simulated FFC imaging inversion recovery data. The reference C map is shown at the top. In the left column, the multi-field C maps were obtained by pixel-wisefitting of each evolution field separately (standard), and the right column results from joint model-basedreconstruction of all three evolution fields together (proposed). The noise level increases from top tobottom from 1 % to 4 %. Simulated phase on C amounts to 0 rad. . . . . . . . . . . . . . . . . . . . . . . 28S3 Absolute value of multi-field α maps obtained from simulated FFC imaging inversion recovery data.The reference α maps for three different evolution fields (200 mT, 21.1 mT and 2.2 mT) are shown atthe top. In the left column, the multi-field α maps were obtained by pixel-wise fitting of each evolutionfield separately (standard), and the right column results from joint model-based reconstruction of allthree evolution fields together (proposed). The noise level increases from top to bottom from 1 % to4 %. Values next to each figure represent the mean value within the simulated phantom in a.u.. . . . . 29S4 Phase of multi-field α maps obtained from simulated FFC imaging inversion recovery data. The refer-ence phase maps for three different evolution fields (200 mT, 21.1 mT and 2.2 mT) are shown at the top.In the left column, the multi-field α phase maps were obtained by pixel-wise fitting of each evolutionfield separately (standard), and the right column results from joint model-based reconstruction of allthree evolution fields together (proposed). The noise level increases from top to bottom from 1 % to4 %. Values next to each figure represent the mean value within the simulated phantom in radiant. . . 30S5 Regions of interest selected to extract the dispersion profiles in Figure 8 in patients I (left) and II (right).The regions for white matter are delineated in light blue dashed lines, grey matter in solid dark bluelines, fat in yellow dotted lines and lesions in red dot-dashed lines. . . . . . . . . . . . . . . . . . . . . . 31 TA B L E 1
Evolution times and fields used to generate the simulated images and corresponding timings for the invivo acquisitions. Application Variable Field (mT) Evolution time (ms) t E n
200 455 242 129 68 36 t E n t E n t E n
200 455 242 129 68 36 t E n t E n t E n
200 455 196 84 36 t E n
37 338 145 63 27 t E n t E n TA B L E 2
Parameter values selected to generate the simulated images.ROI 1 2 3 4Parameter a 5.6 4.4 2.6 3.8power law b -0.1 -0.15 -0.3 -0.08at 200 mT 152 178.5 237.3 231.4at 21 mT 121.3 127.3 120.7 193.2 T (ms) at 2.2 mT 96.8 90.8 61.3 161.3at 200 mT 1/0.5236at 21 mT 0.75/0.6981 α abs (a.u.)/phase (rad) at 2.2 mT 0.6/0.8727 C (a.u.) at all field 1 1/3 2/3 2.03/3 F I G U R E 1
Pulse sequence diagram of an inversion recovery FFC imaging pulse sequence. A pre-polarization ofthe sample magnetization can be applied by cycling the main magnetic field to B P (dashed line) or by setting thepolarization field to the detection field B D , i.e., B P = B D (solid line). After an inversion pulse at t the longitudinalmagnetization M z evolves at the desired evolution field B E for a given evolution time. The following magnetization M z ( t evo ) can be detected by any MRI acquisition module. Note that the MRI signal is both inverted and detected at B D . M and M E represent the equilibrium magnetization for B D and B E , respectively. F I G U R E 2
FFC images of a stroke patient from the PUFFINS study (patient I). Image obtained at 200 mTevolution field show good signal after inversion (b) and after 455 ms evolution time (a) but low contrast in the lesion,while low-field images at 21 mT (c) show good lesion contrast but most other tissue shows little to no signal. F I G U R E 3
Multi-field T maps obtained from simulated FFC imaging inversion recovery data. The reference T maps for three different evolution fields (200 mT, 21.1 mT and 2.2 mT) are shown at the top. In the left column, themulti-field T maps were obtained by pixel-wise fitting of each evolution field separately (standard), and the rightcolumn results from joint model-based reconstruction of all three evolution fields together (proposed). The noiselevel increases from top to bottom from 1 % to 4 %. F I G U R E 4
Pixel-wise relative absolute difference to the ground truth T values. Numbers next to the differenceimages show mean relative absolute error within the phantom. All values are given in percent. F I G U R E 5
2D histogram evaluation for T maps in Figure 3, which were obtained from synthetic FFC imagingdata by pixel-wise fitting (standard) and joint model-based reconstruction (proposed). The dashed line representsidentity. Shown are reference values on the ordinate versus results obtained with NLLS or TGV on the abscissa.Points below the identity line correspond to under-estimation, points above to over-estimation, respectively. Allvalues are given in ms. F I G U R E 6
In vivo multi-field T maps of a transverse slice of the brain of stroke patient I. T maps were obtainedat three different evolution fields B E = { , . , . } mT by pixel-wise fitting of the signal model for each B E separately (top row) and by the proposed multi-field model-based reconstruction approach utilizing the jointinformation of all three evolution fields (bottom row). F I G U R E 7
In vivo multi-field T maps of a transverse slice of the brain of stroke patient II. T maps wereobtained at four different evolution fields B E = { , , . , . } mT by pixel-wise fitting of the signal model foreach B E separately (top row) and by the proposed multi-field model-based reconstruction approach utilizing thejoint information of all four evolution fields (bottom row). F I G U R E 8 T dispersion profiles obtained from the patients I (solid lines) and II (dashed lines) in the regions ofsubcutaneous fat (SC fat, orange) measured between the scalp and the brain, grey matter (GM, grey) measured overtwo centimeters of the cortical region, white matter (WM, blue) measured over the inner region of the lobes and thelesion (red). The error bars stand for twice the standard deviation of the T values measured across the ROIs. C | SUPPORTING INFORMATION
S U P P O R T I N G I N F O R M AT I O N F I G U R E S 1
Absolute value of C map obtained from simulated FFCimaging inversion recovery data. The reference C map is shown at the top. In the left column, the multi-field C mapswere obtained by pixel-wise fitting of each evolution field separately (standard), and the right column results fromjoint model-based reconstruction of all three evolution fields together (proposed). The noise level increases from topto bottom from 1 % to 4 %. S U P P O R T I N G I N F O R M AT I O N F I G U R E S 2
Phase value of C map obtained from simulated FFC imaginginversion recovery data. The reference C map is shown at the top. In the left column, the multi-field C maps wereobtained by pixel-wise fitting of each evolution field separately (standard), and the right column results from jointmodel-based reconstruction of all three evolution fields together (proposed). The noise level increases from top tobottom from 1 % to 4 %. Simulated phase on C amounts to 0 rad. S U P P O R T I N G I N F O R M AT I O N F I G U R E S 3
Absolute value of multi-field α maps obtained fromsimulated FFC imaging inversion recovery data. The reference α maps for three different evolution fields (200 mT,21.1 mT and 2.2 mT) are shown at the top. In the left column, the multi-field α maps were obtained by pixel-wisefitting of each evolution field separately (standard), and the right column results from joint model-basedreconstruction of all three evolution fields together (proposed). The noise level increases from top to bottom from1 % to 4 %. Values next to each figure represent the mean value within the simulated phantom in a.u.. S U P P O R T I N G I N F O R M AT I O N F I G U R E S 4