Optimal Design with Automatic Differentiation: Exploring Unbiased Steady-State Relaxometry

1. Optimal Design with Automatic Differentiation:Exploring Unbiased Steady-State Relaxometry Jan 13, 2014 Jason Su

2. Outline • Concept and theory • Practical usage with automatic differentiation (AD) • Exploring relaxometry with SPGR and bSSFP

3. Motivation • I started with the goal of trying to optimize mcDESPOT and the estimate of MWF • Lankford et al. criticized the theoretical accuracy of mcDESPOT based on CRLB analysis • Interestingly however, this criticism inspired the approach to finding the solution • Since then, I’ve discovered that this project can have greater possibilities • The following provides a framework that is suitable for the analysis, comparison, and optimization of a broad range of parameter mapping/estimation methods • T1, T21, B1+, B0, MT, diffusion coefficient2, etc. [1] Jones et al. JMR 1996 Oct;113(1):25-34. [2] Brihuega-Moreno et al. MRM 2003 Nov;50(5):1069-76.

4. Idea • The Cramer-Rao Lower Bound (CRLB) provides a “best case” limit on the variance of an estimator of parameter • In quantitative MR, the estimator is defined by the applied pulse sequence and protocol details • Typically this is used as an analysis tool, but what if we instead use it as a cost function? • We want to solve the problem: • Find the protocol that gives the lowest variance estimate of θ • This is hard to solve for non-linear equations, relax it to: • Find the protocol that gives the best lower bound on the variance of the estimate of θ

5. CRLB: Fisher Information Matrix • Typically calculated for a given tissue, θ • Interpretation • J captures the sensitivity of the signal equation to changes in a parameter • Its “invertibility” or conditioning is howseparable parameters are from each other,i.e. the specificity of the measurement

6. CRLB: How does it work? • A common formulation • Unbiased estimator • A signal equation with normally distributed noise • Measurement noise is independent

7. CRLB: Computing the Jacobian

8. Automatic Differentiation • An algorithm for generating computer programs that calculate derivatives of other functions based on their source code/graph • Uses repeated application of the chain rule • Fast and accurate to machine precision • Can differentiate extremely complicated functions • Many packages exist for Matlab, C, and Python • Typically you will write the function that computes your signal equation • Then your AD package provides the derivatives essentially for free, i.e. no analysis or coding effort • YMMV depending on the language or package

9. What does this buy us?

10. Application • To vet this approach and my code, which is ultimately intended for open-source use • Went back to analyze a number of classic designs to see if I could reproduce them • Diffusion coefficient, VFA/DESPOT1, DESPOT2 • In all cases, the answer was yes • The DESPOT2-cast was particularly interesting

11. Steady-State Relaxometry with bSSFP and SPGR (aka DESPOT2) • DESPOT2 seeks to estimate T1 and T2 • First do VFA/DESPOT1 with a pair of SPGR images at ~0.71 signal level to get a T1 map • Then acquire a pair of π-phase-cycled bSSFP images at ~0.71 signal level and back out T2, knowing T1 from step 1 • Constraints: • Phase-cycling is π • Acquisition time is equal between the angles in a pair, only optimize the relative time between the SPGR and bSSFP parts of the experiment

12. Methods • Find α’s and acquisition time fraction • Minimize the sum of the coefficient of variations (CoV) of T1 and T2: σT1/T1+σT2/T2 • Where σ was provided by the CRLB which uses the Jacobians of the SPGR and bSSFP signal equations • Exhaustively search for the # SPGR and bSSFP images • Fix the total acquisition time so we find the best protocol per unit time • Optimize for representative T1 and T2 values of WM and GM at 3T • WM (T1=1100ms, T2=60ms) and GM (T1=1645ms, T2=85ms)

13. Solver • Used a general solver for this problem • This is the weakest part of the framework as the solver sometimes fails to converge or encounters inversion errors • May be improved with insights from optimal design theory, e.g. using determinant of F instead of actual CRLB

14. Results • To within 5 decimal places, framework says to acquire the pairs of flip angles at 1/√2 of the max signal for SPGR and π-phase-cycled bSSFP • Compared to the original literature estimate of 0.71 (Deoni 2003) • I later found that 1/√2 was shown to be the analytically correct solution (Schabel and Morell 2009) • ~75% of time should be spent on SPGRs and T1 map • Optimal DESPOT2 protocols achieve a sum of T1 and T2 CoVs of 45.7 and 53.0 for WM and GM

15. Unrestricted Problem • Can we do better? • What if we remove the constraints from the DESPOT2 acquisition? • Allow any combination of SPGRs and bSSFPs to be used, forget about the 2-stage estimate of T1 then T2 and reconstruct with a joint nonlinear fit • Allow any phase-cycle to be used • Allow the acquisition time for any image to be freely adjustable

16. PCVFA Solution • The 0 and π phase cycles should be collected, each with a pair of flip angles that attain 1/√2 of the max signal for that phase-cycle • Equal time fraction should be devoted to each of the 4 images • Achieves a sum of CoVs of 20.8 and 21.6 for WM and GM • More than 2.1x improvement over DESPOT2

17. Performance of Optimal GM Protocol over a Range of Tissues DESPOT2 PCVFA

18. Experimental Validation • In progress with silicone oil phantom at 3T to mitigate B1 inhomogeneity • Harder to acquire DESPOT2 optimal sequence due to time difference between SPGR and bSSFP parts • Fractional angles by adjusting ia_rf1, <1 deg. Needed for PCVFA • Had difficulties doing linear and nonlinear fit for DESPOT2 and PCVFA with acquired data • Maybe made an error in the acquisition, should try again

19. One step further • What if we used complex images for the reconstruction? • 21.6 -> 10.81 • Another 2x gain?

20. Complex PCVFA vs. Magn. PCVFA

21. Summary • This framework gives the ability to study complex pulse sequences and their efficacy analytically • A way to formally compare and explore potential mapping methods • Gets closer to answering the question of, what is the best way of mapping X? • But may not include considerations like immunity to B1+ inhomogeneity or ease of fitting

22. Next Steps • Robust protocol design • Have been looking at B0 robustness, where cost function evaluates CRLB at many B0 values and try to optimize the worst case estimate of θ among them • B1+ robustness • Broad comparison of efficiency of mapping methods • T1 or B1+ mapping is probably a good candidate • mcDESPOT and MWF • MPnRAGE, what should n be? • Analyzing bloch simulations, MRF?