Optimal control theory for applications in magnetic resonance imaging

We apply innovative mathematical tools coming from optimal control theory to improve theoretical and experimental techniques in Magnetic Resonance Imaging (MRI). This approach allows us to explore and to experimentally reach the physical limits of the corresponding spin dynamics in the presence of typical experimental imperfections and limitations. We study in this paper two important goals, namely the optimization of image contrast and the maximization of the signal to noise per unit time. We anticipate that the proposed techniques will find practical applications in medical imaging in a near future to help the medical diagnosis.


Introduction
Optimality with respect to a given criterion is vital in many applications, but it presents a complexity that requires a lot of ingenuity to provide a solution. In this context, optimal control tackles the question of bringing a dynamical system from one state to another with minimum expenditure of time and resources [1]. Optimal control theory was born in its modern version with the Pontryagin Maximum Principle (PMP) in the late 1950's. Its development was originally inspired by problems of space dynamics, but it is now a key tool to study a large spectrum of applications extending from robotics to economics and biology [2,3,4,5,6,7,8,9,10]. Optimal control problems can be solved by two different types of approaches, geometric [2,3,5] and numerical methods [8,11,12] for dynamical systems of low and high dimension, respectively. The geometric techniques lead to a complete mathematical and qualitative understanding of the control problem, from which we can deduce the structure of the optimal solution and the physical limits of a dynamical process, such as the minimum time to achieve a given task [2,3,5,13,14,15,16,17,18]. The advantages of the numerical approach are complementary to the ones of the geometric method. The relative simplicity of the application of the numerical algorithms makes it direction both from the geometric and numerical approaches. These examples will allow us to describe and discuss the efficiency and the limitations of this method in order to solve applied and concrete issues in MRI.
The paper is organized as follows. Section 2 deals with the optimization of the signal to noise per unit time (SNR) by geometric control techniques. The unbounded case is investigated, showing that the Ernst angle solution is the optimal solution of the control problem. This result is an essential prerequisite for the extension of this analysis to realistic situations, which are numerically much more intricate. Section 3 focuses on the optimization of the contrast in in vivo MRI. We numerically maximize the contrast of the brain of a rat. The experimental results show the efficiency of the applied magnetic fields. Conclusion and prospective views are given in Sec. 4. 2 Optimal control of the signal to noise ratio per unit time This paragraph deals with an application of geometric optimal control theory to MRI. Recently, this method has been successfully applied to a question of fundamental and practical interest in MRI, the maximization of the achievable signal-to-noise ratio per unit time (SNR) of a spin 1/2 particle [24,25,26]. The SNR is practically enhanced in spin systems by using a multitude of identical cycles. In this periodic regime, the SNR increases as the square root of the number of scans. Each elementary block is composed of a detection time and of a control period where the spin is subjected to a radio-frequency magnetic pulse, this latter being mandatory to guarantee the periodic character of the overall process. A first solution to this problem was established in the sixties by R. Ernst and his co-workers [48]. In this proposition, the control law is made of a δ-pulse, featured by a specific rotation angle, and known as the Ernst angle solution. This pulse sequence is routinely used in magnetic resonance spectroscopic and medical applications. Note also that optimal control theory was not used in Ref. [48]. We have revisited recently this question by applying the powerful tools of geometric control theory. In [49], we have shown in the case of unbounded controls that the Ernst angle solution is the optimal solution of this control problem. This analysis was generalized in [50] to spin dynamics in the presence of radiation damping effects and crusher gradients. We present in this section a brief review of these results.

The model system
We describe in this paragraph the model system used in the theoretical analysis. We consider an ensemble of spin 1/2 particles subjected to a radio-frequency magnetic field, which is assumed to be homogeneous across the sample. The system is described by a magnetization vector M of coordinates (M x , M y , M z ), whose dynamics is governed, in a given rotating frame, by the standard Bloch equations [24,25,26]: where T 1 and T 2 are respectively the longitudinal and transverse relaxation constants. M 0 is the thermal equilibrium of the magnetization and (ω x , ω y ) are the two control amplitudes of the radio-frequency magnetic field. Normalizing the time with respect to the detection time T d (see below for a description of this parameter), τ = t/T d , and the amplitude of the magnetization vector with respect to M 0 , (x, y, z) = (M x /M 0 , M y /M 0 , M z /M 0 ), we get: where γ = 2πT d /T 1 , Γ = 2πT d /T 2 and u x,y = T d ω x,y . In the relevant physical case where 0 ≤ T 2 ≤ 2T 1 , it can be shown that x 2 + y 2 + z 2 ≤ 1 at any time, which defines the Bloch ball [51]. In the new coordinates, the dynamics only depends on two parameters γ and Γ and the detection time is fixed to 1. The system admits a symmetry of revolution around the z-axis which allows us to consider by setting u y = 0 only a planar projection onto the (y, z)-plane [52,53]. The final equations used in this paper are then given by: where u stands for u x . In this ideal dynamical system, we introduce a simple scenario to describe the optimization of the SNR per unit time (see the schematic description in  problem is defined through the introduction of the following cost functional R: y m being the strength of the signal at the beginning of the detection period (only the transverse magnetization can be measured). A white noise is assumed, which explains the √ N factor in the denominator of R. Simple algebra leads to: Since the total time T and the detection time T d are fixed,we can introduce a normalized SNR Q with respect to the measurement time T d : which only depends on the position of the M point. The factor of quality Q is the figure of merit to maximize.

How to move in the Bloch ball
Before entering into the details of the optimization of the SNR, the first question to solve is the optimal control of the magnetization vector in the Bloch ball.
Here, we only summarize the main results, which have been recently established in a series of papers [51,53,54,55,56]. We recall that the dynamics of the system is governed by Eq. (3), with no constraint on the control field, u(t). At this stage, it is illuminating to introduce the polar coordinates (R, θ) such that y = R sin θ and z = R cos θ. Equations (3) can then be written as follows: We immediately note that the angular speedθ can be directly control by u, while the radial one,Ṙ, does not explicitly depend on the field. The radial coordinate can only be controlled in a two-step process by a judicious choice of θ. In the unbounded case, we have therefore a complete control over the angular degree of freedom, but we still have to understand how to manipulate the system along the radial direction. For that purpose, we have plotted in Fig. 2 the evolution of the speedṘ as a function of the angle θ. The characteristic points of this dynamics, maximum, minimum and zero points, are underlined and reported in the (y, z)-plane. Such specific points allow us to determine the optimal solution to manipulate in minimum time the system along the radial direction [51]. From a control point of view, this analysis allows to recover the singular trajectories, which are the paths where the control field is not equal, in absolute value, to the maximum bound [54]. We obtain here two types of singular controls: the horizontal one along the plane of equation z = z 0 = −γ 2(Γ−γ) , denoted S h and the vertical one along the z-axis denoted S v>0 and S v<0 , for z > z 0 and z < z 0 , respectively. Note that the plane of equation z = z 0 intersects the Bloch ball if In the standard case where T 1 ≥ T 2 , this leads to the condition 2T 1 ≥ 3T 2 . The time-optimal solution is thus the concatenation of bang arcs, denoted B where the intensity of the field is maximum (here infinite) and of singular arcs. We refer the reader to previous works for a thorough mathematical and physical analysis of this control process [51,53,54,55], but according to the constraint given below, the horizontal singular arcs can only be used if the condition 2T 1 ≥ 3T 2 is satisfied. An example is given in Fig. 3 for the saturation process where the goal is to reach as fast as possible the center of the Bloch ball. According to the values of the relaxation parameters T 1 and T 2 , we observe that two different control sequences are used. In the first situation, the optimal control law is composed of a bang pulse followed by two singular arcs, a horizontal and a vertical one. In the second case, an inversion pulse is used to reach the south pole before the application of a zero control along the z-axis. It can be shown that the analytical expression of the corresponding minimum times are given by: in the first and second cases, respectively.  Figure 2: (Color online) Contour plot of the trajectoriesṘ(θ) for two different sets of relaxation parameters. The red line is the line for whichṘ = 0. The blue and green lines are the set of points whereṙ is minimum, while the points of the yellow one correspond to a maximum. The small inserts display in the (y, z)-plane the trace of the corresponding sets of points. Circles of different radii are also plotted to help the reading. A straightforward generalization of this study allows us to design the optimal field bringing in minimum time the system from any initial point to any target state of the Bloch ball, leading to a complete description and understanding of this optimal control problem. Such a classification of all the optimal trajectories is called an optimal synthesis [3]. In each case, we recall that the time-optimal solution can only be the concatenation of different bang arcs and of singular paths along the z-axis or along the line z = z 0 .

Optimization of the SNR in the unbounded case
We have now all the tools in hand to optimize the SNR. This optimization is a non-standard and difficult control problem for which a cost functional Q has to be maximized, together with the determination of the initial S and final M points of the control period. Only the position of one of the two points needs to be computed since the S and the M points are connected by a free evolution. Following Refs. [49,50], we adopt a brute-force strategy to solve this problem, which consists in replacing this global optimization by a two-step procedure. For any M point of coordinates (y m , z m ) of the Bloch ball, we first determine as explained in Sec. 2.2 the time-optimal path going from S to M and we compute the corresponding value of the figure of merit Q. Only five types of control sequences can be optimal. According to the values of the relaxation parameters, only three different optimal syntheses are obtained, as displayed in Fig. 4. We observe the occurring of the Ernst ellipsoid, which is the set of points where R m = R s (the radii of the M and S points), i.e. the control sequence where only a δ-pulse is used to bring the system from S to M . In this unbounded situation, the functional Q can also be determined analytically as a function of the coordinates of M [49,50]. Mathematically, Q is a non-smooth but continuous surface. This surface is represented in Fig. 5 for the four cases of Fig. 4. A gradient method gives the maximum of this surface, which belongs to the Ernst ellipsoid. This point can be characterized by the well-known Ernst angle solution: which was first derived in the original paper by Ernst and his co-workers [48]. The angle θ denotes here the angle characterizing the δ-pulse of the Ernst sequence. This first result of optimality on the SNR per unit time paves the way in a near future to a systematic analysis of the optimization of this parameter in spin systems. This method can be generalized to the case of bounded control amplitudes, the use of crusher gradient pulses, to an inhomogeneous ensemble of spins with magnetic field broadening, or to any other experimental constraint. Other aspects will be discussed in the conclusion of the paper.

Numerical optimal control of the contrast in MRI
In this section, a numerical solution of the contrast optimization problem is presented. Contrast in MRI is a key parameter to visualize and interpret the anatomy and bio-chemical properties of potentially malign tissues. In clinical practice, contrast is usually empirically handled by tuning few acquisition parameters. This approach is however strongly radiologist dependent, only provides limited contrast combinations and offers no guarantee about the optimal nature of the created contrast. The objective of this study is to compute the radio-frequency (RF) magnetic pulse which optimally combines the effect of excitation and relaxation of the spins to produce the maximal signal difference between 2 tissues with different relaxation times. We recall that signal in MRI is directly linked to the norm of the transverse magnetization vector M , denoted M ⊥ = (M x , M y , 0). Geometric approaches have been proposed to solve this problem for different tissues for spins precessing at the nominal Larmor frequency [39,42,57]. For the application on a realistic MRI system, it is essential to account for magnet imperfections which induce various resonance frequency offsets. In this case, each isochromat is characterized by a magnetization vector M (ω) = (M x (ω), M y (ω), M z (ω)) whose dynamics is governed by the following differential system: where the offset term ω belongs to a given interval [ω min , ω max ]. From a numerical point of view, this control problem is solved by digitizing the distribution of the isochromats. However, the resulting dimensionality increase (several thou-sands of spins) makes the geometric solution inapplicable, and requires the use of numerical approaches to solve the problem. We also mention that a similar robustness approach can be applied for T 1 and T 2 deviations as well as control field inhomogeneities. This section details an example of a contrast experiment, from the optimal control formulation to the in vivo acquisition of a rat brain image.

Numerical Optimization
Several numerical schemes have been proposed to solve optimal control problems, including shooting methods [58,59], Krotov methods [43] and gradient ascent (or descent) based approaches [20,22]. This work illustrates a practical application of the gradient ascent pulse engineering (GRAPE) algorithm. The application of the algorithm requires the temporal discretization of the control field, where each time step represents an optimization variable. In its simplest version, this algorithm updates each control field time step by minimizing the user-defined cost function at each iteration, while fulfilling the constraints imposed by the Pontryagin Maximum Principle (PMP). The cost function minimization is performed by applying a step in the direction of the cost function gradient with respect to the control field. The reader is referred to [20,22] for explicit details on the gradient computation. Convergence is reached when the gradient and/or the step norms are bellow a user-defined threshold. Note that unlike the geometric methods, the global nature of the minimum cannot be claimed and global or local minimum is only reached within a certain tolerance. This tolerance is set in accordance with the experimental requirements so that additional improvements of the control field have a negligible impact on the experimental result. The control field initialization also impacts the nature of the minimum. Multiple trials or prior insights on the solution can be used to find the global minimum and to improve the convergence of the algorithm. The cost function (C) that optimizes the contrast must maximize the difference between the transverse magnetization norms of 2 samples a and b, at the end of the control time (t = t f ). The two species are characterized by different relaxation times T 1 and T 2 . It is arbitrarily chosen here that the signal of the samples a and b are respectively maximized and minimized. Note that other definitions of the contrast could be chosen. Static field inhomogeneities lead to different resonance frequency offsets which in turn result in different trajectories in the Bloch ball, that must be accounted for in the definition of the cost function. This can be done by sampling the frequency offset interval in N samples, where each point corresponds to a given trajectory. The cost function thus enforces all trajectories in the considered interval to have similar behaviors:

Example of the rat muscle/brain contrast
To illustrate the impact of optimal contrast pulses, a RF pulse is computed to optimize the contrast between the rat brain (minimize) and surrounding muscles (maximize). This example was chosen as a proof-of-concept because it cannot be created with standard contrast strategies, due to the short transverse relaxation time (T 2 ) of the muscle tissues [26]. The following experiment is performed in agreement with the UCBL's ethic committee on animal experimentation, on a 4.7 T Bruker MR system. Average relaxation times of both structures, as well as the offset inhomogeneity range are estimated by standard MR sequences [25]. Longitudinal (T 1 ) and transverse (T 2 ) relaxation times are respectively estimated to be [920, 60] ms for the brain and [1011, 30] ms for the muscle. The resonance offset range is of the order of 1 kHz. In order to handle slice selectivity, i.e. the location of the image section plane [26], standard MRI excitation schemes are used [26,60]. The optimal pulse is thus used as a contrast preparation pulse, implying that the contrast is prepared on the longitudinal axis (M z ). Only a slight change in the cost function is required to apply this modification [60]. The prepared magnetization is subsequently flipped into the transverse plane with a slice-selective π/2-pulse. We recall that the MRI signal is proportional to the amount of transverse magnetization, i.e. the tissue with the highest |M z | at the end of the preparation will produce higher intensity in the MR image. The magnitude of the resulting optimized pulse is shown in Fig. 6, together with the magnetization trajectories of the contrasted tissues. Magnetization trajectories are represented in a plane formed by the normalized transverse (M ⊥ ) and longitudinal (M z ) magnetization. In this figure, different trajectories represent a specific frequency offset. Notice how all trajectories of a given tissue start from the thermal equilibrium ( − → M 0 = (0, 0, 1)) and reach the same final magnetization state despite important trajectory disparities. This validates the pulse robustness against resonance offset variations. It can be observed that the brain is saturated at the end of the pulse, i.e. that trajectories reach the center of the Bloch sphere, while a significant amount of longitudinal magnetization is left for the muscle tissues. This implies that unlike muscle tissues, brain tissues will have a very low contribution to the acquired MRI signal. The corresponding acquired MR image is shown in Fig. 7. As expected, the average tissue intensity within the brain is much lower than the surrounding corpus callosum brain muscle Figure 7: Transverse section of a rat brain showing the result of the optimal pulse. The aim of the control is to saturate the brain while maximizing the intensity of the surrounding muscle tissues.
muscles. Interestingly, some internal brain structures can be distinguished. It corresponds to structures whose relaxation times slightly differ from the major brain components which are mostly gray matter. In particular, the corpus callosum, the largest white matter structure in the brain, can be observed as an hypersignal element. As white matter has a lower T 2 than gray matter [61], this validate the ability of the optimal pulse to enhance short T 2 tissues. This particular contrast opens interesting perspectives in both pre-clinical and clinical neuro-imaging applications since white matter usually appears in hyposignal due to its short T 2 . Such a contrast, keeping in mind its optimal nature, could help to improve the visualization of white matter in the human brain, which is a key factor in actual clinical challenges such as the effect of aging or Alzheimer disease [62]. This approach could also be used to highlight the differences between the structures of the brain.

Conclusion
In MRI, there exists an immense potential for improvement as well as for reduction of the energy deposit and time a patient has to spend in a scanner during an examination. In this work, we have presented two benchmark examples showing the efficiency of this technique. In the ideal case of a homogeneous ensemble of spin systems, we have demonstrated the optimality in the limit of unbounded controls of the Ernst angle solution, derived in the sixties, to enhance the SNR. Numerical optimization techniques were used to illustrate the possibility of maximizing the contrast of the image in MRI by using a particular pulse sequence. This approach was applied with success theoretically and experimentally on the brain of a rat. We emphasize that one of the main advantages of this contrast enhancement is its general character since the optimal control magnetic fields can be computed with standard routines published in the literature and implemented on a MRI scanner without requiring specific materials and process techniques. An important issue is the use of this method in clinics to optimize the contrast of human body imaging. We anticipate that this technique could be a complementary tool to contrast agents whose successful application has been an important aspect of the development of MRI in recent years. The joint use of these two methods could therefore limit the concentration of contrast agents needed for the imaging, which could be beneficial to the patient.
These two examples validate the potential of optimal control as an accurate pulse design tool to solve non-trivial practical problems in MRI. This ability can be extended to other contexts in MRI, as shown very recently in [63] with the introduction of optimal control pulses for magnetization phase control in MRI. There is a huge interest for accurate phase control in many important applications in MRI, such as diffusion imaging, thermometry and elastography, to improve the signal sensitivity to specific bio-markers.
The use of optimal control for these applications relies on further developments in robust optimal control analytic and numerical tools, accurate modeling of the MR physics and tailored acquisition sequences, that will further improve existing MR methods and approach the system's physical limits.