Optimal Experiment Design for Magnetic Resonance Fingerprinting: Cramer-Rao Bound Meets Spin Dynamics

IEEE TRANSACTIONS ON MEDICAL IMA GING 1 Optimal Experiment Design for Magnetic Resonance Fingerprinting: Cram ´ er -Ra o Bound Meets Spin Dynamics Bo Zhao, Member , IEEE , J ustin P . Ha ldar , Se nior Member , IEEE , Con g yu Liao, Dan Ma, Y un Jiang, Mark A . Griswold, Kawin Setsompop, a nd Lawrenc e L. W ald, Memb er , IEEE Abstract —Magnetic resonance (MR) ﬁngerprinting is a n ew quantitative imaging paradigm, which si multaneously acqu i res multiple MR tissue parameter maps in a single experiment. In this p aper , we present an estimation-theoretic framework to perfo rm experiment design for M R ﬁn gerprinting. Speciﬁ cally , we describe a discrete-time dynamic system to model spin dynamics, and deri ve an estimation-theoretic bound , i.e., the Cram ´ er -Rao bound (CRB), to characterize the signal-t o-noise ratio (SNR) efﬁciency of an MR ﬁngerprinting experiment. W e then formulate an optimal experiment design problem, which determines a sequence of acquisition p arameters to encode M R tissue param eters with the maximal S NR efﬁciency , while respect- ing the physical constraints and other constraints from the image decoding/reconstruction process. W e ev aluate the performance of the p roposed approa ch with numerical simulations, p hantom experiments, and in viv o experiments. W e demonstrate that the opti mized experiments substantially reduce data acquisition time and/or improv e p arameter estimation. For example, the optimized experiments achiev e about a factor of two impro vement in t h e accuracy of T 2 maps, whi le keeping similar or slightly better accuracy of T 1 maps. Finally , as a remarkable observ ation, we ﬁn d that the sequence of opt i mized acqu i sition parameters appears to be highly structured rather th an randomly/pseudo- randomly varying as is prescribed in the con ventional MR ﬁngerprinting experiments. Copyri ght (c) 2017 IEEE. Personal use of this material is permitt ed. Ho wev er , permission to use this material for any other purposes must be obtaine d from the IEEE by sending a request to pubs-permissions@ieee .org This work was supported, in part, by the National Institutes of Health under Grant No. NIH-R01-EB017219, NIH-R01-EB017337, and NIH-P41- EB015896, and the National Science Foun dation under Grant No. NSF-CCF- 135063. B. Z hao was also supported by a Ruth L. Kirschstein Nation al Researc h Service A ward Postdoctoral Fello wship from the National Institutes of Health under Grant No. NIH-F32-EB024381. B. Zhao is with the Athinoula A. Martino s Center for Biomedical Imaging, Massachuset ts General Hospital, Charlest own, MA 02129 USA, and also with the Department of Radiology , Harvard Medical School, Boston, MA 02115 USA (email: bzhao4@mgh.ha rvard .edu). J. P . Haldar is w ith the Signal and Image Processing Institute and Ming Hsieh Department of Electri cal Engineeri ng, Univ ersity of Southern Califor- nia, Los Angeles, CA 90089 USA (jhaldar@usc.edu ). C. Liao is with the Athinoula A. Martinos Center for Biomedica l Imaging, Massachuset ts General Hospital, Charlest own, MA 02129 USA, and also with the Department of Biomedica l E ngineeri ng, Zhejiang Univ ersity , Hangzhou, Zhejiang Provin ce 310027 China (email: cliao2@mgh.harv ard.edu). D. Ma, Y . Jiang, and M. A. Griswold are with the Department of Radi- ology , Case W estern Reserve Univ ersity , Cle vel and, OH 44106 USA (email: dan.ma@case .edu, yun.jiang@case.ed u, m ark.griswold@case .edu). K. Setsompop and L . L. W ald are with the Athinoul a A. Martinos Center for Biomedic al Imaging, Massachusetts General Hospital, Charlest own, MA 02129 USA, and also with the Department of Radiol ogy , Harv ard Medical School, Boston, MA, 02115 USA, and also with the Harvard-MIT Divi - sion of Health Sciences and T echnology , Massachusett s Institute of T ech- nology , Cambridge, MA 02139 USA (email: kawin@n mr .mgh.harva rd.edu, wald @nmr .mgh.harva rd.edu). Index T erms —Optimal experiment design, Cram ´ er -Rao bound, spin dynamics, system identiﬁcation, statistical inference, q u an- titative imaging. I . I N T RO D U C T I O N M A GNETIC resonanc e (MR) ﬁnger p rinting is a novel quantitative magnetic resonan ce imag ing (MRI) fram e- work [1]. It enables simultaneous acquisition of multiple MR tissue parameter s (e .g., T 1 , T 2 , and spin d ensity) in a single imag ing experime nt. Compare d to co n ventional MR relaxometr y technique s (e.g ., [2] –[5]), th e or iginal pr oof-o f- principle MR ﬁngerp rinting implementatio n featur ed several key innovations in its enco ding and decod ing pr ocesses [ 1]. Speciﬁcally , d u ring enc o ding, it applies a series of tim e- varying, ran d om or q uasi-rando m data acqu isition para m eters (e.g., ﬂip angles and rep etition times) to probe the sp in system. This g enerates uniq ue transient-state signal evolutions, or ﬁng erprints, f o r different MR tissue parameter s. I t f u rther applies incoh erent spatial encod ing (e.g., variable density spiral ac q uisition) to collect k-spac e data. During decodin g, it perfo rms a simple g r idding reconstructio n to reco nstruct contrast-weig h ted imag es, after wh ich it a pplies a diction ary- based pattern ma tch ing to obtain MR tissue pa r ameter maps of intere st. Giv en its m any de p artures from conventional qu antitative imaging, so m e fundamental questions still remain unclear about the me chanism of MR ﬁngerprinting. For example, from a theoretical perspec tive, the op tim ality of the encodin g and de c o ding pro cesses for MR ﬁnger printing has not b een examined in [1]. While a rece n t pa p er has perfo r med a th eo- retical analysis of the MR ﬁnger p rinting p roblem from a low- dimensiona l man ifold r e covery v iewpoint [6], this analysis is asymptotic and p robab ilistic in n ature, and does n ot pr ovide results that co uld be u sed to ensur e the quality o r gu id e the design of ﬁnite-dur ation MR ﬁngerpr inting expe riments. From a practical p erspective, while existing M R ﬁngerprin t- ing metho ds work well in so me scenarios, th e r e ar e o ther importan t cases where the perf ormanc e of the se metho ds can be much worse. For example, the accu racy of T 2 maps fr om MR ﬁngerp rinting often dep ends critically o n the length of data acq uisition, and is mu ch worse than that of T 1 maps, especially with short acquisition lengths [7]– [ 11]. In these set- tings, it would b e highly desirable to either improve paramete r mapping quality with out in creasing experiment du ration, or IEEE TRANSACTIONS ON ME DICAL IMA GING 2 reduce exper iment du r ation without com p romising parameter mapping q uality . Recently , we introdu ced a novel statistical imag ing fram e - work for MR ﬁngerprin tin g [7], [8 ]. On th e decodin g side, we propo sed a maximu m likelihood (ML) a pproac h that directly reconstru cts MR tissue p arameter maps from highly- undersam pled, noisy k -space data. Mor eover , we showed analytically th at the co n ventiona l reconstru ction ap proach [1 ] is sub -optimal fro m a statistical estimatio n perspective [7], [8]. The use of ML reco nstruction has dramatically improved the estimation accuracy and /or redu ced acqu isition time. Apart from the qua lity of th e decodin g scheme s, th e perfor mance of MR ﬁng erprintin g inheren tly d epends on the q uality of the available d ata. Th is motivates the optimizatio n of MR ﬁngerpr in ting on the encod ing side, which natur ally falls into the dom ain of optim a l experim ent de sig n [12], [13] in the statistical ima g ing framework. In this paper, we address the o ptimal experim ent de sig n problem for MR ﬁng erprin tin g. Our goal is to en c o de M R tissue p arameters into the m ost in formative measure m ents in the presence o f n o ise. Speciﬁcally , we ﬁrst mo d el spin dyn am- ics with a d iscrete-time dy namic system. W e then calculate an estimation-theo retic boun d, i. e ., the Cra m ´ er-Rao bou nd (CRB) [14], as a m easure o f th e sign al-to-no ise (SNR) efﬁciency for MR ﬁngerprin ting experiments. W e fur ther utilize this bo und to for mulate a n optim al experiment design prob lem to cho ose MR ﬁngerp rinting acquisition para m eters for the ma ximal SNR efﬁciency , while respecting both phy sical constrain ts an d other con straints fro m image deco ding. W e show represen ta- ti ve results fro m numerical simulations, pha n tom exp eriments, and in viv o experiments to illustrate the per forman ce of the propo sed f r amework. As a remar kable ob servation, we ﬁnd that the sequence of optimized acquisition parame ters app ears to be highly struc tu red rath er than ra n domly /p seudo-r andomly varying as used in th e existing MR ﬁngerpr in ting experiments (e.g., [1], [15]–[1 8 ]). A pr eliminary account of this work was presented in o ur early co nferenc e p apers [1 9]–[21], which, to the b est o f o ur knowledge, ﬁrst introduced th e u se o f estimation-theore tic bound s f or MR ﬁn gerprin ting exper iment design. Rece n tly , similar pro blems have also been investigated by o ther re- searchers. For example, Assl ¨ ander et al. applied our CRB- based experimen t desig n framework to a new MRF imag- ing sequence [ 16] u nder th e po lar coo rdinates of the Bloch equation [22]. Maid e ns et al. provided an optimal co ntrol interpretatio n of the expe riment design p roblem, for wh ich they developed a dyn amic progra mming b a sed algo rithm [23]. Howe ver , du e to the curse of dimen sionality f or d y namic progr amming [24], th e feasibility of their app roach was only demonstra te d in a high ly simpliﬁed scenario that d oes not reﬂect the full complexities of rea l MRF ap plications. Besides CRB-based ap p roache s, the experimen t d esign p roblem was also ad dressed f rom other perspectives. For examp le, Cohen et al. optimize d MR ﬁn gerprin ting acqu isition parame te r s by maximizing the d iscrimination p ower between d ifferent tissue types [25]. Note that the CRB ha s been previously used in an alyzing and d esigning exper iments for conv entio n al MR relax o metry (e.g., [26]–[34]). Here we extend these ap p roache s to MR ﬁngerpr in ting, which is a uniq ue tran sient-state q u antitative imaging technique that en codes MR tissue p arameters into spin dyna m ics. In this con text, the CRB calculation can b e more comp lex, wh ich requires hand ling a dy namic system, rather than an analytical signal model as in c o n ventiona l MR relaxo metry [2 6]–[34]. In this paper, we introd uce a n efﬁcient way of ca lc u lating the CRB for dyna m ic systems described b y the Bloch eq uation. In particular, we show that the CRB calculatio n can be done b y iteratin g a set of difference equations. For easy r eference, we summar iz e he r e th e key n otations and symbols used in the paper . W e use R to denote the ﬁeld of real nu mbers, and use R n and R m × n to respectively den ote the space o f real vectors of length n and the space o f real m × n matrices. W e u se bold letter s (e.g., x or X ) to d enote vecto rs o r matrices. W e respectively u se X T , X − 1 , and tr ( X ) to deno te the transpo se, in verse, and trace o f X . For a scalar-v alued function f : R N → R with the argumen t x ∈ R N , we deﬁne its g r adient, i.e., ∂ f /∂ x , as an N × 1 vector with [ ∂ f /∂ x ] n = ∂ f /∂ x n ; f o r a vector-valued function F : R N → R M with the argument x ∈ R N , we deﬁne its Jacob ian matr ix, i.e., ∂ F /∂ x , as an M × N m atrix with [ ∂ F /∂ x ] m,n = ∂ F m /∂ x n ; and for a matr ix-valued fun ction F : R → R M × N with the argumen t x ∈ R , we deﬁne its deriv ativ e as an M × N matrix with [ ∂ F /∂ x ] m,n = ∂ F m,n /∂ x . For a r andom vector x ∈ R N , we deno te its expectation as E [ x ] ∈ R N , and its covariance matrix as Cov ( x ) = E h ( x − E ( x ) ) ( x − E ( x )) T i ∈ R N × N . The re st of th e pap er is o rganized as follows. Section II p r esents the prop osed framework in detail, wh ic h starts with a state-space mo del for spin d ynamics, followed by the calculation of the Cram ´ er-Rao boun d a n d the optim al design of MR ﬁnge r printing experiments. Section I II demon strates the perfor mance of the pr o posed ap p roach with n umerical simulations, p hantom experim ents, a nd in viv o experiments. Section IV discusses th e related issues and future work, followed by the con cluding remark s in Section V. I I . P RO P O S E D F R A M E W O R K A. S ignal Model A numbe r of MR ﬁnge rprinting imaging sequences have re- cently b een developed (e.g., [1], [ 15]–[18]). I n this subsection, we star t by form alizing a gene r ic state-space mo del [3 5] f or spin dynam ics that unde rlie all MR ﬁngerprin ting sequences, and then giv e an example b y specializin g th is m odel to a representative MR ﬁnger printing sequence . 1) State- Space Mod el: Spin dynamics are governed by the Bloch equ ation [36], which is a system of ﬁrst-or der ordin ary differential eq uations. While the mag netization ev olves in con - tinuous tim e, w e are mainly interested in its values a t a ﬁn ite set of time p oints in an MR ﬁng erprintin g experiment. H e re we c o nsider a d iscrete-time state-space model for simplicity , which com pletely capture s th e fea tures of interest fo r spin dynamics in co ntinuou s time . M o reover , we focu s on the magnetizatio n ev olution with respec t to a small sample o f tissue (i.e . , a voxel); other issues related to spatial encodin g will be d iscussed later . IEEE TRANSACTIONS ON ME DICAL IMA GING 3 In the presence of magne tic ﬁeld inho m ogeneity (e.g., with the use of g radients), in trav oxel spin deph asing occurs. T o ac- count for this effect, we mo del mu ltiple isochromats in a voxel, each o f which repr esents an ensemble of spin s that have the same re so nance freq u ency [37]. Here we use the isochro mat- summation approac h [ 3 8], [39] to ap prox im ate the mag ne- tization evolution f or the voxel of in terest. This appro ach perfor ms Bloch simulations for each individual isochroma t, and then sums up their magn etization evolutions to obtain the ensemble m agnetization ev olution for a voxel. From a dynamic system viewpoint, th is a p proach essentially emp loys an in depend ent state equ ation to describ e the magn e tiza tion ev olution of each isoch romat, an d then obtains the ensemb le magnetizatio n throu gh an ob servation eq u ation. In the f o llowing, we start by describing th e state eq uation for a single isoch r omat. Here we divide th e voxel of interest, i.e., △ , into a set o f sufﬁciently small, eq ual-sized subvoxels {△ r } , where △ r is centered a t r and △ r ⊂ △ . W e assume that there exists a sing le isochromat f or each △ r . Let M r [ n ] ∈ R 3 denote the magn etization associated with th e isoch r omat at △ r at the end of the n th re p etition time (immed iately before the next signal excitation). The ma g netization ev olution can be describ e d by th e following state equ ation: M r [ n ] = A r ( u [ n ] , θ ) M r [ n − 1] + B r ( u [ n ] , θ ) , (1) for n = 1 , · · · , N , wher e θ ∈ R p contains the unk n own parameters in an MR ﬁngerprinting experimen t, includ ing the tissue-speciﬁc p arameters (e.g., T 1 , T 2 , and sp in den sity) and experiment-spe ciﬁc p arameters (e.g., off-resonan ce fr e - quency); u [ n ] ∈ R q contains the data acqu isition pa r ameters applied d u ring the n th repetition time , includ ing the ﬂip an g le α n , the phase of the radio freque ncy (RF) pu lse φ n , th e ech o time T E n , and the repetition time T R n ; A r ( u [ n ] , θ ) ∈ R 3 × 3 and B r ( u [ n ] , θ ) ∈ R 3 × 1 respectively denote the system matrix and in p ut matrix fo r the n th repetition time. No te that in (1), we implicitly assume th at all th e subvoxels shar e the same set of parameters θ . Moreover, assuming that the imaging experiment starts from thermal equilibirium , the in i- tial condition fo r (1) is given by M r [0] = [0 , 0 , M 0 ( r )] T , where M 0 ( r ) = M 0 / N v , M 0 denotes th e magnitud e o f th e magnetizatio n fo r the voxel △ at thermal equ ilibrium, and N v denotes the number o f sub-voxels in △ . Next, we describe th e o bservation equ ation. Note that the magnetizatio n is m easured at the co rrespond ing echo tim e ( i.e., T E n ) after each RF excitation , and th e signal detected by the receiver coil is propor tional to th e transverse c o mpon ent of the magnetization . Denoting m r [ n ] ∈ R 2 as the transverse magnetizatio n associated with the iso c h romat at △ r at the n th echo time, we have m r [ n ] = C r ( u [ n ] , θ ) M r [ n − 1] , (2) for n = 1 , · · · , N , where C r ( u [ n ] , θ ) ∈ R 2 × 3 denotes the output matrix. Summing up the magn etizations o f all th e isochroma ts, we can o btain th e transverse magnetization f or △ as m [ n ] = X r : △ r ⊂△ m r [ n ] . (3) Putting tog ether (2) an d (3), we have m [ n ] = X r : △ r ⊂△ C r ( u [ n ] , θ ) M r [ n − 1] . (4) The equ a tions (1) and (4) tog ether form a state-space model, which can describe spin d ynamics fo r various MR ﬁngerpr in ting imagin g sequences (e.g ., [1], [15], [16], [ 18]). Note that this mo del is n onlinear an d time -varying in nature. 2) Examp le S e quence: A s an example, we illustrate the above state-space model using a widely-used MR ﬁng er- printing sequ ence, i.e. , inv ersion recovery fast im aging with steady-state precession (IR-FISP) seq uence [15]. Th is se- quence is robust to off-resonance effects with the use of spoiler gra d ients [15]. Here the unkn own parameters are θ = [ T 1 , T 2 , M 0 ] T , and the acquisition parame ters are u [ n ] = [ α n , φ n , T E n , T R n ] T . In the IR-FISP sequ e nce, th ree physical processes drive the magnetizatio n e volutions: ( 1) RF excitation; (2) spin relax- ation; and (3) spin d ephasing . In the f ollowing, we u se the above isoch romat-su mmation appr oach to m odel spin dyn a m - ics. Speciﬁcally , und er Cartesian coordinates in the rotating frame [37], 1 we can form the system m atrix A r as A r ( u [ n ] , θ ) = G ( β r ) R ( T 1 , T 2 , T R n ) Q ( α n , φ n ) , (5) where Q ( α n , φ n ) ∈ R 3 × 3 models the RF excitation , i.e. , Q ( α n , φ n ) =   cos( φ n ) sin( φ n ) 0 − sin( φ n ) cos( φ n ) 0 0 0 1     1 0 0 0 cos( α n ) sin( α n ) 0 − sin( α n ) cos( α n )     cos( φ n ) − sin( φ n ) 0 sin( φ n ) cos( φ n ) 0 0 0 1   ; R ( T 1 , T 2 , t ) ∈ R 3 × 3 models the spin re laxation, i.e., R ( T 1 , T 2 , t ) =   e − t/T 2 0 0 0 e − t/T 2 0 0 0 e − t/T 1   ; and G ( β r ) ∈ R 3 × 3 models the spin d ephasing, i.e., G ( β r ) =   cos( β r ) sin( β r ) 0 − sin( β r ) cos( β r ) 0 0 0 1   , where β r denotes the pha se dispersion a ssoc ia ted with th e isochroma t at △ r . No te th at in the IR-FISP seq uence, β r is dominated by the impact o f spoiler grad ients, althoug h th ere are various o ther factors (e.g., diffusion effects) that can also contribute to spin d ephasing . Moreover , we can form the inpu t matrix B r as follows: B r ( u [ n ] , θ ) = M 0 ( r ) b ( T 1 , T R n ) , (6) where b ( T 1 , t ) =  0 , 0 , 1 − e − t/T 1  T . Note th at B r ( u [ n ] , θ ) models the recovery of the lon gitudinal magnetizatio n. Finally , we can f orm the output matrix C r as C r ( u [ n ] , θ ) = PR ( T 1 , T 2 , T E n ) Q ( α n , φ n ) , (7) 1 Alterna tiv ely , the model can be describe d in other coordinate systems (e.g., the polar coor dinates [22]); howe ver , note that the Cram ´ er-Ra o bound does not change under an y transform of the coordinate system for the Bloch equati on. IEEE TRANSACTIONS ON ME DICAL IMA GING 4 where P ∈ R 2 × 3 is the projectio n matrix tha t extracts th e transverse m agnetization , i. e., P =  1 0 0 0 1 0  , and R ( T 1 , T 2 , T E n ) an d Q ( α n , φ n ) ar e deﬁned in the same way as befo re. Here C r ( u [ n ] , θ ) m odels the RF excitation as well as the spin r elaxation. Note that this matrix does not account for spin d ephasing , since a spoiler grad ient is p laced after the echo time within each rep e titio n time. B. Cram ´ er -Ra o Boun d W e proceed to describ e the data mo del, with which we cal- culate the CRB an d perf orm experiment design . Note that the spin dy namics d escribed ab ove are utilized to perf orm contrast encodin g in MR ﬁngerp rinting experiments. Besides co ntrast encodin g , the d ata generatin g pr o cess also en compasses spatial encodin g and noise contamin ation. In [7], [8], we described a d ata mo del for this data gen e rating process. In princip le, we can use it to c a lc u late the CRB and p erform experime n t design. In practice, h owever , this is often com putationa lly very expensiv e with the use of the non -Cartesian Fourier tran sform in spatial en coding . Here we describ e a pr actical approa c h, in which we ignor e the spatial encod in g, and use th e following simpliﬁed data mode l: s [ n ] = m [ n ] + z [ n ] , (8) for n = 1 , · · · , N . Here s [ n ] ∈ R 2 is a vector that contains the magnetizatio n s collected at the n th echo time, and { z [ n ] } N n =1 denotes inde p endent, identica lly d istributed Gau ssian noise with z [ n ] ∼ N ( 0 , σ 2 I ) . Note tha t (8) correspon ds to the data model for a single voxel nuclear magn etic reso nance (NMR) experiment. Alterna tively , it can also be v iewed as the data model for a fu lly -sampled imag ing exp e riment, 2 in which there is no “crosstalk” be tween ma g netization ev olution s at different voxels. Despite such simpliﬁcation, we will demon stra te later that the exper iment design with (8) can be very effectiv e for highly- u ndersam p led M R ﬁngerp rinting experiments. Next, we derive the CRB f or the data mo del (8). From estimation th eory , the CRB provides a lower bound on the covariance o f any unbiased estimator under mild regula r ity condition s, and this bound can b e asymp totically achieved by the ML e stima tor [14]. Mathematically , the CRB can be expressed as the following informatio n ineq uality [14]: E   θ − ˆ θ   θ − ˆ θ  T  ≥ V ( θ ) = I − 1 ( θ ) , (9) for any unbiased estima to r ˆ θ , where I ( θ ) ∈ R p × p denotes the Fisher in f ormation matrix (FIM) deﬁned a s I ( θ ) = E "  ∂ ln p ( { s [ n ] } ; θ ) ∂ θ   ∂ ln p ( { s [ n ] } ; θ ) ∂ θ  T # , V ( θ ) ∈ R p × p denotes the CRB ma tr ix, and ln p ( x ; θ ) d enotes the log-likeliho o d fu nction of th e observation x param e terized 2 More precisely , this is equi v alent to a single-cha nnel N yquist-sampled Cartesi an Fourier acquisition with a discrete Fouri er transform based image reconstru ction. by θ . In (9), the matrix inequality A ≥ B means th at A − B is a p o siti ve semideﬁn ite m a trix. Note that bo th the CRB and FIM depen d on the underlyin g tissue par ameter θ , given that the data mode l ( 8) is nonlinear with respect to θ . M o reover , we can obtain the bound on the variance of individual tissue parameter estimate by extracting th e co rrespon ding diagon a l entry of the CRB matrix, i.e. , V a r  ˆ θ i  ≥ [ V ( θ )] i,i . (10) T o calculate the CRB in (9), we n eed to co m pute the FIM I ( θ ) . For the ad ditiv e Gaussian data mod el in (8), th e FIM has the particu larly simple form, which can be written as follows [14]: I ( θ ) = 1 σ 2 N X n =1 J T n ( θ ) J n ( θ ) , (11) where J n ( θ ) = ∂ m [ n ] / ∂ θ ∈ R 2 × p is th e Jacobian matrix. Finally , we describe the calculation of J n ( θ ) fo r (11). Giv en the state- space mo d el in ( 1) an d (4), su ch calcu latio n is equ iv alent to iterating a set o f difference equations. More speciﬁcally , notin g that J n ( θ ) = ∂ m [ n ] ∂ θ = h ∂ ∂ θ 1 m [ n ] · · · ∂ ∂ θ p m [ n ] i , (12) we can compute ∂ m [ n ] /∂ θ i for eac h entry o f θ . This ca n be done as follows. First, we take the derivati ve with respec t to θ i on bo th sides o f (4), which yie ld s ∂ m [ n ] ∂ θ i = X r : △ r ⊂△ ∂ C r ( u [ n ] , θ ) ∂ θ i M r [ n − 1] + X r : △ r ⊂△ C r ( u [ n ] , θ ) ∂ M r [ n − 1] ∂ θ i , (13) for n = 1 , · · · , N . Th en w e inv oke the der i vati ve with respect to θ i on bo th sides o f (1), which yie ld s ∂ M r [ n ] ∂ θ i = ∂ A r ( u [ n ] , θ ) ∂ θ i M r [ n − 1] + A r ( u [ n ] , θ ) ∂ M r [ n − 1] ∂ θ i + ∂ B r ( u [ n ] , θ ) ∂ θ i . (14) Now we can iter ate the two difference equ ations (13) and (14) to calculate ∂ ∂ θ i m [ n ] . Note that the initial condition s are g iv en by M r [0] = [0 , 0 , M 0 ( r )] T and ∂ M r [0] /∂ θ i = [0 , 0 , ∂ M 0 ( r ) /∂ θ i ] T . For the sake of concr eteness, we illus- trate the above pro cedure with the exam p le IR-FISP seque nce in the Append ix. C. Optimal Experimen t Design Giv en that the CRB provides a lower bound on the sm a llest possible variance f or any unb iased estimator, it can be used to charac ter ize the SNR efﬁciency o f an imagin g exper iment. This help s und e rstand the potential reliability of an MR ﬁn ger- printing expe riment, and ﬁgure out h ow much acquisition time is necessary to ach iev e a certain level of q uantitative accuracy . More imp ortantly , we c a n u se the CRB as a princip led tool to op timize the en coding pro cess of an M R ﬁngerp rinting experiment. For example, gi ven a set of representative MR IEEE TRANSACTIONS ON ME DICAL IMA GING 5 tissue parameter s { θ ( l ) } L l =1 , we cou ld o ptimize the data ac- quisition par ameters of an MR ﬁngerpr inting exper iment to maximize its SNR efﬁciency . Mathematically , such an optima l experiment d esign pr oblem can b e fo rmulated as fo llows 3 : min u L X l =1 Ψ  V ( θ ( l ) )  s.t. u ⊂ U , (15) where Ψ( · ) den otes the design criterion , which is a scalar function o f the CRB matr ix; u = [ u [1] , · · · , u [ N ]] ∈ R q × N denotes the acq uisition paramete r s for an M R ﬁngerp rinting experiment with N tim e p oints; and U ⊂ R q × N denotes the constraint set for f easible data acq uisition pa r ameters. Note that, for the sake of co ncrete illustration , we assume that the total n u mber of time po in ts, i.e., N , is g i ven in (15). In pr actice, N can be speciﬁed acco rding to the desired experiment du ration and th e constraints o n TRs. Mo reover , giv en that the CRB m atrix depen ds on th e u nderlyin g tissue parameter, we assume tha t, for a speciﬁc applicatio n of interest (e.g., neuro imaging ) , we have the k nowledge of the r ange of MR tissue parameter values p rior to ou r experimen t design. While it is d esirable to design experimen ts that are u niversally optimal fo r all possible param eters with in th e ran ge, this is often not f e asible. As such, we select a few representative tissues as a p ractical com promise. As a n example, we specia lize (15) to optimize MR ﬁnger- printing experimen ts with the I R-FI SP sequenc e . First, we specify the desig n criter io n for Ψ( · ) . No te that there are various infor m ation criteria that can be used, in cluding the A-optimality , D-optim ality , and E-op timality cr iteria (see [13] for a compr e hensive survey). Here we cho o se the A-o ptimality criterion [1 3], which minimize s the tra c e of the CRB matrix (i.e., the total variance of tissue pa r ameter estimates). Fur ther, we incorpo rate weigh tin gs into the design criter ion, wh ich is motiv ated by: (1 ) the CRBs for different tissue p arameters are often at very different scales; and (2) we may w ant to tailor a design to the parameter s that ar e most r elev ant to speciﬁc applications of interest. Accord ingly , we have Ψ( · ) = tr ( WV ( θ )) , where W is a diagona l matr ix whose entries conta in weighting s fo r different tissue param eters. Second, we specify the data acquisition param e te r s u = [ u [1] , · · · , u [ N ]] . Note that the acq u isition par ameters for the IR-FISP seque n ce in clude the ﬂip an gles, RF pu lse phases, echo times, and rep etition tim es. Following the early work [15], we assume the ﬂip angles and rep etition times to be the design parameters, wh ile ﬁx ing th e RF p ulse phases an d echo times. Acc o rdingly , we h ave u [ n ] = [ α n , T R n ] T , for n = 1 , · · · , N . Lastly , we specify the constraint set U for the acquisition pa- rameters. T aking into acco unt various ph ysical consideratio n s (e.g., speciﬁc absorption rate (SAR), and/or total acquisition time), we im pose upper bou n ds an d lower bou nds f or the 3 The term “optimal experimen t design” refers to the goal of the problem formulati on, i.e., maximizing the SNR ef ﬁcienc y of an experiment . The use of this term follo ws the con vent ion in statistic s [12], [13], [40] and in MR imaging [41]–[43], and is unrelated to whether a speciﬁc solution algorithm produces a globally optimal design or not. acquisition parame ters, i.e. , T R n ∈  T R min n , T R max n  , 4 and α n ∈  α min n , α max n  , for n = 1 , · · · , N . Ac c ording ly , we can formu late the optimal experim ent design prob lem as follows: min { α n ,T R n } P L l =1 tr  WV  θ ( l )  s.t. T R min n ≤ T R n ≤ T R max n , 1 ≤ n ≤ N , α min n ≤ α n ≤ α max n , 1 ≤ n ≤ N . (16) Besides the physical constrain ts, it is often useful to take into accoun t other constrain ts from image deco d - ing/recon struction, especially when dealing with highly- undersam pled MR ﬁngerp rinting experimen ts. While different reconstruc tion metho ds may use different strategies for image decodin g (from highly- undersam pled data) , th ey often bene ﬁt from magne tiza tion evolutions being smoothly varying. 5 There are a nu mber o f ways to enforce th is p roperty . Here we inc o r- porate an add itional set of constraints into (16), which restricts the max imum ﬂip ang le variations betwee n con secutive time points. As will be demo nstrated later, such constraints are effecti ve in p romotin g smooth mag netization evolutions, which can yield better perfo rmance fo r imag e reconstructio n with highly- u ndersam p led data. Accord ingly , we re f ormulate the optimal exper iment desig n pro b lem as follows: min { α n ,T R n } N n =1 P L l =1 tr  WV ( θ ( l ) )  s.t. T R min n ≤ T R n ≤ T R max n , 1 ≤ n ≤ N , α min n ≤ α n ≤ α max n , 1 ≤ n ≤ N , | α n +1 − α n | ≤ ∆ α max n , 1 ≤ n ≤ N − 1 . (17) where ∆ α max n speciﬁes th e maxim um ﬂip angle variations between α n and α n +1 . The p roposed for mulations in (1 6) and (17) re sult in n on- linear and no nconve x o ptimization pr o blems. A numb er of numerical algor ithms can b e employed to solve the opti- mization prob lems, including standard n o nlinear optimization methods [4 6] and stochastic optimization meth ods [ 4 7]. As an example, we use a state-of-the - art nonline a r optimization method, i.e., seque n tial qu a dratic p rogram ming ( SQP) [ 46], to seek lo cal m inima for ( 16) a nd (17). T h e SQP algorithm is an iterative algo rithm, and at each iteratio n, it perf orms a quadra tic appro ximation of th e cost function at the cu r rent 4 Alterna tiv ely , we could replac e this constraint with a constraint on the total acquisition time as in our pre vious work [19], [20 ]. In principle, such a formulati on could allo w larger variat ions of T Rs, which could potential ly improv e the CRB of the optimized expe riments. Howe ver , note that this formulati on has to simultane ously determine both N and { α n , T R n } N n =1 , which often leads to a highly non-con vex (or mixed intege r) optimiz ation problem. Wit hout loss of general ity , we hav e decid ed to focus this paper on a simpler optimizat ion formulation that av oids this issue. 5 For example , for the conv entiona l reconstruction method [1] that utilizes direct pattern matching, smoother magnetiz ation evol utions often have less correla tion with noise-like aliasing artifa cts. This causes the underlying magnetiz ation ev olution to be better dif ferentiat ed from aliasing artifact s, leadi ng to improv ed accurac y in dictionary matching. As another example, for the low-ra nk/subspace reconstruc tion method [11 ], [44], [45], smooth magnetiz ation e volutio ns often result in a sm aller low-ra nk approximati on error giv en the same rank va lue. This in turn enables better reconstructio n accura cy . Finally , for statistical reconstruction methods [8], [45] that in volv e solving noncon ve x optimizati on problems, an initi alizati on from the impro ved patte rn-matching recon struction or lo w-rank reconstruc tion method ofte n leads to better performance . IEEE TRANSACTIONS ON ME DICAL IMA GING 6 solution, and also linear izes the con straints. Then it solves a constraine d quad ratic optimizatio n p r oblem to update the solution. As with other nonlin e ar optimiz a tion method s, the perfor mance o f th e SQP a lg orithm is gener ally depen dent on initialization. Her e we in itialize the algorithm with th e ac- quisition p arameters f rom the conventional MR ﬁnge r printing experiment [15], with which the algo rithm consistently yields good perfo rmance, althou gh other initialization schemes (e.g ., a mu lti-start strategy) may lead to b etter perf ormanc e . I I I . R E S U LT S In this sectio n, we show representative results fr om numeri- cal simulations, phanto m experim ents, and in vivo exper iments to illustrate the p erform ance of the pro posed a p proach . A. S imulations 1) General S etup: W e created a nu m erical brain phan tom to simulate single- channel MR ﬁngerp rinting experiments. W e took th e T 1 , T 2 , an d M 0 maps from the Brainweb d atabase [48] as the grou nd truth, as shown in Fig. 1. W e set the experimental ﬁeld-of-v iew (FO V) as 300 × 300 mm 2 , an d the matrix size as 25 6 × 256 . W e simulated MR ﬁngerpr in ting ex- periments with th e IR-FI SP sequ e nce [1 5], which is r obust to main mag netic ﬁeld inh omog eneity . Mo reover , for simp licity , we assum ed that the transmit RF ﬁeld was homog eneous in the simu lations. W e per f ormed Bloc h simu lations to gen erate contra st- weighted images. Here we used th e isochrom a t-summation approa c h [38], [39], in which we simulated magnetization ev olutions with 400 isochr omats f o r each v oxel. In Bloch simulations, w e c o nsidered th r ee different sets of acquisi- tion p arameters: (1) the conventional scheme [15], (2) the optimized scheme with (16), and (3) the optimized scheme with (17). In Fig. 2 , we sh ow the acquisition p arameters from the con ventiona l scheme with N = 100 0 , as well as the resulting ma g netization evolutions for the two r egions-of- interest (R OIs), respec ti vely , in the white m atter tissue a nd a gray m a tter tissue ( a s marked in Fig. 1 (a)). W e gener a ted k -space data from the contrast-weig hted images using the non-u niform Fourier tran sform [4 9]. I n the simulations, we con sidered two setups of MR ﬁngerp rinting experiments: • Fully- sampled experiment: we acquired fu lly-sampled Cartesian k -space data for each time point (or TR in d ex), as in [50]. • High ly-und ersampled experimen t: we acquired h ighly- undersam pled spiral k -space data a t each time poin t, as in [15]. For this, we used the sam e spiral trajector y as [15], and acqu ired only one spiral interleaf for ea c h time point, whereas a full set o f spiral tr ajectory consists of 48 interle aves. W e add e d complex white Gaussian noise to the measured k - space data accor ding to the p re-speciﬁed no ise level σ 2 . Her e we deﬁne the signal-to-n oise ratio as SNR = 2 0 log 10 ( s/σ ) , where s denote s the average value of M 0 in a region of white matter . T his de ﬁn ition measures the SNR in decib e ls (dB). W e perfor med the ML reco nstruction [7], [8] for th e above experiments. Note that fo r th e fully-samp led Cartesian experiments, the ML appr oach is equivalent to the d irect Fourier reconstru ction, followed by the d ictionary - based pat- tern m atching [51]. For the highly- under sam pled experimen ts, we so lved th e r econstructio n problem with the algorithm in [7], [8], which we initialized with the g r idding recon struction. Her e the dictionary used in the ML reconstru ction was constructed based on the following p arameter discretization scheme: we set the T 1 value in th e range [20 , 30 00] ms , in which we used an increment of 10 ms for [20 , 15 00] ms and an incremen t of 30 ms for [1501 , 3000 ] ms ; we set the T 2 value in th e range [30 , 500 ] ms , in which we used an incr ement o f 1 m s for [30 , 200] ms and an incremen t of 5 m s for [2 01 , 50 0] ms . T o assess the reconstructio n accu r acy , we u sed th e following two metrics: (a) overall error, i.e., k I − ˆ I k 2 / k I k 2 , where I and ˆ I respectively d enote the tru e parame te r m a p an d reco n structed parameter m ap, and (b) voxel wise relative er ror, i.e., | I v − ˆ I v | / | I v | , whe re I v and ˆ I v respectively den ote the values o f I and ˆ I at the v th voxel. 2) Implem e n tation of (16) and (17) : Her e w e describe the detailed imp lem entation of the p ropo sed appr o ach fo r the above a p plication example. For convenience, we r e fer to the optimized schemes with (16) and (1 7) as Optimized- I and Optimized-I I, respectively . I n this work, we assumed that T 1 and T 2 were o f primary interest, an d chose thr ee represen - tati ve tissues, i.e., θ (1) = [700 ms , 60 ms , 0 . 6] , θ (2) = [850 ms , 5 0ms , 0 . 6] , and θ (3) = [1100 ms , 102 ms , 0 . 6 ] , for ( 16) and (17). 6 Moreover , we m anually chose the weightin g matrix W = diag ([2 . 0 × 10 − 5 , 5 . 0 × 1 0 − 4 , 3 . 0 × 1 0 1 ]) for b oth (16) and ( 1 7) to en su re the goo d perfo rmance of the optimized experiments. Further we specify the constraints for the acqu isition p a- rameters in ( 16) and (17). Spec iﬁcally , we set the maximum ﬂip an g le a s α max n = ( 180 ◦ , if n = 1 , 60 ◦ , if 2 ≤ n ≤ N . Note that this allows an 18 0 ◦ in version pulse imposed at the beginning o f imaging experimen ts, which is ofte n ad vanta- geous for the T 1 estimation [20]. W e set the minimum ﬂip angle as α min n = 10 ◦ for 1 ≤ n ≤ N . W e respectively set the maximum and minim um repetition times as T R max n = 15 ms and T R min n = 11 ms fo r 1 ≤ n ≤ N . Here, it is worth mentionin g that the ab ove co nstraints rough ly match the range of the acqu isition par ameters in the conv ention al scheme . For (17), we have th e additio nal constraints on the ﬂip ang le variations, which were set as ∆ α max n = ( + ∞ , if n = 1 , 1 ◦ , if 2 ≤ n ≤ N − 1 . Note th at the ab ove co nstraint does not r estrict the ﬂip a n gle variation b e twe en the inversion p ulse and the second RF pulse. 6 The T 1 and T 2 v alues of these tissues were arbitrarily chosen from the range of tissue parameter val ues that are rele v ant to neuroimagin g. T o avoid the inv erse crime, we did not take tissue paramete r val ues directly from the brain phantom. IEEE TRANSACTIONS ON ME DICAL IMA GING 7 Fig. 1. Ground truth parameter maps for the brain phantom: (a) T 1 map, (b) T 2 map, and (c) M 0 map. Note that the two R OIs (respect iv ely in the white matter and gray matter) are m arked in (a). Fig. 2. Acquisition parameter s of the con ventiona l scheme with N = 1000 as well as the resulting m agneti zation e voluti ons. (a) Flip angle train. (b) Repeti tion time train. (c) Magnetizat ion ev olutions for the white matter and gray matter R OIs marked in Fig. 1 (a). Note that the ﬁrst RF pulse is 180 ◦ , which excee ds the scale of the vertica l axes in (a). W e applied the SQP algo rithm to so lve the op timization problem s associated with (16) and ( 17), and in itialized the al- gorithm with the acq u isition para m eters { α n , T R n } N n =1 from the co n ventiona l scheme. W e termin ated the algorith m , w h en the cha nge of th e solution was less than the pre-spec iﬁed tolerance (i. e ., ǫ = 1 e − 4 ) or the max im um iteration (i.e., J max = 5 × 10 4 ) was reach ed. The runtime o f the alg orithm depend s on the length of acq uisition. For example, with N = 400 , solving (1 6) and (17) respectively took about 29 0 min and 140 min o n a Linux workstation with 2 4 Intel Xe on E5-26 4 3, 3.4 0 GHz pr o cessors and 128 GB RAM ru nning Matlab R201 5b. W e optimiz e d th e acquisition paramete r s ind ependen tly for se veral ch o ices of N ( i.e., N = 30 0 , 400, 500 , 600 , 7 0 0, and 800). As an example, Fig. 3 shows th e op timized acquisition parameters with N = 400 fro m Optimized-I and Optimized-II . As can be seen, th e o ptimized acquisition parame te r s appear to be hig hly structur ed, which are remark ably different from the acqu isition parameter s fr om the co n ventional sche m e. In particular, th e optimized repetition times turn out to be binar y (i.e., switch in g between T R max n and T R min n ). Besides N = 400 , we also had similar observations for all o ther acquisition lengths. W e show one mo re e xam ple, i.e. , the optim ized schemes w ith N = 60 0 , in the Sup plementar y Ma terial. In Section IV, we will discuss this interesting o bservation. Fig. 3 also shows the resulting magn etization ev olutions from Optimized-I and Optimized-II. As can b e seen, the magnetizatio n evolutions from Op timized-I exhibit signiﬁcant oscillation d ue to the dramatic chan ge of acquisition p a- rameters (with in the ﬁrst 100 time points). I n co ntrast, by enforcin g the a d ditional constraints o n the ﬂip ang le variations in Optimized-II, the o scillation behavior has been signiﬁcantly suppressed, and the resultin g magnetiza tio n evolutions beco me much smoother . Later we will demon stra te that this is advan- tageous for ima ge rec onstruction from hig h ly-und ersampled data. 3) Evalu ation of CRB: W e ev aluated the CRB f or the con- ventional scheme, Optimized-I , an d Optimized- II. Note th at the CRB is a lower bo und on the variance (or equ iv alently , the mean-squ are error) o f an u nbiased estimator, w h ich c h aracter- izes th e SNR prop erty of a n im aging exper iment. Speciﬁcally , we c a lculated the CRB associated with the three sets of acquisition parameter s over different acqu isition length s. For all the experimen ts, we set SNR = 33 dB . I n or der to sh ow the CRB f or b o th T 1 and T 2 at th e sam e scale, we u sed the normalized CRB deﬁned as nCRB i = p CRB( θ i ) / θ i . Fig. 4 shows the norma lized CRB versus the acquisition length for the wh ite matter R OI. As expected, for all three schemes, the n CRBs of T 1 and T 2 reduce as th e acquisition IEEE TRANSACTIONS ON ME DICAL IMA GING 8 Fig. 3. T wo optimized acqui sition schemes with N = 400 , and the resulting magnetiza tion e voluti ons for the white m atter and the gray matter R OIs marked in Fig. 1 (a). (a) Flip angle train in Optimiz ed-I. (b) Repeti tion time trai n in Optimized -I. (c) Magnetizati on ev olutions from Optimize d-I. (d) Flip angle train in Optimize d-II. (e) Repetition time train in Optimized-II. (f) Magneti zation evo lutions from Optimize d-II. Note that the ﬁrst RF pulses are 180 ◦ for both Optimize d-I and Optimized-II, w hich exceed the scale of the vert ical axes in (a) and (d). Fig. 4. N ormalized CRB versus acquisit ion length N for the R OI in the w hite matter tissue. (a) Normalized CRBs for T 1 . (b) Normalized CRBs for T 2 . length becom es lon ger . For T 1 , the nCRB r apidly appr oaches its asym p totic limit with short a cquisitions for all thr ee schemes. However , for T 2 , the conventional scheme needs much longer acq u isition to attain a go od nCRB . This is consistent with th e previous o bservations in [7]–[11], and indicates that from an estimation-theo retic pe r spective, there is sign iﬁcan t roo m for improvement over the co n ventional scheme. In contrast, th e two optim ized schemes signiﬁcantly improve the nCRB o f T 2 , especially with shor t acquisition lengths. For example, the optimized experiments reduc e the nCRB o f T 2 by abo ut a factor o f two fo r N = 400 . Also note th at the nCRB f or Optimized-I is better than for Optimized-I I over all the acqu isition len gths. This is as expected, since, with a smaller set of constraints, Optimized-I searches over a larger feasible space o f ac q uisition parame - ters. No netheless, the difference between the two optimized schemes is quite small. In the Supp lementary Material, we also show the CRB ev aluation with respect to the gray matter R OI (mar ked in Fig. 1 (a)), from which we had similar observations. 4) Evalu ation of fully-sample d experiments: W e evaluated the fully-samp led MR ﬁngerprin ting experiments described in Section II I.A. Note that this scenario exactly correspo nds to the data model in (8), w ith which we calcu lated the CRB and perf ormed experimen t desig n . Speciﬁcally , w e simulated the experimen ts at N = 4 00 and SNR = 33 dB, using the conv ention al sch e me, Optimized-I, and Optimized-I I. Fig. 5 shows the reconstru cted T 1 and T 2 maps. As can be seen, the two o p timized schemes improve the accu r acy of b oth T 1 and T 2 , altho ugh the improvemen t is more signiﬁcan t for T 2 . Moreover , th e performa nce of Optimized-I is slightly better than tha t of Optimized- II, c onsistent with the CRB prediction shown b efore. W e furth er in vestigated the bias-variance p roperty of the reconstruc ted parameter maps. Speciﬁcally , we perf ormed Monte Carlo (MC) simulations with 100 trials to calculate the bias, stand ard deviation, and roo t-mean-sq uare erro r of the re- constructed p arameter maps. For convenience, we norm a lized these quan tities as f ollows: (1 ) normalize d bias: Nbias v = ˆ E h    I v − ˆ I v    i / I v , (18) IEEE TRANSACTIONS ON ME DICAL IMA GING 9 Fig. 5. Reconstruc ted parameter maps from the fully-sampl ed MR ﬁngerprinting exper iments ( N = 400 and SNR = 33 dB ), using the acquisition paramete rs from the con ventio nal scheme, Optimize d-I, and Optimized -II. (a) Reconstruc ted T 1 maps and associated relati ve error maps. (b) Reconstructed T 2 maps and associate d relati ve error maps. Note that the overa ll error is labele d at the lower right corner of each error map, and the regions associate d with the backgrou nd, skull, scalp, and CSF were set to be zero. Fig. 6. Bias-var iance analysis of the reconstruc ted paramete r maps from the fully-sample d MR ﬁngerprinti ng experi ments ( N = 400 and SNR = 33 dB ), using the acquisiti on parameters from the con ventiona l scheme, Optimize d-I, and Optimized-II. (a) Normalized bias, standar d devi ation, and root-mea n-square error for (a) T 1 maps and (b) T 2 maps. The regions associa ted with the background, skull, scalp, and CSF were set to be zero. (2) no rmalized standard d eviation: Nstd v = s ˆ E     ˆ I v − ˆ E ( ˆ I v )    2  / I v , (19) and (3 ) norm alized ro ot-mean - square erro r: NRMSE v = s ˆ E     I v − ˆ I v    2  / I v , (20) where ˆ E ( · ) de notes the empir ical me an ev aluated for the MC simulations, and I v and ˆ I v respectively den ote the v th voxel IEEE TRANSACTIONS ON ME DICAL IMA GING 10 Fig. 7. Overall error versus the acquisition length N for the fully-sampl ed MR ﬁngerprinti ng expe riments with SNR = 33 dB . (a) Overa ll error of T 1 map. (b) Overa ll error of T 2 map. Note that the overal l error is calcul ated with respect to the whole brain, excluding the s kull, scalp, and CSF . Fig. 8. Overa ll error versus the SNR le vel for the fully-sampl ed MR ﬁngerprint ing experiment s with N = 400 . (a) Overal l error of T 1 map. (b) Overa ll error of T 2 map. Note that the ov erall error is calcula ted with respect to the whole brain, excludi ng the skull, s calp, and CSF . from the true param eter map and reconstructed p arameter map. Note that NRMSE v = q Nbias 2 v + Nstd 2 v . Fig. 6 shows the nor malized bias, standard deviation, and root-me a n-square er ror maps for the reconstruc te d T 1 and T 2 maps. As can be seen, Op timized-I an d Optimized-I I reduce the n ormalized standard d eviations f or both T 1 and T 2 , com- pared to the co n ventional schem e. Co n sistent with the CRB prediction and the resu lts shown in Fig. 5, the improvement for T 2 is mo re substan tial than fo r T 1 . M o reover , for all three acquisition schemes, the norm alized standard deviation is much larger than the nor malized b ias, and the n o rmalized roo t- mean-squ are erro r is do minated b y the norma lize d standard deviation. The ab ove behavior ca n be expected , given that th e ML r e c onstruction is asym ptotically unb iased [ 14]. Moreover , we evaluated the fu lly-sampled experiments with different acqu isition lengths, with N rangin g fr o m 300 to 800 . Here we set SNR = 33 dB fo r th e exper iments. Fig. 7 shows the overall errors of T 1 and T 2 versus the acq uisition length. Clearly , the two optimized acqu isition schem es outp erform the conv ention al scheme over all the acquisition len gths. A s one more example, we show the rec onstruction results fo r N = 600 in the Su p plementar y Material. Finally , we ev aluated th e fu lly-sampled experimen ts with different SNR lev els, ran ging f rom 2 8 dB to 38 dB. He r e we set at N = 4 00 f or all the experiments. Fig. 8 shows the overall errors of T 1 and T 2 versus the SNR. As can b e seen, the optimized acquisition schem es outp e rform the conventional scheme over all the SNR levels. 5) Evalu ation of highly- undersampled experiments: W e repeated the same ev aluation s but applied to the h ighly- undersam pled ca se. Note that this more closely matches the way that MR ﬁnge rprinting is ap plied in p ractice. Fig. 9 shows the recon structed T 1 and T 2 maps fro m the hig hly- undersam pled exper iments at N = 400 and SNR = 33 d B, u s- ing the co n ventional sch eme, and the two optimized schem es. As can be seen, Optimized -I im p roves th e accuracy of the T 2 map over the conventional scheme, but at the expense of degrading th e accu racy of the T 1 map. In con trast, Optimized- II pr ovid es better accu racy for both T 1 and T 2 maps, which is h ighly desirable. Note that the ML reco nstruction inv olves solving a non linear an d n onconve x optimization p r oblem, for which a g ood in itializatio n is o ften req u ired. By enf orcing the additional c o nstraint on the ﬂip angle variations, Optimized-I I results in mu ch smoo ther m agnetization evolutions (as shown in Fig. 3). W ith th e highly -under sampled data, this often leads to better pattern matching results for th e conventional r econ- struction, which in turn provides an improved initialization for the ML reconstru ction. Fig. 10 sh ows the norma lize d bias, stand ard deviation, and root-me a n-square erro r maps fo r the re c onstructed T 1 and T 2 maps fr om the MC simulations (with 1 00 trials). Clearly , Optimized-I I r educes the n ormalized standard deviation and root-me a n-square error f o r both T 1 and T 2 maps, com pared to the conv ention al scheme. Mor eover , with smooth magneti- zation ev olutions, Op timized-II reduces the bias comp ared to Optimized-I . T his f urther illustrates the merit of introdu c in g the constraint on the ﬂip angle variations for the highly - undersam pled exper iments. Fig. 11 shows the overall error s of T 1 and T 2 versus the acquisition length, for the highly-un dersample d experiments with SNR = 3 3 dB. Clearly , Optimized -II o u tperfor ms the conv ention al scheme and Optimized- I for all the acq u isition lengths. As a furth e r illustratio n, we show the re construction results fo r N = 600 in the Su pplementa r y Material. Fig. 12 shows the overall error s of T 1 and T 2 versus the SNR, for the highly- undersam pled experiments with N = 40 0 . As can be seen, Optimized-II provides better a c c uracy than the conv ention al scheme and Optimized-I over all the SNR le vels. B. Ph antom Expe rimen ts W e evaluated the p ropo sed ap proach with phantom exper i- ments. He r e we foc us on the scenario of high ly -und e r sampled MR ﬁngerp rinting experimen ts, which is of the most practical interest for quantitative MR imaging. Speciﬁcally , we created a physical phanto m that consists of 9 plastic tubes, each one ﬁlled with a solution of Gadolinium an d Agar at d ifferent concentr a tions. T h is created different comb in ations of T 1 and T 2 values that are relevant to the neuroimag ing app lication [52]. W e carried out the exp e riments on a 3T Siemens Tim T rio scann er (Siemens Medic a l So lutions, Erlan gen, Germany) equippe d with a 32-chann el head array coil. The relev ant imaging parameters include: FO V = 300 × 300 mm 2 , matrix size = 256 × 256 , and slice thickn ess = 5 mm . IEEE TRANSACTIONS ON ME DICAL IMA GING 11 Fig. 9. Reconstructed parameter m aps from the highly-under sampled MR ﬁngerprinting experime nts ( N = 400 and SNR = 33 dB ), using the acquisiti on paramete rs from the con ventio nal scheme, Optimized-I, and Optimize d-II. (a) Reconstru cted T 1 maps and associated relati ve error maps. (b) Reconstructe d T 2 maps and associated relat iv e error maps. Note that the overal l error is labe led at the lower right corner of each error map, and the regio ns associated with the backgrou nd, skull, scalp, and CSF were set to be zero. Fig. 10. Bias-v arianc e analysis of the reconstru cted parameter maps from the highly-undersampl ed MR ﬁngerprinting expe riments ( N = 400 and SNR = 33 dB ), using the acquisiti on parameter s from the con ventio nal scheme, Optimized-I , and Optimized -II. (a) Normalized bias, standard deviat ion, and root- mean-square error for (a) T 1 maps and (b) T 2 maps. The regi ons associate d with the backgr ound, skull, scalp, and CSF were set to be zero. 1) General evalua tion: W e p erform ed three sets of ex- periments with N = 40 0 , r e spectiv ely , using the acqu isi- tion par a meters from the conventional schem e, Optimized -I, and Optimized -II. W e used the same spiral tr ajectory and sampling pattern as in the nume rical simulations. Here the acquisition times for the conventional schem e, Optimize d- I, an d Op timized-II were 5.2 8 sec, 5.2 4 sec, and 5.22 sec, respectively . T o evaluate th e perfo rmance of the above exper- iments, we also acq uired a set of r e ference T 1 and T 2 maps, IEEE TRANSACTIONS ON ME DICAL IMA GING 12 Fig. 11. Ov erall error ver sus the acquisition length N for the highly- undersample d MR ﬁngerprin ting experi ments with SNR = 33 dB . (a) Ove rall error of T 1 map. (b) Overal l error of T 2 map. Note that the overa ll error is calcu lated with respect to the whole brain, excludi ng the skull, scalp, and CSF . Fig. 12. Overall error versus the SNR lev el for the highl y-undersampled MR ﬁngerprint ing experiment s with N = 400 . (a) Overal l error of T 1 map. (b) Overa ll error of T 2 map. Note that the ov erall error is calcula ted with respect to the whole brain, excludi ng the skull, s calp, and CSF . by perf orming a fully-samp led MR ﬁnger printing experiment 7 with N = 1000 using the acquisition parame te r s from th e conv ention al scheme. Th e acquisition tim e for this exp eriment was ab out 18 m in. Additionally , we calibrated the spiral trajectory with a specialized pulse seque nce [5 3] to av oid the potential trajectory distortion (caused by eddy cu r rents an d gradient delay ). Finally , we perfor m ed an aux iliary scan with a gradien t echo (GRE) sequ ence, from which we e stimated the co il sensitivity ma p s. Th is acquisition to ok 1.28 sec. W e pe rforme d the ML recon struction f or the above exper- iments, and in corpo rated the slice-pro ﬁle co r rection [54] into the dictionar y . Fig. 13 shows the reconstru cted T 1 and T 2 maps, the relativ e error maps (evaluated with respect to the referenc e data) , and the correspo nding reconstru ction errors over each tu be. I t is clear that Optim ized-II sign iﬁcantly improves the accu racy of th e T 2 map over the conv ention al scheme, while pr oviding similar accur a cy f or the T 1 map. Compared to Optimized-I, Optimized- II also provides better accuracy , par ticularly for T 1 maps. This conﬁr ms the be neﬁts of in troducin g the constraint on the ﬂip angle variation for highly- u ndersam p led M R ﬁngerp rinting experiments. 7 The fully-sampl ed expe riment was performed by repeating a highly- undersample d acquisition 48 times. For each acqui sition, we switched to a dif ferent spiral interlea f at ev ery time point. Note that a short time delay was added between consecuti ve acquisit ions to ensure that the magnetizati on starts at thermal equilib rium. 2) Evalu ation o f cr oss-scan variance: Given that variance reduction is a direct beneﬁt of the CRB b a sed exper iment design [3 3], we evaluated the variance associated with the three acqu isition schem e s. Speciﬁcally , we co n ducted each acquisition 15 times, an d calculated the no rmalized standard deviation associated with each acq uisition scheme. W ith the absence of th e gro und tru th for the phanto m experiments, the standar d d eviations were n ormalized with respect to th e parameter m aps reco n structed from the fully- sampled data with N = 1000 . Fig. 14 shows the norm alized stand ard deviation m aps for the three acquisition sch emes, and Fig. 15 shows the norm a lized stand ard d eviation averaged over each tube. It is e viden t that the two optimized exp e r iments sub- stantially red uce th e standard deviation for T 2 maps, wh ile providing similar standard deviation for T 1 maps. In pa rticular, Optimized-I I achiev es a factor o f two re duction in the standard deviation of T 2 maps for all the tubes. C. In V ivo Exp eriments W e ev aluated the per forman ce of the p roposed appro ach f o r in viv o experiments. W e con ducted th e im aging experiments on a healthy v olun teer with the ap proval fr om the local Institutional Re view Boar d and the informed consent was obtained fro m the subject. W e per f ormed highly- u ndersamp led MR ﬁngerp rinting exper iments with the co nventional scheme, Optimized-I , and Optimiz e d-II on the same scanner with the acquisition len gth N = 40 0 . W e perfor med th e M L reconstruc tion to estimate the T 1 and T 2 maps from hig h ly- undersam pled data . Moreover , we repeated each acquisition 15 times, and evaluated the sample variance associated with each acqu isition scheme. Fig. 16 shows the recon structed T 1 and T 2 maps from the three acquisition schem e s, as well as the normalized standard deviation maps from the repe titio ns of the imaging experiments. Here the standar d deviations were all nor malized with respect to the p a rameter maps r econstructed from the con - ventional sche me with N = 40 0 . Fig. 17 sh ows the norm alized standard deviations averaged over the two representative R OIs respectively in the gray matter and white matter . Additio nal results comparin g the tissue par a meter stan dard deviations for the whole brain can be found in the Supplem entary Material. First, note tha t the two optimized acquisition schem es pro- vide th e T 1 and T 2 maps in a similar rang e as the conventional scheme, conﬁrming the fe a sibility of th e p roposed experiment designs. Seco nd, the two optim iz e d acqu isition sch emes im- prove th e stand ard d eviation fo r T 2 , wh ile p roviding similar standard deviation for T 1 . T he above results a re consistent w ith those for th e numer ical simulations and p h antom experimen ts. Lastly , comparin g with Optimized-I, Optimized-II reduces the artifacts associated with the T 1 reconstruc tion in the white matter , and also enables better variance for b oth T 1 and T 2 estimation. This fu rther illustrates the beneﬁts of intro ducing the ﬂip variation con straint in th e experimen t design. I V . D I S C U S S I O N It is worth discussing several r elated aspects as well as the po tential extensions of the pr o posed framework. First, IEEE TRANSACTIONS ON ME DICAL IMA GING 13 Fig. 13. Reconstructe d parameter maps for the phantom ex periments ( N = 400 ) using the acquisiti on parameters from the con ventiona l scheme, Optimized-I, and Optimized-II. (a) Reconstr ucted T 1 map and associated relat iv e error m ap. (b) Overall error of T 1 reconstru ction ov er each tube. (c) Reconstruc ted T 2 map and associated relati ve error map. (d) Overall error of T 2 ov er each tube. Note that the highly-un dersampled MR ﬁngerprinti ng experi ments using the con venti onal scheme, Optimized-I, and O ptimized-II respecti vely took 5.28 sec, 5.24 sec, and 5. 22 sec, whereas the fully-sample d MR ﬁngerprinting expe riment for the reference data took 18.3 m in. Fig. 14. Normalized standard de viation maps associate d with the con ventiona l acquisi tion, Optimized-I, and Optimized -II, estimated from 15 indepen dent imaging expe riments. (a) Normalized standard devi ation maps for T 1 . (b) Normalize d standard devia tion maps for T 2 . regarding the m odel of sp in dynamics, we introd u ced th e state- space mod el b ased on the isochro mat sum m ation appro ach [38], [39]. Alternatively , w e can derive the m o del from th e expanded p hase graph (EPG) for malism [1 5], [55]. Similar to the early work [3 9], we have observed that the two models are equivalent in terms of simulating ma gnetization ev olutio ns (see the Supplemen tary Ma terial for the num e rical ju stiﬁca- tion). Howev er, the state-space mod el from the isochrom at summation app roach is co n ceptually simple, a nd easy to describe math ematically , which facilitates the subsequent CRB deriv ation. Also note th at ou r work used a simpliﬁed data mo d el (8) to calculate the CRB and to optimize acquisition parameters. Th is yields very small-scale FIMs, enab ling efﬁcient computatio n and storage. T o tailor f or highly-u ndersamp led experiments, we f urther in corpor ated th e constraint on th e ﬂip a n gle vari- ation to b alance the imp rovement of encod ing efﬁciency a n d the con siderations for decod in g. A s a genera liza tion, it would be interesting to formulate the p roblem that d irectly opti- mizes high ly-und ersampled MR ﬁng erprintin g experiments. This could allow a jo int design of acq uisition parameter s and k -space trajectories (e.g., [56], [57]). Howev er, no te that the CRB calculation in this case r equires a signiﬁcan tly h igher- dimensiona l fo rmulation th at mo dels b oth co n trast encod ing and spatial enco ding [30]. This often results in far mo re IEEE TRANSACTIONS ON ME DICAL IMA GING 14 Fig. 15. Normalized s tandard devia tion ave raged over each tube for (a) T 1 and (b) T 2 . expensiv e co mputation , which mo tiv ates the in vestigation of advanced com p utational algor ithms (e.g., [58], [59]). For th e pro posed f ormulatio n , we m anually selected the parameters of th e weighting matrix and the constraint sets, which p rovided the goo d performance for th e app lication example. In g eneral, we sho uld choose these parameter s accordin g to application -speciﬁc requirem ents. For example , it may b e desirab le to co nstrain ﬂip an gles to b e small, an d apply multiple in version p ulses during acq uisition for cardia c MR ﬁnger printing experimen ts (e.g ., [ 60]). Another r e lated aspect is that the example im plementation of our pr oposed approach u tilizes on ly a f ew r e p resentative tissues to design MR ﬁngerp r inting experim ents. This is prac- tical, since the range of tissue pa rameter values is often k nown a p riori for a spec iﬁc imagin g applicatio n . Th is strategy gener- alizes well, as demonstrate d in the numer ical simulation s and real experiments. In the Supplemen tary Material, we include additional results to furth er con ﬁr m the good generalizatio n capability o f th is strategy . Alternatively , we can incorpo rate prior tissue para m eters into the experimen t design in oth er ways. For exam ple, we can perfor m a Bay e sian exper iment design by inco rporatin g a probability d istribution of tissues parameter values [32], and perf orm op timal experiment design with the Bayesian CRB [61]. O th er alternatives inclu de the min-max expe r iment d esign th a t min imizes the worst-case perfor mance [33] or average variance exper iment design that minimizes the average-case per forman ce acr oss an inter val of parameter values [62]. Regarding the num erical optim iza tion, it is worth reiter- ating that the pr oposed fo rmulation s in ( 16) an d (17) r esult in non-co nvex op timization problems, an d that the so lution algorithm is only guar anteed to have the local co nvergence. In this work , we initialized the algo rithm with the acquisition parameters from co nventional MR ﬁnger printing experiments, which demo nstrated the g ood pe r forman ce. In the Su pple- mentary Material, we ev aluated the impact of in itialization on the experimen t design. Simply pu t, we ﬁn d that different initializations fo r (16) can lead to different local m inima, although the cost fun ction values o f these local minim a are very close. Mor e interesting ly , we ﬁnd th at for (17), several very distinct initialization s pro duce the same optimized ac- quisition param eters. For a non-co n vex optimization pro blem, this behavior is rather remar kable, whic h is worth an in-depth theoretical study in the future. By solving (16) a n d (1 7), we observed that the o ptimized acquisition param eters appeared to be hig hly structured . In particular, the op tim ized repetition times a ppeared to be binary . Such an inter esting behavior is worth an in-de pth study . No te that as shown in [23], [63], cer tain CRB based exper iment design problem inv olving a dy namic system can be cast as an op timal control proble m . It is po ssible that our emp irical observations co uld b e explained by the princ ip les of op timal control the ory [24], [ 64]. 8 Although a control-th eoretic ch ar- acterization of (16) and (17) is beyond the scope o f this paper, it is an inter e sting topic to follow up on. I n pa r ticular, if we know in ad vance that (16) and (17) adm it structured solutions, this fact cou ld be leveraged to simp lify th e algorithms u sed to derive an optim al experiment design . From an o ptimal con trol viewpoint, the SQP algorith m f or solving (16) an d (17) ta rgets a necessary cond ition of an opti- mal contro l p roblem (i.e., the Po ntryagin ’ s minimum prin ciple [24]). As with othe r nonlin e a r p rogram ming algo rithms, the local minima obtaine d from the SQP algor ith m are gen erally not as goo d as those that would be obtain ed from d ynamic progr amming [23], [24], which, in con trast, deals with a sufﬁcient con dition of an op tim al con tr ol prob le m . However , the complexity of the SQP algorithm grows polyn omially with the nu mber of the d esign param eters, an d thus does not have the curse- o f-dimen sionality issue as does dy namic progr amming. In pr actice, the SQP solutions show substantial improvement over the conventional acqu isition, demonstra tin g their pra c tica l utility . W ith respect to the experimen tal validation, we demon - strated the key m erit of the pr oposed framework, i.e., variance reduction in par a m eter estimation , especially for T 2 estima- tion. For exam ple, we have demo nstrated th a t th e propo sed framework enables ab out a factor of two reductio n of the 8 T o our kno wledge, this work is the ﬁrst that reports Bang-Bang structure for the optimize d acquisi tion parameters in MR ﬁngerprinting. Although Maidens et al. formulat ed the experi ment design problem as an optimal cont rol problem and deri ved an approximat e dynamic programming algorit hm [23], the y only considered a highly-simpliﬁe d problem setup and did not report the Bang-Ba ng behavior . IEEE TRANSACTIONS ON ME DICAL IMA GING 15 Fig. 16. Reconstru cted paramete r maps and normalized standard devi ation maps for the in vi vo expe riments ( N = 400 ) using the acquisit ion paramet ers from the con vent ional scheme, Optimized-I, and Optimize d-II. (a) Reconstructe d T 1 maps. (b) Reconstruc ted T 2 maps. (c) Normalized standard de viation maps for T 1 . (d) Normalize d standard devia tion maps for T 2 . Fig. 17. Normaliz ed standard devi ation ave raged over the R OIs in the gray matter and white matter (marked in Fig. 16) for (a) T 1 and (b) T 2 . standard deviation f or T 2 . Such an improvement can be sub- stantial, which could potentially corr espond to a q uadru p ling of averaging time, or an upgr ade to a h igher ﬁeld scan ner . W e expect that it will have a po sitive impa c t o n the utility of MR ﬁngerpr in ting. While the CRB optim ization focuses on minimizin g th e variance with the assumption of unbiased estimation, the real par a m eter estima tio n error will also inclu de contributions from b ias. For numer ica l simulations, we demonstrated that Optimized-I I provides a similar level of bias a s th e c o n ven- tional sche me. Howev er, for real experim ents, p articularly in viv o exper iments, assessing the bias can be more inv olved. First, note that there is n o a n un derlying g round tru th f or in viv o experiments [65]. Sec o nd, various model mismatches can often comp licate the bias ev aluatio n a nd compar ison. For example, the par tial volume effect o f ten yie ld different mag- netization ev olution s unde r different signa l excitations, which can result in signiﬁcant difﬁculty when co m paring different acquisition schemes [ 33], [65]. A thorou gh assessment o f the bias for parameter estimation still re mains an open problem for MR ﬁngerprintin g, and for quantitative M R imaging in general. This pa p er sy stem atically investigated one speciﬁc im- plementation of our p ropo sed framework, but many o ther implementatio ns are p ossible and some of these altern ativ es may b e preferab le d e pending on th e imag ing co ntext. For example, we used the weighted A-op tim ality as a cr iterion to design MR ﬁngerpr inting experimen ts an d demonstrated its good perf ormance . Th ere ar e alter native de sig n criteria (e.g., the D-optima lity) that could b e inco r porated into the propo sed fra mew ork , and the explo ration of these criteria m ay be b eneﬁcial in cer tain app lications. As ano th er example, we can also include other acq uisition para meters (e.g., RF p ulse phases, ech o times, e tc) into the optimization . W ith a larger search spac e o f acquisition p arameters, we co u ld ach iev e b etter perfor mance, altho ugh the resultin g optimization pro blem will be computatio nally more expensive. Add itionally , w e cou ld optimize MR ﬁng erprintin g expe riments to in clude other tis- sue param eters, such a s appar ent diffusion coefﬁcients [66], magnetizatio n tran sf e r [67], [68], etc. Lastly , altho ugh the implementatio n in this p aper fo cused on op timizing the CRB for a ﬁxed n umber of TRs, there are many other design con - siderations (e.g., acoustic characteristics of th e pulse sequ ence [69], SAR, total acq uisition time , etc), which are importan t and may be worth in cluding within the optimization framework. Explorin g these possibilities may fu rther improve the practical signiﬁcance of our app r oach. IEEE TRANSACTIONS ON ME DICAL IMA GING 16 V . C O N C L U S I O N In this work, we presented a novel estimation- th eoretic framework to o ptimize acqu isition param eters of MR ﬁnge r- printing exper iments. W e formu lated th e optimal experiment design problem that maximizes the SNR e fﬁciency of a n MR ﬁngerprin ting experiment, while incorp orating addition a l constraints that are a d vantageous to th e reco nstruction p rocess. The optimize d e xpe riments enables substantially imp roved accuracy f or T 2 maps, wh ile providin g similar o r slightly better acc uracy for T 1 maps. Rema r kably , we fo u nd that the optimized acquisition parameters a p pear to be highly structured, rather than r a ndom/p seudo-r a n domly varying as used in the conventional MR ﬁn gerprin ting exper im ents. V I . A C K N O W L E D G E M E N T The au thors would like to thank the anonym ous r evie wers for their constructive c o mments, which help improve this paper . B. Zhao would like to than k Dr . P . Polo for the help with cr e ating the p h ysical phan tom, and Drs. D. Park and Y . Chang f or the discussions on the corr ection of a cquisition imperfectio ns. V I I . A P P E N D I X In this append ix, we der iv e the Jaco bian m atrix J n ( θ ) in (12) for the IR-FISP sequen c e. This provides an exam ple to illustrate the proced ure describ ed in ( 1 3) and (14). Recall that in th e IR-FISP sequ ence, θ = [ T 1 , T 2 , M 0 ] T and J n ( θ ) = ∂ m [ n ] ∂ θ =  ∂ ∂ T 1 m [ n ] ∂ ∂ T 2 m [ n ] ∂ ∂ M 0 m [ n ]  . For clarity , we ﬁrst summarize the state-spa c e mod el for the IR-FISP sequ ence as fo llows: M r [ n ] = G ( β r ) R ( T 1 , T 2 , T R n ) Q ( α n , φ n ) M r [ n − 1] + M 0 N v b ( T 1 , T R n ) , (21) m [ n ] = X r PR ( T 1 , T 2 , T E n ) Q ( α n , φ n ) M r [ n − 1] , (22) for n = 1 , · · · , N . Next, we calculate the der ivati ves o f m [ n ] with re sp ect to T 1 , T 2 , an d M 0 based on (13) a n d (14). A. Deriva tive of m [ n ] with r espect to T 1 In voking th e derivati ve with respe c t to T 1 on bo th sides o f (22), we have ∂ m [ n ] ∂ T 1 = X r  P ∂ R ( T 1 , T 2 , T E n ) ∂ T 1 Q ( α n , φ n ) M r [ n − 1] + PR ( T 1 , T 2 , T E n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ T 1  , (23) where ∂ R ( T 1 , T 2 , T E n ) ∂ T 1 = T E n T 2 1 exp ( − T E n T 1 )   0 0 0 0 0 0 0 0 1   . Noting tha t P ∂ R ( T 1 , T 2 , T E n ) ∂ T 1 = 0 , (23) can b e simp liﬁed as ∂ m [ n ] ∂ T 1 = X r PR ( T 1 , T 2 , T E n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ T 1 . (24) W e the n take the derivati ve with respe ct to T 1 on b oth sides of (21), i.e., ∂ M r [ n ] ∂ T 1 = G ( β r ) ∂ R ( T 1 , T 2 , T R n ) ∂ T 1 Q ( α n , φ n ) M r [ n − 1] + G ( β r ) R ( T 1 , T 2 , T R n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ T 1 + M 0 N v ∂ b ( T 1 , T R n ) ∂ T 1 , (25) where ∂ R ( T 1 , T 2 , T R n ) ∂ T 1 = T R n T 2 1 exp ( − T R n T 1 )   0 0 0 0 0 0 0 0 1   , and ∂ b ( T 1 , T R n ) ∂ T 1 = − T R n T 2 1 exp ( − T R n T 1 )   0 0 1   . Here we can calculate ∂ m [ n ] /∂ T 1 by iterating (24) and (25) with th e initial con ditions M r [0] = M 0 N v [0 0 1] T and ∂ M r [0] /∂ T 1 = [0 0 0] T . B. Deriva tive of m [ n ] with r espect to T 2 In voking the derivati ve respect to T 2 on both sides of (22), we have ∂ m [ n ] ∂ T 2 = X r  P ∂ R ( T 1 , T 2 , T E n ) ∂ T 2 Q ( α n , φ n ) M r [ n − 1] + PR ( T 1 , T 2 , T E n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ T 2  , (26) where ∂ R ( T 1 , T 2 , T E n ) ∂ T 2 = T E n T 2 2 exp ( − T E n T 2 )   1 0 0 0 1 0 0 0 0   . W e the n take the derivati ve with respe ct to T 2 on b oth sides of (21), i.e., ∂ M r [ n ] ∂ T 2 = G ( β r ) ∂ R ( T 1 , T 2 , T R n ) ∂ T 2 Q ( α n , φ n ) M r [ n − 1] + G ( β r ) R ( T 1 , T 2 , T R n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ T 2 , (27) where ∂ R ( T 1 , T 2 , T R n ) ∂ T 2 = T R n T 2 2 exp ( − T R n T 2 )   1 0 0 0 1 0 0 0 0   . Here we can calculate ∂ m [ n ] /∂ T 2 by iterating (26) and (27) with th e initial con ditions M r [0] = M 0 N v [0 0 1] T and ∂ M r [0] /∂ T 2 = [0 0 0] T . IEEE TRANSACTIONS ON ME DICAL IMA GING 17 C. Derivative of m [ n ] with respect to M 0 In voking the d e riv ativ e with respect to M 0 on both sides of (22), we have ∂ m [ n ] ∂ M 0 = X r PR ( T 1 , T 2 , T E n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ M 0 . (28) W e then take th e derivati ve with respect to M 0 on bo th sid e s of (21), i.e., ∂ M [ n ] ∂ M 0 = G ( β r ) R ( T 1 , T 2 , T R n ) Q ( α n , φ n ) ∂ M r [ n − 1] ∂ M 0 + 1 N v b ( T 1 , T R n ) . (29) Here we can calculate ∂ m [ n ] / ∂ M 0 by iterating (2 8) and (29) with the initial co ndition ∂ M r [0] /∂ M 0 = 1 N v [0 0 1] T . R E F E R E N C E S [1] D. Ma, V . Gulani, N. Seiberlich, K. L iu, J. L. Sunshine, J. L. Duerk, and M. A. G riswold, “Magne tic resonance ﬁngerprinting, ” Nature , vol . 495, pp. 187–192, 2013. [2] K. Schefﬂer and J. Hennig, “T1 quantiﬁca tion with in version recov ery True FISP, ” Magn. Reson. Med. , vol. 45, pp. 720–723, 2001. [3] S. C. Deoni, B. K. Rutt, and T . M. Peters, “Rapid combined T 1 and T2 mapping using gradient recall ed acquisition in the steady s tate, ” Magn. Reson. Med. , vol. 49, pp. 515–526, 2003. [4] P . Schmitt, M. A. Grisw old, P . M. Jakob , M. Kotas, V . Gulani , M. Flent je, and A. Haase, “In version recov ery TrueFISP: Quantiﬁcation of T1, T2, and spin density, ” Magn . Reson. Med. , vol. 51, pp. 661–667, 2004. [5] J. W arntjes, O. Dahlqvist, and P . Lundberg, “Nov el method for rapid, simultane ous T1, T2*, and proton density quantiﬁcat ion, ” Magn. Reson. Med. , vol. 57, pp. 528–537, 2007. [6] M. Davies, G. Puy , P . V anderg heynst, and Y . Wia ux, “ A compressed sensing frame work for magnetic resonance ﬁngerprinting , ” SIAM J. Imagin g Science s , vol. 7, pp. 2623–2656, 2014. [7] B. Zhao, F . Lam, B. Bilgic, H. Y e, and K. Setsompop, “Maximum lik elihood reconstruction for magnetic resonance ﬁngerprint ing, ” in Pr oc. IEEE Int. Symp. Biomed. Imaging , 2015, pp. 905–909. [8] B. Z hao, K. S etsompop, H. Y e, S. Caule y , and L. L. W ald, “Maximum lik elihood reconstructi on for m agneti c resonance ﬁngerprinting, ” IEEE T rans. Med. Imaging , vol. 35, pp. 1812–1823, 2016. [9] E . Y . Pierre, D. Ma, Y . Chen, C. Badve, and M. A. Griswold, “Multiscale reconstru ction for MR ﬁngerprint ing, ” Magn. Reson. Med. , vol. 75, pp. 2481–2492, 2016. [10] X. Cao, C. L iao, Z. W ang, Y . Chen, H. Y e, H. He, and J. Zhong, “Rob ust sliding-wi ndow reconstru ction for accel erating the acquisit ion of MR ﬁngerprint ing, ” Magn . Reson. Med. , vol. 78, pp. 1579–1588, 2017. [11] B. Zhao, K. Setsompop, E. Adalsteinsson , B. Gagoski, H. Y e, D. Ma, Y . Jiang, P . E . Grant, M. A. Griswold , and L. L. W ald, “Improv ed magnetic resonanc e ﬁngerprinting reconstruct ion with lo w-rank and subspace modeling, ” Magn. Reson. Med. , vol. 79, pp. 933–942, 2018. [12] A. Pazman, F oundat ions of Optimum Experimental Design . Nethe r- lands: Springer , 1986. [13] F . Pukelshe im, Optimal Design of Experiments . New Y ork, NY : John W iley & Sons, 1993. [14] S. M. Kay , Fundamentals of Statistical Signal Proc essing: Estimation Theory . Upper Saddle Riv er , NJ: Printice Hall, 1993, vol. I. [15] Y . Jiang, D. Ma, N. Seiberl ich, V . Gulani, and M. A. Griswold, “MR ﬁngerprint ing using fast imaging with steady state precession (FISP) with spiral readout , ” Ma gn. R eson. Med. , vol . 74, pp. 1621–1631, 2015. [16] J. Assl ¨ ander , S. J. Glaser , and J . Hennig, “Pseudo s teady-state free precessio n for MR-ﬁngerprint ing, ” Magn . Reson. Med. , vol. 77, pp. 1151–1161, 2017. [17] Y . Jiang, D. Ma, R. Jerecic, J. Duerk, N. Seiberlich, V . Gulani, and M. A. Griswold, “MR ﬁngerpr inting using the quick echo splitt ing NMR imaging techniq ue, ” Magn. Reson. Med. , vol. 77, pp. 979–988, 2017. [18] B. Rieger , F . Zimmer , J. Zapp, S. W eing ¨ artner , and L. R. Schad, “Magnet ic resonance ﬁngerprinti ng using echo-plana r imaging: Joint quantiﬁc ation of T1 and T 2 * relaxa tion times, ” Magn. Reson. Med. , vol. 77, pp. 979–988, 2017. [19] B. Z hao, J. P . Haldar , K. Setsompop, and L. L. W ald, “T oward s optimize d expe riment design for magneti c resonanc e ﬁngerprint ing, ” in Pr oc. Int. Symp. Magn. Reson. Med. , 2016, p. 2835. [20] ——, “Optimal experi ment design for magnetic resonance ﬁngerprint- ing, ” in Proc. IEEE Eng. Med. Bio. Conf. , 2016, pp. 453–456. [21] B. Zhao, J. P . Haldar , C. Liao, D. Ma, M. A. Griswol d, K . Setsompop, and L. L . W ald, “Optimal experiment design for magnetic resonanc e ﬁngerprint ing: New insights and further improvement s, ” in Pr oc. Int. Symp. Magn. Reson. Med. , 2018, p. 674. [22] J. Ass l ¨ ander , D. Sodickson, R. Lattanzi, and M. Cloos, “Relaxa tion in polar coordinates: Analysis and optimizati on of MR-ﬁngerprintin g, ” in Pr oc. Int. Symp. Magn. Reson. Med. , 2017, p. 127. [23] J. Maidens, A. Packa rd, and M. Arcak, “P arallel dynamic programming for optimal experiment design in nonlinear systms, ” in P r oc. 55th IEEE Conf . on Decision and Contr ol , 2016, pp. 2894–2899. [24] D. P . Bertsekas, Dynamic Pr ogrammin g and Optimal Contr ol , 4th ed. Belmont, MA: Athena Scientif c, 2017, vol . I. [25] O. Cohen and M. S. Rosen, “ Algorithm comparison for schedule optimiza tion in MR ﬁngerprinti ng, ” Magn. Reson. Imaging , vol. 41, pp. 15–21, 2017. [26] J. A. Jones, P . Hodgkinson, A. L. Barker , and P . J. Hore, “Optimal sampling strate gies for the measurement of spin-spin relaxat ion times, ” J . Magn . Reson. Ser B , vol. 113, pp. 25–34, 1996. [27] Y . Z hang, H. N. Y eung, M. O’Donnel l, and P . L. Carson, “Deter mination of sample time for T 1 m easurement , ” J. Magn. Reson. Imaging , vol. 8, pp. 675–681, 1998. [28] J. P . Haldar , D. Herna ndo, and Z. -P . Liang, “Super -resolution reconstruc- tion of MR image sequenc es with contrast m odeling , ” in Pr oc. IEEE Int. Symp. Biomed. Imaging , 2009, pp. 266–269. [29] A. Funai and J. A. Fessler , “Cram ´ er Rao bound analysis of joint B1/T1 mapping methods in MRI, ” in Proc. IEEE Int. Symp. Biomed. Imaging , 2010, pp. 712–715. [30] B. Zhao, F . Lam, and Z.-P . Liang, “Model-base d MR parameter m ap- ping with sparsity constraints: Parame ter estimation and performance bounds, ” IEE E T rans. Med. Imaging , vol. 33, pp. 1832–1844, 2014. [31] M. Akc ¸ akaya, S. W eing ¨ artner , S. Roujol, and R. Nezafat, “On the select ion of sampling points for myocardi al T1 mapping, ” Magn. Reson. Med. , vol. 73, pp. 1741–175 3, 2015. [32] C. M. Lewis, S. A. Hurley , M. E. Meyera nd, and C. G. Ko ay , “Data- dri ven optimize d ﬂip angle selectio n for T1 estimation from spoiled gradien t echo acquisitions, ” Magn. Reson. Med. , vol. 76, pp. 792–802, 2016. [33] G. Natara j, J. F . Nielsen, and J. A. Fessler , “Optimizing MR scan design for model-based T1, T2 estimati on from steady-st ate sequence s, ” IEEE T rans. Med. Imaging , vol. 36, pp. 467–477, 2017. [34] R. P . A. T eixe ira, S. J. Malik, and J. V . Hajnal, “Joint system relax- ometry (JSR) and Cram ´ er-Rao lowe r bound optimization of sequence paramete rs: A frame work for enhanced precision of DESPOT T1 and T2 estimatio n, ” Magn. Reson. Med. , 2017, in press. [35] D. Luenber ger , Intr oduction to Dynamic Systems: Theory , Models, and Applicat ions . Ne w Y ork, NY : John W iley & Sons, 1979. [36] F . Bloch , “Nuclear induction, ” P hys. Rev . , vol. 70, pp. 460–474, 1946. [37] Z. -P . Liang and P . C. Lauterb ur , Principles of Magnet ic Resonance Imagin g: A Signal P r ocessing P erspectiv e . N ew Y ork, NY : IEE E Press/W ile y , 1999. [38] P . Shkarin and R. G . Spencer , “Time domain simulati on of Fourier imaging by sum m ation of isochromats, ” Int J . Imag. Syst. T ech . , vol. 8, pp. 419–426, 1997. [39] S. J. Malik, A. Sbrizzi, H. Hoogduin, and J. V . Hajnal , “Equi va lence of E PG and isochromat -based simulati on of MR signals, ” in P r oc. Int. Symp. Magn. Reson. Med. , 2016, p. 3196. [40] V . V . Fedorov , Theory of Optimal E xperiment s . Ne w Y ork, NY : Academic Press, 1972. [41] J. Xie, D. Gallichan, R. N. Gunn, and P . Jezzard, “Optimal design of pulsed arteria l spin labeling MRI expe riments, ” Magn. Reson. Med. , vol. 59, pp. 826–834, 2008. [42] D. H. Poot, A. J. den Dekker , E. Achten, M. V erhoye , and J. Sijbers, “Optimal expe rimental design for diffusio n kurtosis imaging, ” IE EE T rans. Med. Imaging , vol. 29, pp. 819–829, 2010. [43] B. Lampinen, F . Szczepanki ewic z, D. van W esten, E. Englund, P . C Sundgren, J. L ¨ att, F . St ˚ ahlberg, and M. Nilsson, “Optimal ex- periment al design for ﬁlter excha nge imaging: Apparent exch ange rate measurement s in the healthy brain and in intracrania l tumors, ” Magn. Reson. Med. , vol. 77, pp. 1104–1114, 2017. [44] B. Zhao, W . Lu, T . K. Hitchens, F . Lam, C. Ho, and Z .-P . Liang, “ Accel- erated MR parameter m apping with lo w-rank and sparsity constra ints, ” Magn. Reson. Med. , vol. 74, pp. 489–498, 2015. IEEE TRANSACTIONS ON ME DICAL IMA GING 18 [45] B. Zhao, “Model-based iterati ve reconstruction for magnetic resonance ﬁngerprint ing, ” in P roc. IEEE Int. Conf. Image Proc ess. , 2015, pp. 3392– 3396. [46] J. Nocedal and S. J. W right, Numerical Optimizatio n . Ne w Y ork, NY : Springer , 2006. [47] J. C. Spall, Intr oduction to Stoc hastic Searc h and Optimization: Esti- mation, Simulation, and Contr ol . Ne w Y ork, NY : John W iley & Sons, 2003. [48] D. L. Collins, A. P . Zijdenbos, V . Kollo kian, J. G. Sled, N. J. Kabani, C. J. Holmes, and A. C. Evans, “Design and construct ion of a realistic digita l brain phantom, ” IEE E T rans. Med. Imaging , vo l. 17, pp. 463–468, 1998. [49] J. A. Fessler and B. P . Sutton, “Nonuniform fast Fourier transforms using min-max interpolatio n, ” IEEE Tr ans. Signal Pr ocess. , vol. 51, pp. 560–574, 2003. [50] Y . Gao, Y . Chen, D. Ma, Y . Jiang, K. A. Herrmann, J. A. V incent , K. M. Dell, M. L . Drumm, S . M. Brady-Ka lnay , M. A. Griswold, C. A. Flask, and L. Lu, “Preclini cal MR ﬁngerprint ing (MRF) at 7 T: Effecti ve quantit ati ve imaging for rodent disease models, ” NMR Biomed. , vol. 28, pp. 384–394, 2015. [51] J. P . Haldar , J. Anderson, and S. W . Sun, “Maximum likel ihood estimati on of T 1 relaxa tion parameters using V ARPR O, ” in Proc. Int. Symp. Magn. Reson. Med. , 2007, p. 41. [52] K. Hattori, Y . Ikemoto, W . T akao, S. Ohno, T . Harimoto, S. Kanaza wa, M. Oita, K. Shibuya, M. Kuroda, and H. Kato, “De velopmen t of MRI phantom equi va lent to human tissues for 3.0-T MRI, ” Med. Phys. , vol. 40, p. 032303, 2013. [53] H. T an and C. H. Meyer , “Estimation of k-space trajectorie s in spiral MRI, ” Magn. Reson. Med. , vol. 61, pp. 1396–1404, 2009. [54] D. Ma, S. Coppo, Y . Chen, D. F . McGivne y , Y . J iang, S. Pahwa, V . Gulani, and M. A. Griswold, “Slice proﬁle and B1 corrections in 2D magneti c resonanc e ﬁngerpri nting, ” Magn. Reson. Med. , vol. 78, pp. 1781–1789, 2017. [55] M. W eigel, “Extended phase graphs: Dephasing , RF pulses, and echoes - pure and simple, ” J. Magn. Reson. Imaging , vol. 41, pp. 266–295, 2015. [56] X. Dan, M. Jacob, and Z.-P . Liang, “Optimal sampling of k-space on Cartesi an grids for parall el imaging, ” in Proc . Int. Symp. Magn . Reson. Med. , 2005, p. 2450. [57] C. Ma, D. Xu, K. F . King, and Z.-P . Liang, “Joint design of spoke trajec tories and RF pulses for paral lel exc itation, ” Magn. Reson. Med. , vol. 65, pp. 973–985, 2011. [58] A. Hero and J. A. Fessler , “ A recursiv e algori thm for computing Cramer- Rao-type bounds on estimato r cov ariance, ” IEEE T rans. Inf. Theory , vol. 40, pp. 1205–1210 , 1994. [59] P . Tune , “Computing constraine d Cram ´ er-Rao bounds, ” IEEE T rans. Signal Pr ocess. , vol. 60, pp. 5543–5548, 2012. [60] J. I. Hamilton, Y . Jiang, Y . Chen, D. Ma, W .-C. Lo, M. A. Griswold, and N. Seibe rlich, “MR ﬁngerprin ting for rapid quantiﬁcati on of myocardial T1, T2, and proton spin density, ” Magn. Reson. Med. , vol. 77, pp. 1446– 1458, 2017. [61] H. L. V an Tree s and K. L. Bell, Bayesian Bounds for P arameter Estimation and Nonlin ear F iltering/ T rac king . Ne w Y ork, N Y : Wi ley- IEEE Press, 2007. [62] J. P . Haldar , Q. Gao, X. J. Zhao, and Z. -P . L iang, “Optimized measure- ment of anomalous diff usion, ” in P roc. Int. Symp. Magn. R eson. Med. , 2009, p. 3570. [63] R. T . N. Chen, “Input design for aircraft parameter identiﬁcat ion: Using time-opt imal control formulation, ” in Methods for Airc raft State and P arameter Identiﬁcat ion, Advisory Group for Aerospa ce Resear ch and Deve lopment (AGARD), Confer ence P r oceeding s , no. 172, 1975. [64] N. P . Osmolovskii and H. Maurer , Application s to Re gular and Bang- Bang Contr ol: Second-Or der Necessary and Sufﬁcien t Optimali ty Con- ditions in Calculus of V ariations and Optimal Contr ol . Philade lphia, P A: SIAM, 2012. [65] N. Stiko v , M. Boudreau, I. R. Le vesque, C. L . T ardif, J. K. Barral, and G. B. Pike, “On the accurac y of T1 mapping: Searchin g for common ground, ” Magn. Reson. Med. , vol. 73, pp. 514–522, 2015. [66] Y . Jiang, D . Ma, K. Wright, N. Seiberlich, V . Gulani , and M. A. Gris- wold, “Simultaneous T1, T2, diffusion and proton density quantiﬁc ation with MR ﬁngerprinti ng, ” in P r oc. Int. Symp. Ma gn. Reson. Med. , 2014, p. 28. [67] C. Y . W ang, Y . Liu, S. Huang, M. A. Griswold, N. Seiberlic h, and X. Y u, “ 31 P magneti c resonance ﬁ ngerprinting for rapid quantiﬁca tion of creatine kinase react ion rate in viv o, ” NMR in B iomed. , vol. 30, p. e3786, 2017. [68] T . Hilbert, T . Kob er , T . Zhao, T . K. Block, Z. Y u, J.-P . Thiran, G. Kruege r , D. K. Sodickson, and M. Cloos, “Mitigati ng the ef fect of magnetiz ation transfer in magnetic resonance ﬁngerprinting, ” in Proc. Int. Symp. Magn. Reson. Med. , 2017, p. 74. [69] D. Ma, E. Y . Pierre, Y . Jian g, M. D. Schluc hter , K. Setsompop, V . Gulani, and M. A. Griswold, “Music-based magnetic resonanc e ﬁngerprinting to improve patient comfort during MRI examin ations, ” Magn. Reson. Med. , vol. 75, pp. 2303–231 4, 2016.

Optimal Experiment Design for Magnetic Resonance Fingerprinting: Cramer-Rao Bound Meets Spin Dynamics

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment