Abstract:
In a method for temperature control in MR-guided administration of ultrasound, ultrasound therapy is administered to a patient at an in vivo site by emitting focused ultrasound into the site at multiple foci with a multi-focus ultrasound therapy device. The temperature is monitored in a localized region of an examination subject in which the site is located firing said therapy in real-time by MR thermometry. From the MR thermometry, characteristics of the temperature distribution in the monitored region of the examination subject are automatically identified. Temperature control is implemented by regulating the energy output of the ultrasound therapy device, according to a rapidly converging master equation.

Description:
BACKGROUND OF THE INVENTION 
       [0001]    1. Field of the Invention 
         [0002]    The present invention concerns a method for controlling temperature in a patient undergoing ultrasound therapy under the guidance of magnetic resonance thermometry, and in particular to a rapidly convergent, multi-input, multi-output, non-parametric temperature controller for implementing such a method. 
         [0003]    2. Description of the Prior Art 
         [0004]    The localized application of ultrasound to treat certain types of local cancer often takes place with the patient being located, during the therapy, in a magnetic resonance imaging system, so that the temperature of the patient can be monitored by magnetic resonance thermometry. Magnetic resonance thermometry is a known manner of monitoring the temperature of an examination subject while the subject is located in a magnetic resonance data acquisition unit (scanner), and is based on the known phenomenon that the magnetic resonance signals emitted by nuclear spins in the examination subject, which have been forced out of a steady-state condition by the application of energy thereto, are temperature-dependent. The precessing nuclear spins that have been forced out of their equilibrium state emit signals that are detected and are used for conventional magnetic resonance imaging. The temperature-dependency of these signals also allows a temperature map of the examination subject to be generated, that very precisely shows temperature variations along multiple axes. 
         [0005]    High intensity ultrasound is effective for use in treating localized cancers or other pathologies. For this purpose, the high intensity focused ultrasound (HIFU) is administered. 
         [0006]    In order to steer and focus the therapeutic ultrasound, such an ultrasound head or device typically has multiple ultrasound emitters (transducers) arranged and operated in an array. 
         [0007]    In general terms, magnetic resonance guiding or monitoring of such therapy takes place by administering HIFU to the patient in the examination region of a magnetic resonance scanner. Magnetic resonance data are acquired in a known manner from the patient in the scanner while the ultrasound therapy is in progress. In real-time with the therapy, magnetic resonance thermometry images of the patient are generated, showing the temperature distribution within a designated region of the patient. The resulting MR thermometry image can be shown on a monitor in real-time during the therapy for visual review and manual control of the therapy by a physician or a technician, or known image processing techniques can be used in order to generate appropriate extractions of information from the thermometry image for use in automatic control of the therapy. 
         [0008]    An overview of MR-guided focused ultrasound hyperthermia is described in “Hyperthermia by MR-guided Focused Ultrasound: Accurate Temperature Control Based on Fast MRI and A Physical Model of Local Energy Deposition and Heat Conduction,” Salomir et al., Magnetic Resonance in Medicine, Vol. 43 (2000) pages 342-347. Various automatic control techniques are described, for example, in “Automatic Spatial and Temporal Temperature Control for MR-Guided Focused Ultrasound Using Fast 3D MR Thermometry and Multispiral Trajectory of the Focal Point,” Mougenot et al., Magnetic Resonance in Medicine, Vol. 52 (2004) pgs 1005-1015; “Three-Dimensional Spatial and Temporal Temperature Control with MR Thermometry-Guided Focused Ultrasound (MRgHIFU),” Mougenot et al., Magnetic Resonance in Medicine, Vol. 61 (2009) pgs 603-614 and “Curvilinear Transurethral Ultrasound Applicator for Selective Prostrate Thermal Therapy,” Ross et al., Medical Physics, Vol. 32, No. 6 (2005) pgs 1555-1565. A prevalent technique for temperature control in such therapy is to make use of a proportional-integral-derivative (PID) controller. In general terms, PID controller calculates an error value as the difference between a measured process variable and a desired reference point. As indicated by its name, a PID controller makes use of three feedback reaction terms, respectively referred to as proportional, integral and derivative values. The proportional value is dependent on the current error, the integral value represents an accumulation of past error, and the derivative value represents a prediction of future errors, based on a current rate of change. Typically, a weighted sum of these three values is used to adjust the process in question, by adjusting some type of control element that has the capability of changing the monitored parameter of the process in question. A study of the use of PID temperature control in the context of MRI-guided phased-array contact ultrasound is provided in “Endocavitary Thermal Therapy by MRI-guided Phased-Array Contact Ultrasound: Experimental and Numerical Studies on the Multi-input Single-Output PID Temperature Controller&#39;s Convergence and Stability,” Salomir et al., Medical Physics, Vol. 36, No. 10 (2009) pgs 4726-4741. 
         [0009]    The PID-based controller in general exhibit good stability and robustness to noise, but its convergence is rather slow. When the number of temperature sampling points during the controlled sonication is low, the PID algorithm does not have sufficient time to adjust the power level so as to compensate for local variability in heat deposition or tissue cooling by diffusion/perfusion. The condition is worsened when the total time base is short, at the scale of the tissue response time to an elementary heat source. Moreover, when the static tuning of the physical parameters in the underlying model is far from the true parameters, overshooting of the temperature curves occurs with transient but significant and possibly longer heating above the prescribed temperature. 
       SUMMARY OF THE INVENTION 
       [0010]    It is an object of the present invention to provide an improved temperature control algorithm for use in MRI-guided phased-array contact ultrasound that is particularly suited for cycled sonication along a volumetric pattern that is repeated for each MR dynamic scan. 
         [0011]    It is a further object of the present invention to provide such a temperature control algorithm that delivers an equivalent thermal dose in every focus in a time interval of 20 to 30 times the duration of the MR-thermometry measurement. 
         [0012]    The above object is achieved in accordance with the present invention by a method for temperature control in MR-guided administration of ultrasound wherein the temperature in a localized region of an examination subject is monitored in real-time by MR thermometry and, from the MR thermometry, characteristics of the temperature distribution in the monitored region of the examination subject are automatically identified and are employed in a computerized processor to implement temperature control, by regulating the ultrasound therapy device. 
         [0013]    The temperature controller is designed in two steps: 1) the classic PID equation is modified for improved convergence, and 2) the algorithm is rendered non-parametric by considering the relative offset of temperature of a given focus with respect to the average temperature calculated over all the foci. 
         [0014]    The classic equation of a PID controller is modified into: 
         [0000]    
       
         
           
             
               P 
                
               
                 ( 
                 t 
                 ) 
               
             
             = 
             
               
                 1 
                 
                   μ 
                    
                   
                     ( 
                     
                       
                         r 
                         → 
                       
                       o 
                     
                     ) 
                   
                 
               
                
               
                 { 
                 
                   
                     
                        
                       
                         Θ 
                          
                         
                           ( 
                           t 
                           ) 
                         
                       
                     
                     
                        
                       t 
                     
                   
                   + 
                   
                     
                       D 
                        
                       
                         ( 
                         
                           
                             r 
                             → 
                           
                           o 
                         
                         ) 
                       
                     
                     · 
                     
                       
                         ∇ 
                         2 
                       
                        
                       
                         T 
                          
                         
                           ( 
                           
                             
                               
                                 r 
                                 → 
                               
                               o 
                             
                             , 
                             t 
                           
                           ) 
                         
                       
                     
                   
                   + 
                   
                     w 
                     · 
                     t 
                     · 
                     
                       Δ 
                        
                       
                         ( 
                         
                           
                             
                               r 
                               → 
                             
                             o 
                           
                           , 
                           t 
                         
                         ) 
                       
                     
                   
                   + 
                   
                     w 
                     · 
                     
                       Ω 
                        
                       
                         ( 
                         
                           
                             
                               r 
                               → 
                             
                             o 
                           
                           , 
                           t 
                         
                         ) 
                       
                     
                   
                 
                 } 
               
             
           
         
       
     
         [0000]    wherein
   μ({right arrow over (r)} o ) is the ultrasound absorption point {right arrow over (r)} o      D({right arrow over (r)} o ) is the heat diffusivity point {right arrow over (r)} o      Θ(t) is the prescribed (or target) evolution of the temperature,   T({right arrow over (r)} o ,t) is the experimental (measured) temperature point {right arrow over (r)} o  time t,   Δ({right arrow over (r)} o ,t)=Θ(t)−T({right arrow over (r)} o ,t) is the current error at location {right arrow over (r)} o ; and   
 
         [0000]    
       
         
           
             
               Ω 
                
               
                 ( 
                 
                   
                     
                       r 
                       → 
                     
                     o 
                   
                   , 
                   t 
                 
                 ) 
               
             
             = 
             
               
                 ∫ 
                 0 
                 t 
               
                
               
                 
                   [ 
                   
                     
                       Θ 
                        
                       
                         ( 
                         
                           t 
                           ′ 
                         
                         ) 
                       
                     
                     - 
                     
                       T 
                        
                       
                         ( 
                         
                           
                             
                               r 
                               → 
                             
                             o 
                           
                           , 
                           
                             t 
                             ′ 
                           
                         
                         ) 
                       
                     
                   
                   ] 
                 
                  
                 
                     
                 
                  
                 
                    
                   
                     t 
                     ′ 
                   
                 
               
             
           
         
       
     
         [0000]    is the integrated error at location {right arrow over (r)} o . 
         [0020]    The characteristic differential equation in accordance with the invention is: 
         [0000]    
       
         
           
             
               
                 
                   
                      
                     2 
                   
                    
                   Ω 
                 
                 
                    
                   
                     t 
                     2 
                   
                 
               
               + 
               
                 w 
                 · 
                 t 
                 · 
                 
                   
                      
                     Ω 
                   
                   
                      
                     t 
                   
                 
               
               + 
               
                 w 
                 · 
                 Ω 
               
             
             = 
             0 
           
         
       
     
         [0021]    The time domain convergence solution is: 
         [0000]    
       
         
           
             
               Ω 
                
               
                 ( 
                 t 
                 ) 
               
             
             ∝ 
             
               exp 
                
               
                 ( 
                 
                   - 
                   
                     
                       w 
                       · 
                       
                         t 
                         2 
                       
                     
                     2 
                   
                 
                 ) 
               
             
           
         
       
     
         [0022]    The convergence time is (2/w) 1/2,  typically for value being √{square root over (8)} to √{square root over (10)} sampling points. A value for the negative reaction w can be selected in a range between 0.01 divided by the square of the sampling time, and 10.0 divided by the square of the sampling time. Recommended values for the negative reaction (denoted w) are in a range between 0.20 divided by the square of the sampling time, and 0.25 divided by the square of the sampling time. 
         [0023]    In accordance with the present invention, a feedback control algorithm is employed wherein the negative reaction from the actual error (e.g., deviation or offset) of the measured temperature relative to a target is increased linearly with time. The temperature controller in accordance with the invention is based on a second order differential equation that is linear but with non-constant coefficients, wherein the solution is in the form of exp(−w·t 2 /2), wherein t=time and w is a constant. Such a solution is more stable and rapidly converging than the solution of the PID equation, of the form (1+a·t)·exp(−q·t/2), wherein a and q are different constants. 
         [0024]    The reference value of temperature, that is used to compare if the local temperature at a given focus has a positive or negative offset at a given time point, is defined by the following expression (instead of a continuous predefined target curve as it used to be the case for previous PID-based reports): 
         [0000]    
       
         
           
             
               Θ 
               n 
             
             = 
             
               { 
               
                 
                   
                     
                       
                         1 
                         
                           N 
                           f 
                         
                       
                        
                       
                         
                           ∑ 
                           
                             k 
                             = 
                             1 
                           
                           
                             N 
                             f 
                           
                         
                          
                         
                             
                         
                          
                         
                           T 
                           
                             k 
                             , 
                             n 
                           
                         
                       
                     
                   
                   
                     
                       before 
                        
                       
                           
                       
                        
                       the 
                        
                       
                           
                       
                        
                       steady 
                        
                       
                           
                       
                        
                       state 
                        
                       
                           
                       
                        
                       regimen 
                     
                   
                 
                 
                   
                     
                       Θ 
                       top 
                     
                   
                   
                     
                       during 
                        
                       
                           
                       
                        
                       the 
                        
                       
                           
                       
                        
                       steady 
                        
                       
                           
                       
                        
                       state 
                        
                       
                           
                       
                        
                       regimen 
                     
                   
                 
               
             
           
         
       
     
         [0025]    Note that Θ n  is the spatial average of temperature values T k,n  over the N f  heated points (i.e. the sonicated foci) before the steady-state regimen was established and that Θ n  becomes identical to the prescribed target temperature after the steady-state regimen has been reached (Θ top , i.e. the flattop segment of the target curve). The temperature control algorithm disclosed herein aims to minimize or to null the relative offsets 
         [0000]    
       
         
           
             
               Δ 
               
                 k 
                 , 
                 n 
               
             
             = 
             
               1 
               - 
               
                 
                   
                     T 
                     
                       k 
                       , 
                       n 
                     
                   
                   
                     Θ 
                     n 
                   
                 
                 . 
               
             
           
         
       
     
         [0000]    The offset and the integral of the offset are calculated in each control point and these calculi are used to modulate the amount of energy delivered to each focus. The method in accordance with the present invention thereby homogenizes the temperature history of all heated locations. 
         [0026]    The convergence of the controller is considered powerful enough to perform accurately even without a priori knowledge on heated tissue parameters. Therefore any tissue dependent parameters are discarded from the controller&#39;s equation. This condition simplifies the practical procedure of treatment and avoids the need for identification of parameters (subject to errors, time-consuming etc). The controller is therefore non-parametric, because no a priori hypothesis is made on the ultrasound absorption, the beam propagation (other than the fact that a focus must exist), heat diffusion, perfusion rate, etc. In other words, the master equation that is used to adjust the energy to be delivered at a given location uses no a priori parameter. 
         [0027]    The feedback compensation of delivered energy at each focus is updated for each new dynamic scan, based on either the relative offset in temperature at that location as compared to the average temperature overall foci, during the rising portion of the heating curve, or the relative offset in temperature at that location as compared to the prescribed target value, during the steady state portion of the heating curve. 
         [0028]    The method in accordance with the invention achieves faster convergence for a spatially homogenous temperature considered over all of the control points, namely approximately 10-15 temporal sampling points. The method according to the invention achieves higher accuracy, which means a reduced amplitude of any overshooting or undershooting of the experimental temperature curve, even when confronted by significant changes in the energy deposition from one heating location to another, or when confronted with dynamic changes of tissue properties. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0029]      FIG. 1  is an image of a multi-focus pattern for volumetric sonication of a subject of the type that is controlled in accordance with the present invention. 
           [0030]      FIG. 2  schematically illustrates the initial amount of energy E k ,1 prescribed for each focus. 
           [0031]      FIG. 3  is a comparison of the controller modification in accordance with the present invention compared to conventional PID control. 
           [0032]      FIG. 4  is a comparison of the controller modification in accordance with the present invention compared to parametric and continuous time control. 
           [0033]      FIG. 5  is a graph illustrating the threshold limitation of the negative reaction of the proportional term. 
           [0034]      FIG. 6  is a graph illustrating the hyperbolic tangent limitation of the negative reaction of the proportional term. 
           [0035]      FIG. 7  illustrates temperature maps on magnitude images for three orthogonal planes during magnetic resonance guided high intensity focused ultrasound (HIFU) sonication, with control in accordance with the present invention. 
           [0036]      FIG. 8  illustrates the temperature evolution at 5 individual locations indicated by pixel coordinates, with temperature control in accordance with the present invention. 
           [0037]      FIGS. 9A ,  9 B,  9 C and  9 D are surface representations of a 1 D temperature profile scrolled over time, for the patient for which the temperature maps of  FIG. 7  were generated. 
           [0038]      FIG. 10  schematically illustrates acoustic energy delivered per target focus in accordance with the present invention. 
           [0039]      FIG. 11  is a temperature map overlaid on a magnitude image in the coronal plane during six-focus sonication in a rabbit thigh, in accordance with the present invention. 
           [0040]      FIG. 12  shows the temperature evolution in three of the six foci in the sonication of the rabbit thigh shown in  FIG. 11 . 
           [0041]      FIG. 13  schematically illustrates the basic components for implementing the method in accordance with the present invention. 
       
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENTS 
       [0042]      FIG. 13  schematically illustrates the basic components of MR-guided ultrasound therapy, such as high intensity focused ultrasound (HIFU). For this purpose, a patient P is placed in a magnetic resonance imaging data acquisition unit  1 , which is operated in a known manner by control and imaging components  2 . A high intensity focused ultrasound device, such as a multi-focus device, and is operated to implement sonication of the patient P at a localized site at which, for example, cancer tissue had been determined to be located. The MR data acquired from the patient P in a localized region in which the ultrasound has been administered are processed in a known manner in an MR thermometry unit  4 , in order to obtain an electronic representation of the temperature distribution in the localized region of the patient P produced by operation of the ultrasound device  3 . An image representing this temperature distribution can be displayed at a monitor (not shown) for visual examination and use by a technician or a physician. Additionally, relevant information can be electronically obtained or extracted from the thermometry representation, and supplied to pay controller  5  that operates the ultrasound device  3  by controlling ultrasound emission from the respective foci thereof. The controller  5 , in accordance with the invention, is designated as a P 2 ID controller, in order to distinguish it from the conventional PID controller that is commonly used. 
         [0043]    The volumetric sonication that is implemented in accordance with the present invention is accomplished by rapidly iterated (activated) foci, the activation taking place at the time scale of heat diffusion in the tissue being insonified. The thermal build-up (distribution) is constructed as a spatially continuous function, such as by repeated cycles of sonication along the same sonication pattern, with one cycle duration being equal to the cycle of one MR dynamic scan (i.e., data acquisition or measurement cycle). The number of foci forming the pattern is N f ≧2, and they are sonicated nearly simultaneously using a phased-array electronic steering of the beam. The distance between foci (inter-focus spatial gap) must be larger than the intrinsic size of the focal spot for individual identification, but should not be so large as to permit the formation of a continuous thermal build-up by exploiting heat diffusion. 
         [0044]      FIG. 1  illustrates a multi-focus pattern for volumetric pattern having ten individual foci N f =10 as an example. Index k references the foci 1:N f and E k,1  is the initial amount of energy prescribed to each focus (scan=1), compensated for (a-priori) the energy loss by steering and for the preferential heat diffusion at borders. During the procedure, E k,1  will be dynamically modulated by the temperature controller at each scan number n to become E k,n =E k,1 +δE k (n), to insure similar temperature elevation at each foci, as schematically illustrated in  FIG. 2 . 
         [0045]    The principles of treatment are as follows:
       1. a sonication-cycle represents the sequential displacement of the focus trough all prescribed foci; the coordinates of the prescribed foci are fixed trough the full duration of the treatment;   2. the foci are successively sonicated by electronic steering, such as one sonication-cycle is completed in the same duration as that of one measurement (or, equivalent, with the temporal resolution of MR thermometry);   3. the controller tries to compensate the local temperature deviation at each focus as compared to the average temperature over all foci   4. the controller updates the energy (product power×time per point) to be delivered each time a new MR data set becomes available;   5. the dwell time (or: sampling time) of the feedback loop is the same for each focus location and is therefore set by the MR thermometry:       
 
         [0000]      sampling_time=MRT_temporal_resolution       6. the distance between foci must be larger than the intrinsic size of the focal spot (for individual identification) but not too large, to enable the formation of a continuous thermal build up by exploiting heat diffusion; a value of 4 mm gap is suggested as optimal   7. the number of cycles is decided a-priori or adjusted on the fly based on cumulated thermal dose;         
         [0053]    In  FIGS. 3 and 4 , for easier understanding, the controller is considered as modulating the power level at each focus. In reality, what is important is the energy delivered at each focus location. The first choice is to modulate the duration of sonication, and to additionally reinforce the delivered energy by changing the power level (if the modulation of duration was not enough to assure local control of temperature). Modulating the energy deposition at the sampling point n for the focus k by the amount δE k (n) relative to the initial value E k,1  is achieved based on following double-branch equation: 
         [0000]    
       
         
           
             
               
                 δ 
                  
                 
                     
                 
                  
                 
                   
                     E 
                     k 
                   
                    
                   
                     ( 
                     n 
                     ) 
                   
                 
               
               
                 E 
                 
                   k 
                   , 
                   1 
                 
               
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           w 
                           * 
                         
                         · 
                         n 
                         · 
                         
                           ( 
                           
                             1 
                             - 
                             
                               
                                 T 
                                 
                                   k 
                                   , 
                                   n 
                                 
                               
                               
                                 Θ 
                                 n 
                               
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           w 
                           * 
                         
                         · 
                         
                           Ω 
                           
                             k 
                             , 
                             n 
                           
                         
                       
                     
                   
                   
                     
                       
                         if 
                          
                         
                             
                         
                          
                         n 
                       
                       ≤ 
                       
                         n 
                         switch 
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             
                               - 
                               1 
                             
                             + 
                             
                               tanh 
                                
                               
                                 ( 
                                 
                                   α 
                                   · 
                                   
                                     n 
                                     switch 
                                   
                                 
                                 ) 
                               
                             
                             + 
                             
                               
                                 w 
                                 * 
                               
                               · 
                               n 
                               · 
                             
                           
                         
                       
                       
                         
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   
                                     T 
                                     
                                       k 
                                       , 
                                       n 
                                     
                                   
                                   
                                     Θ 
                                     top 
                                   
                                 
                               
                               ) 
                             
                             + 
                             
                               
                                 w 
                                 * 
                               
                               · 
                               
                                 Ω 
                                 
                                   k 
                                   , 
                                   n 
                                 
                               
                             
                           
                         
                       
                     
                   
                   
                     else 
                   
                 
               
             
           
         
       
     
         [0000]    wherein tan h is the analytical function hyperbolic tangent, n switch  denotes the sampling point where the controller switches irreversibly from the rapid temperature elevation regimen to the steady state regimen, α is a positive number in general smaller than 0.1 (recommended to be approximately 0.025), T k,n  is the actual temperature at the sonicated focus k for the measurement n, Θ n  is the average value of measured temperatures T k,n  over the sonicated foci for measurement n, Θ top  is the flattop target value of steady state temperature elevation, w* is the negative feedback parameter a-dimensional and scaled to the square of the temporal resolution (τ) as w*=w/τ 2 , and Ω k,n  is the cumulative relative offset (or, error) of the measured temperature until the measurement n: 
         [0000]    
       
         
           
             
               Ω 
               
                 k 
                 , 
                 n 
               
             
             = 
             
               { 
               
                 
                   
                     
                       
                         ∑ 
                         
                           
                             n 
                             ′ 
                           
                           = 
                           1 
                         
                         n 
                       
                        
                       
                           
                       
                        
                       
                         ( 
                         
                           1 
                           - 
                           
                             
                               T 
                               
                                 k 
                                 , 
                                 
                                   n 
                                   ′ 
                                 
                               
                             
                             
                               Θ 
                               
                                 n 
                                 ′ 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                   
                     
                       
                         if 
                          
                         
                             
                         
                          
                         n 
                       
                       ≤ 
                       
                         n 
                         switch 
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           ∑ 
                           
                             
                               n 
                               ′ 
                             
                             = 
                             1 
                           
                           
                             n 
                             switch 
                           
                         
                          
                         
                             
                         
                          
                         
                           ( 
                           
                             1 
                             - 
                             
                               
                                 T 
                                 
                                   k 
                                   , 
                                   
                                     n 
                                     ′ 
                                   
                                 
                               
                               
                                 Θ 
                                 
                                   n 
                                   ′ 
                                 
                               
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           ∑ 
                           
                             
                               n 
                               ′ 
                             
                             = 
                             
                               
                                 n 
                                 switch 
                               
                               + 
                               1 
                             
                           
                           n 
                         
                          
                         
                             
                         
                          
                         
                           ( 
                           
                             1 
                             - 
                             
                               
                                 T 
                                 
                                   k 
                                   , 
                                   
                                     n 
                                     ′ 
                                   
                                 
                               
                               
                                 Θ 
                                 top 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                   
                     
                       
                         if 
                          
                         
                             
                         
                          
                         n 
                       
                       &gt; 
                       
                         n 
                         switch 
                       
                     
                   
                 
               
             
           
         
       
     
         [0054]    During the steady state regime (starting at the sampling point n switch +1), the term tan h(α·n switch ) takes into account the predicted amount of energy (relative to the initial value at n=1) necessary to maintain the steady state (constant temperature) in homogeneous tissue; inhomogeneities are managed for by the feedback terms proportional and integral. Based on the behavior of the hyperbolic-tangent function: a). if the tissue reaches rapidly the target top temperature, this means the maintaining energy is low; b). when n switch  increases to infinity, the full energy is necessary to maintain the steady state. 
         [0055]    The controller algorithm decides when to switch to the state-state computing equation, meaning when the sampling point is declared equal to n switch . 
         [0056]    This is done in a predictive manner. 
         [0057]    Sampling point “n” becomes “n switch ” if following condition is fulfilled: 
         [0000]      [θ n +(θ n −θ n−1 )]&gt;(θ top −δ) or θ top −2·θ n +θ n−1 &lt;δ
 
         [0058]    Textually, this means that the predicted temperature elevation after one more sonication cycle will exceed the transition level. Here δ=0.5° C. The predicted temperature elevation is the actual value plus the last increment between the previous and actual cycle. 
         [0059]    The transition is decided irreversibly. Once the master equation was switched to the steady-state form, this one will be used until the end of the active heating. 
         [0060]    According to equations provided above, the proportional term brings a negative reaction that linearly increases with time. Because the temperature data are contaminated by experimental noise, over weighted negative reaction can drive the controller into instability. 
         [0061]    Moreover, it is important that the negative reaction increases linearly during the “learning” stage, which corresponds to the time interval of the first 10 to 15 dynamics. 
         [0062]    Once the temperature reached the steady-state level (“cruise altitude”), the heating device only needs to compensate for local cooling by diffusion and perfusion, and the negative reaction can be stabilized. 
         [0063]    Therefore we modify the master equation for the controller output, by limiting the negative reaction of the proportional term (n≦n max ) as follows: 
         [0064]    A. Threshold limitation of the negative reaction of proportional term 
         [0000]    
       
         
           
             
               
                 δ 
                  
                 
                     
                 
                  
                 
                   
                     E 
                     k 
                   
                    
                   
                     ( 
                     n 
                     ) 
                   
                 
               
               
                 E 
                 
                   k 
                   , 
                   1 
                 
               
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           w 
                           * 
                         
                         · 
                         
                           min 
                            
                           
                             ( 
                             
                               n 
                               , 
                               
                                 n 
                                 max 
                               
                             
                             ) 
                           
                         
                         · 
                         
                           ( 
                           
                             1 
                             - 
                             
                               
                                 T 
                                 
                                   k 
                                   , 
                                   n 
                                 
                               
                               
                                 Θ 
                                 n 
                               
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           w 
                           * 
                         
                         · 
                         
                           Ω 
                           
                             k 
                             , 
                             n 
                           
                         
                       
                     
                   
                   
                     
                       
                         if 
                          
                         
                             
                         
                          
                         n 
                       
                       ≤ 
                       
                         n 
                         switch 
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             
                               - 
                               1 
                             
                             + 
                             
                               tanh 
                                
                               
                                 ( 
                                 
                                   
                                     α 
                                     * 
                                   
                                    
                                   
                                     n 
                                     switch 
                                   
                                 
                                 ) 
                               
                             
                             + 
                             
                               
                                 
                                   w 
                                   * 
                                 
                                 · 
                                 min 
                               
                                
                               
                                 
                                   ( 
                                   
                                     n 
                                     , 
                                     
                                       n 
                                       max 
                                     
                                   
                                   ) 
                                 
                                 · 
                               
                             
                           
                         
                       
                       
                         
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   
                                     T 
                                     
                                       k 
                                       , 
                                       n 
                                     
                                   
                                   
                                     Θ 
                                     top 
                                   
                                 
                               
                               ) 
                             
                             + 
                             
                               
                                 w 
                                 * 
                               
                               · 
                               
                                 Ω 
                                 
                                   k 
                                   , 
                                   n 
                                 
                               
                             
                           
                         
                       
                     
                   
                   
                     else 
                   
                 
               
             
           
         
       
     
         [0065]    The plot for a recommended value n max =12 is shown in  FIG. 5 . 
         [0066]    B. Hyperbolic tangent limitation of the negative reaction of the proportional term 
         [0000]    
       
         
           
             
               
                 δ 
                  
                 
                     
                 
                  
                 
                   
                     E 
                     k 
                   
                    
                   
                     ( 
                     n 
                     ) 
                   
                 
               
               
                 E 
                 
                   k 
                   , 
                   1 
                 
               
             
             = 
             
               { 
               
                 
                   
                     
                       
                         
                           w 
                           * 
                         
                         · 
                         
                           n 
                           max 
                         
                         · 
                         
                           tanh 
                            
                           
                             ( 
                             
                               n 
                               
                                 n 
                                 max 
                               
                             
                             ) 
                           
                         
                         · 
                         
                           ( 
                           
                             1 
                             - 
                             
                               
                                 T 
                                 
                                   k 
                                   , 
                                   n 
                                 
                               
                               
                                 Θ 
                                 n 
                               
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           w 
                           * 
                         
                         · 
                         
                           Ω 
                           
                             k 
                             , 
                             n 
                           
                         
                       
                     
                   
                   
                     
                       
                         if 
                          
                         
                             
                         
                          
                         n 
                       
                       ≤ 
                       
                         n 
                         switch 
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             
                               - 
                               1 
                             
                             + 
                             
                               tanh 
                                
                               
                                 ( 
                                 
                                   
                                     α 
                                     * 
                                   
                                    
                                   
                                     n 
                                     switch 
                                   
                                 
                                 ) 
                               
                             
                             + 
                             
                               
                                 w 
                                 * 
                               
                               · 
                               
                                 n 
                                 max 
                               
                               · 
                               
                                 tanh 
                                  
                                 
                                   ( 
                                   
                                     n 
                                     
                                       n 
                                       max 
                                     
                                   
                                   ) 
                                 
                               
                               · 
                             
                           
                         
                       
                       
                         
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   
                                     T 
                                     
                                       k 
                                       , 
                                       n 
                                     
                                   
                                   
                                     Θ 
                                     top 
                                   
                                 
                               
                               ) 
                             
                             + 
                             
                               
                                 w 
                                 * 
                               
                               · 
                               
                                 Ω 
                                 
                                   k 
                                   , 
                                   n 
                                 
                               
                             
                           
                         
                       
                     
                   
                   
                     else 
                   
                 
               
             
           
         
       
     
         [0067]    The plot for a recommended value n max =12 is shown in  FIG. 6 . 
         [0068]      FIG. 7  shows temperature maps overlaid on magnitude images for three orthogonal planes, during MRgHIFU sonication under active control (using the disclosed controller). The time point for each map is indicated on left. During steady state regime, the maximum (central) temperature is constant (±noise) but the edges of the thermal build up are expanding and smoothing because of diffusion (hence the interest to treat as fast as possible, provided a rapidly convergent controller). 
         [0069]      FIG. 8  shows the temperature evolution at five individual locations (pixel coordinates indicated lower right for the same experiment as  FIG. 7 . The similar temperature history of each focus, and the rising-temperature regime and further the steady-state regime can be seen. Foci from the pattern edge are cooling faster than central ones after sonication ended. 
         [0070]      FIGS. 9A ,  9 B,  9 C and  9 D are surface representations of a 1D temperature profile scrolled over time, for the same experiment as for  FIGS. 7 and 8 .  FIG. 10  shows acoustic energies delivered per target focus, as recorded by the SW interface. The edge foci received significantly more energy while undergoing similar temperature regime as central foci. This justifies the need for feedback control. 
         [0000]    
       
         
               
             
               
               
               
               
               
             
               
               
               
               
               
             
           
               
                 TABLE 
               
             
             
               
                   
               
               
                 Listing of average temperature elevation of the 10  
               
               
                 sonicated foci over time during the actively controlled  
               
               
                 heating in one typical experiment. The inter-foci  
               
               
                 standard deviation (at a given time-point) is indicated, 
               
               
                 i.e. the degree of spatial homogeneity, and the  
               
               
                 spatially-integrated acoustic energy for each cycle. 
               
               
                 The prescribed (or, target) temperature  
               
               
                 elevation was 20° C. above the initial baseline. 
               
             
          
           
               
                   
                   
                   
                   
                 Total energy  
               
               
                   
                   
                 Average 
                 STD over 
                 per scan equivalent, 
               
               
                   
                 Scan/ 
                 temperature 
                 foci 
                 per sonication cycle 
               
               
                   
                 or cycle 
                 [° C.] 
                 [° C.] 
                 [J] 
               
               
                   
                   
               
             
          
           
               
                   
                 1 
                 2.92 
                 0.80 
                 465 
               
               
                   
                 2 
                 3.66 
                 0.81 
                 465 
               
               
                   
                 3 
                 5.78 
                 1.21 
                 508 
               
               
                   
                 4 
                 8.08 
                 1.12 
                 524 
               
               
                   
                 5 
                 10.22 
                 0.88 
                 523 
               
               
                   
                 6 
                 11.97 
                 0.71 
                 514 
               
               
                   
                 7 
                 13.91 
                 0.62 
                 504 
               
               
                   
                 8 
                 15.09 
                 0.58 
                 502 
               
               
                   
                 9 
                 16.88 
                 0.65 
                 491 
               
               
                   
                 10 
                 18.20 
                 0.68 
                 337 
               
               
                   
                 11 
                 19.18 
                 0.70 
                 192 
               
               
                   
                 12 
                 19.92 
                 0.61 
                 154 
               
               
                   
                 13 
                 19.76 
                 0.58 
                 166 
               
               
                   
                 14 
                 19.45 
                 0.67 
                 192 
               
               
                   
                 15 
                 19.56 
                 0.78 
                 189 
               
               
                   
                 16 
                 19.66 
                 0.48 
                 180 
               
               
                   
                 17 
                 19.80 
                 0.83 
                 176 
               
               
                   
                 18 
                 19.70 
                 0.57 
                 181 
               
               
                   
                 19 
                 19.82 
                 0.45 
                 173 
               
               
                   
                 20 
                 20.03 
                 0.50 
                 159 
               
               
                   
                 21 
                 19.99 
                 0.48 
                 162 
               
               
                   
                 22 
                 20.02 
                 0.41 
                 160 
               
               
                   
                 23 
                 20.04 
                 0.48 
                 158 
               
               
                   
                 24 
                 20.14 
                 0.33 
                 151 
               
               
                   
                 25 
                 20.09 
                 0.35 
                 153 
               
               
                   
                 26 
                 20.19 
                 0.51 
                 147 
               
               
                   
                 27 
                 20.24 
                 0.30 
                 141 
               
               
                   
                 28 
                 20.18 
                 0.37 
                 144 
               
               
                   
                 29 
                 20.41 
                 0.38 
                 126 
               
               
                   
                 30 
                 20.57 
                 0.31 
                 112 
               
               
                   
                 31 
                 20.16 
                 0.17 
                 139 
               
               
                   
                 32 
                 20.33 
                 0.23 
                 126 
               
               
                   
                 33 
                 20.19 
                 0.34 
                 136 
               
               
                   
                 34 
                 20.59 
                 0.36 
                 105 
               
               
                   
                 35 
                 20.16 
                 0.32 
                 133 
               
               
                   
                 36 
                 20.39 
                 0.35 
                 116 
               
               
                   
                 37 
                 20.40 
                 0.16 
                 113 
               
               
                   
                 38 
                 20.05 
                 0.25 
                 137 
               
               
                   
                 39 
                 20.56 
                 0.32 
                 101 
               
               
                   
                 40 
                 20.16 
                 0.21 
                 125 
               
               
                   
                   
               
             
          
         
       
     
         [0071]      FIG. 11  shows temperature maps overlaid on magnitude images in coronal plane during 6-focus sonication in rabbit thigh in accordance with the invention.  FIG. 12  shows temperature evolution shown for (of the six) foci. 
         [0072]    Although modifications and changes may be suggested by those skilled in the art, it is the intention of the inventors to embody within the patent warranted hereon all changes and modifications as reasonably and properly come within the scope of their contribution to the art.