WO2016049013A1 - Local temperature rise constrained radio frequency pulse design in parallel transmission mri - Google Patents

Local temperature rise constrained radio frequency pulse design in parallel transmission mri Download PDF

Info

Publication number
WO2016049013A1
WO2016049013A1 PCT/US2015/051425 US2015051425W WO2016049013A1 WO 2016049013 A1 WO2016049013 A1 WO 2016049013A1 US 2015051425 W US2015051425 W US 2015051425W WO 2016049013 A1 WO2016049013 A1 WO 2016049013A1
Authority
WO
WIPO (PCT)
Prior art keywords
temperature
recited
computed
temperature rise
sar
Prior art date
Application number
PCT/US2015/051425
Other languages
French (fr)
Inventor
Nicolas Boulant
Pierre-Francois Van De Moortele
Original Assignee
Regents Of The University Of Minnesota
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Regents Of The University Of Minnesota filed Critical Regents Of The University Of Minnesota
Priority to US15/513,042 priority Critical patent/US10545210B2/en
Publication of WO2016049013A1 publication Critical patent/WO2016049013A1/en

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00Arrangements or instruments for measuring magnetic variables
    • G01R33/20Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48NMR imaging systems
    • G01R33/54Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/56Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
    • G01R33/561Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution by reduction of the scanning time, i.e. fast acquiring systems, e.g. using echo-planar pulse sequences
    • G01R33/5611Parallel magnetic resonance imaging, e.g. sensitivity encoding [SENSE], simultaneous acquisition of spatial harmonics [SMASH], unaliasing by Fourier encoding of the overlaps using the temporal dimension [UNFOLD], k-t-broad-use linear acquisition speed-up technique [k-t-BLAST], k-t-SENSE
    • G01R33/5612Parallel RF transmission, i.e. RF pulse transmission using a plurality of independent transmission channels
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00Arrangements or instruments for measuring magnetic variables
    • G01R33/20Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/28Details of apparatus provided for in groups G01R33/44 - G01R33/64
    • G01R33/288Provisions within MR facilities for enhancing safety during MR, e.g. reduction of the specific absorption rate [SAR], detection of ferromagnetic objects in the scanner room
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00Arrangements or instruments for measuring magnetic variables
    • G01R33/20Arrangements or instruments for measuring magnetic variables involving magnetic resonance
    • G01R33/44Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
    • G01R33/48NMR imaging systems
    • G01R33/54Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
    • G01R33/543Control of the operation of the MR system, e.g. setting of acquisition parameters prior to or during MR data acquisition, dynamic shimming, use of one or more scout images for scan plane prescription

Definitions

  • the field of the invention is systems and methods for magnetic resonance imaging ("MRI" ⁇ . More particularly, the invention relates to systems and methods for designing radio frequency (“RF" ⁇ pulses for use in parallel transmission MRI.
  • MRI magnetic resonance imaging
  • RF radio frequency
  • the present invention overcomes the aforementioned drawbacks by providing a method for designing parallel transmission radio frequency (“RF" ⁇ pulses for use in magnetic resonance imaging (“MRI” ⁇ .
  • the method includes providing a set of pre-computed temperature matrices to a computer system, where the temperature matrices describe temperature rises in response to RF pulses based on a linear thermal model.
  • a set of temperature virtual observation points (“TVOPs" ⁇ based on the pre- computed temperature matrices is then selected, and a plurality of RF pulse waveforms for use in MRI are computed by solving an optimization problem that is constrained by a maximum temperature rise evaluated at the set of TVOPs.
  • FIG. 1 is a flowchart setting forth the steps of an example method for designing radio frequency ("RF" ⁇ pulse waveforms subject to a temperature rise constraint.
  • FIG. 2 is a flowchart setting forth the steps of an example method for setting a temperature rise constraint based on determining a set of temperature virtual observation points ("TVOPs" ⁇ at which temperature rise should be evaluated.
  • TVOPs temperature virtual observation points
  • FIG. 3A illustrates an example of specific absorption rate ("SAR" ⁇ virtual observation points, shown as white squares in two different sagittal slices.
  • FIG. 3B illustrates an example of TVOPs, shown as white squares in the same two sagittal slices as depicted in FIG. 3A.
  • FIG. 4 illustrates a validation of TVOPs for a six minute sequence. Twenty thousand random RF waveforms were used to compute a SAR distribution, which yielded a maximum temperature rise over a full model, a maximum temperature over the TVOPs, and a maximum temperature over the TVOPs plus a defined tolerance. The results are scaled for 10 W total RF power.
  • FIG. 5 depicts histograms of the maximum temperature rises for different sequence durations, calculated over 20,000 random scenarios.
  • the waveforms are scaled to always yield 10 W/kg peak 10-g SAR.
  • the global SAR limit of 3.2 W/kg was never reached.
  • FIG. 6 depicts flip angle distributions for different pulse designs, which from top to bottom include SAR-constraint pulse design and 1°C, 1.5°C, 2°C temperature-constraints pulse designs.
  • FIG.7 depicts resulting temperature rises for different pulse designs, which from top to bottom include SAR-constraint pulse design and 1°C, 1.5°C, 2°C temperature-constraints pulse designs.
  • FIG. 8 is a block diagram of an example of a magnetic resonance imaging
  • Described here are systems and methods for utilizing parallel transmission in magnetic resonance imaging ("MRI" ⁇ to mitigate B1+ inhomogeneities while explicitly constraining the temperature rise during parallel transmission.
  • MRI magnetic resonance imaging
  • finite difference time domain simulations can be performed on a numerical human head model and for a 16 channel coil at 10.5 Tesla.
  • a virtual observation point compression model for the temperature rise is derived. This compact representation is then used in a nonlinear programming algorithm for pulse design under explicit temperature rise constraints.
  • radio frequency (“RF" ⁇ pulse performance in some cases is increased by a factor of two compared to specific absorption rate (“SAR" ⁇ constrained RF pulses, while temperature rise is directly and efficiently controlled.
  • RF pulse performance can be gained by relaxing the SAR constraints, but at the expense of a loss of direct control on temperature.
  • the systems and methods of the present invention design RF pulses using a technique that explicitly controls temperature rise via a compression model that is based on virtual observation points ("VOPs" ⁇ .
  • VOPs virtual observation points
  • Thermal pre-calculations are first carried out for a given RF exposure time, coil, and subject model in order to obtain complex temperature matrices, T , after which the compression scheme is implemented.
  • the thermal model employed can be Pennes' bio-heat equation.
  • a nonlinear programming algorithm can be implemented for designing the RF pulses.
  • the algorithm used for RF pulse design is based on the active-set ("A-S" ⁇ method, which has been previously shown to be fast, robust, and powerful for nonlinear programming problems.
  • the method includes selecting at least one region to be imaged within a subject, as indicated at step 102.
  • previously acquired images of the subject can be used for determining or identifying the regions to be imaged, and for which the RF pulses should be designed.
  • an objective function is then initialized using the selected regions.
  • this step can include providing an encoding matrix for each region, and a desired flip angle for each region.
  • the objective function can be based on a ma nitude least squares ("MLS" ⁇ problem having the following form:
  • A is a diagonal matrix whose diagonal entries, A n , correspond to the encoding matrix for each region in units of flip angle per square root of power
  • X is a vector whose entries, X n , correspond to the RF pulses waveforms designed to be applied to the n' h region
  • is a vector whose entries, ⁇ ⁇ , correspond to the flip angle for the RF pulses designed to be applied to the n' h region.
  • Eqn. (1 ⁇ can have the following form:
  • the matrices A j and A 2 have the units of flip angle per square root of power (W ⁇ , and encode the spins' dynamics in the small tip angle approximation for the first and second regions, respectively, of the excitation volume-of-interest.
  • W ⁇ flip angle per square root of power
  • 0 X can equal 20 degrees and ⁇ 2 can equal 60 degrees.
  • a set of RF waveforms is then determined by optimizing the objective function, as indicated at step 106.
  • the optimization problem can be solved using an appropriate nonlinear programming algorithm.
  • the nonlinear programming algorithm can include an active-set ("A-S" ⁇ algorithm. It will be appreciated, however, that other methods can also be used in the pulse design problem. For example, interior point methods could also be used.
  • the A-S algorithm When using the A-S algorithm, it can be implemented via the optimization toolbox of Matlab (The Mathworks, Natick, MA, USA ⁇ .
  • the A-S algorithm is an iterative (Quasi-Newton] procedure where, at each step, the problem is approximated by a quadratic program with the constraints linearized.
  • the number of iterations i.e., the number of times the quadratic problem was solved] for this implementation of the temperature constraint pulse design can be set to a predetermined number of iterations, such as 300; however, it will be appreciated that other stopping criterion could also be suitably implemented.
  • each sub-pulse can be a one millisecond sine pulse apodized with a Hanning window and with a time- bandwidth product equal to four. Further, in this example the duty-cycle can be taken as 20 percent.
  • the optimization problem is constrained using a temperature rise constraint, as indicated at step 108.
  • the temperature rise constraint can be implemented using compressed temperature matrices, thereby reducing the computational complexity of the RF design process.
  • the compressed temperature matrices can be evaluated at a set of temperature virtual observation points ("TVOPs" ⁇ .
  • the temperature rise constraint can be,
  • N TVOPs is the number of TVOPs obtained after compression and c. (x) are the values over the temperature rises calculated over the TVOPs, with the predefined tolerance.
  • the temperature rise constraint is described below in more detail.
  • other constraints can be implemented, such as peak power constraint, as indicated at step 110.
  • a peak power constraint can be defined as,
  • the optimization of the objective function can also take into consideration other factors, including RF frequency, average power, number of RF channels, and so on.
  • the designed RF pulse waveforms can be provided to an MRI system, as indicated at step 112.
  • the MRI system can then be programmed or otherwise directed to generate RF pulses, preferably using parallel transmission, using the designed RF waveforms in order to achieve a desired effect in the selected regions. Because of the constraint on temperature rise, the applied RF pulses will not cause a significant temperature rise in the subject.
  • FIG. 2 a flowchart is illustrated as setting forth the steps of an example method for setting a temperature rise constraint based on TVOPs.
  • the method includes pre-computing a number of temperature matrices and then reducing the set of temperature matrices to a compressed set of TVOPs.
  • the method thus includes pre-computing a plurality of temperature matrices, as indicated at step 202.
  • the temperature matrices described here model a temperature rise, but could also be readily adapted to model absolute temperature.
  • One method for pre-computing these temperature matrices is now described.
  • T is the local temperature (K)
  • f is a linear function of temperature
  • g is an RF source-independent function leading to a nonzero equilibrium temperature in the absence of RF. All terms may be spatially dependent. In the case of Pennes' bio-heat model, the different terms are,
  • variable, g which contains T b and C, drops out of the equation.
  • the temperature rise does not depend on the equilibrium temperature or on the possibly uncertain function, g .
  • SAR can be calculated by using Hermitian local matrices, Q( ) , defined for every voxel as,
  • E x , E , and E z are normalized complex row vectors containing electric field values coming from each channel, and for the corresponding axis, at the location, r .
  • I C being the
  • T r are the temperature matrices that can be pre- computed.
  • Ar(r) x ⁇ T(r)x (13 ⁇ .
  • pre-computations as described above, the computation of the temperature rises can be sped up for fast optimization.
  • pre- computations can be performed using a finite difference time domain ("FDTD" ⁇ method.
  • FDTD-based electromagnetic simulations were used to perform simulation of the B 1 and electric fields. Once computed, the field maps were eventually downsampled to a 4x 4x 4 mm 3 grid to make the SAR and thermal calculations more tractable. After computing the Q matrices for each voxel, Q 10g matrices were constructed to calculate
  • the method for setting a temperature rise constraint next selects an initial set of TVOPs based on the computed temperature matrices, as indicated at step 204. For instance, once the temperature matrices, T(r) , have been computed (e.g., built from the AT. . elements], the number of constraints can be compressed based, in part, on the knowledge that the temperature rise for an arbitrary pulse can be compactly calculated. In particular, a TVOP compression scheme can be implemented.
  • the TVOPs can thus be defined as the matrices among the T matrices obtained via the pre-computations that provide a lower-bound and upper-bound of the maximum temperature rise obtained over the whole model.
  • the T matrices can be computed in a number of different ways.
  • the TVOPs can be computed based on measurements of and B 0 and using model-based E -field values to calculate E- fields. This approach is non-specific to the subject being imaged.
  • the TVOPs can be computed based on measurements of and B 0 and using a E -field algorithm to generate E -field measurements specific to the subject being imaged.
  • the E -field algorithm can be the one described in copending US Patent Application Serial No. 14/245,145, which is herein incorporated by reference in its entirety.
  • the TVOPs can be defined as the set of T matrices, which guarantees that,
  • AT k denotes the temperature rise calculated for the k th TVOP. In general, larger values of ⁇ will result in more compression (i.e., fewer TVOPs ⁇ .
  • the identity matrix can also be used. In this case, the tolerance is expressed in terms of the average RF power.
  • TVOPs will also depend on the duration of the RF exposure. By construction of the TVOPs, it is thereby ensured that the maximum temperature rise over the full model is always upper-bounded by one of the TVOPs plus the defined tolerance, and lower bounded by the maximum over the TVOPs.
  • the T matrices are first sorted according to their maximum eigenvalues in decreasing order.
  • the set of TVOPs is then initialized as the first T matrix in the sorted list.
  • the initial set of TVOPs can be further refined, as indicated at step 206.
  • the set of TVOPs can be reduced based on a search for a set of N positive coefficients, c n , where N is the current number of TVOPs.
  • the set of coefficients that is searched for has the sum,
  • a candidate TVOP, T is added to the initial set of TVOPs unless a set of coefficients that renders the following metric to be nonnegative can be found for that candidate TVOP:
  • T is the candidate TVOP and T G is the global temperature matrix
  • the temperature constraint is defined, as indicated at step 208.
  • T max corresponds to the maximum allowable rise in temperature, and can be selected as any suitable value that constrains the temperature rise to safe values.
  • T max can be in the range of about 1.0-2.0 degrees Celsius.
  • FIG. 3A An example set of SAR virtual observation points ("VOPs" ⁇ computed for a subject is illustrated in FIG. 3A, and an example set of TVOPs determined according to the methods described here for the same subject is illustrated in FIG. 3B.
  • VOPs virtual observation points
  • FIG. 3B An example set of SAR virtual observation points
  • the locations of the TVOPs shown as white squares in FIG. 3B, are not correlated with the locations of the SAR VOPs, shown as white squares in FIG. 3A.
  • the SAR VOPs can be located deep inside the brain, the TVOPs are mostly located in low-perfused regions (e.g., skin, muscle, cerebrospinal fluid ⁇ .
  • FIG. 4 illustrates a validation of an example set of TVOPs determined for a six minute sequence, with the temperature rise calculated over the TVOPs, the same temperature plus the defined tolerance, and the temperature calculated over the whole model.
  • the temperature calculated over the whole model (shown as a dashed line] is upper-bounded by the TVOPs plus the selected tolerance, and is lower-bounded by the TVOPs.
  • the calculation was performed for 20,000 random RF waveforms having a uniform distribution for amplitude and phase.
  • the total power for each waveform was scaled to 10 W.
  • FIG. 5 depicts histograms of the maximum temperature rises for six, ten, and thirty minute RF exposures, computed with the same 20,000 random RF scenarios.
  • FIG. 5 depicts histograms of the maximum temperature rises for six, ten, and thirty minute RF exposures, computed with the same 20,000 random RF scenarios.
  • Table 1 summarizes the results of an example pulse design using the A-S algorithm under various constraints, for six, ten, and thirty minutes time-of-flight sequences.
  • Constraint NRMSE (%) SAR lOg (W/kg) SAR Global (W/kg) Temperature (°C)
  • the A-S algorithm converged to a local minimum after 300 iterations in about 30 seconds for the SAR and temperature-constrained pulse designs.
  • the peak power constraint (2 kW in this example] was always saturated.
  • the pulse design under explicit temperature rise constraints noticeably outperformed the pulse design constrained under SAR, with a lower Normalized Root Mean Square Error ("NRMSE" ⁇ .
  • FIGS. 6 and 7 Examples of flip angle maps and temperature rise maps determined for the 6 minute sequence are illustrated in FIGS. 6 and 7, respectively. While the target flip angles in FIG. 6 cannot be achieved in the SAR-constraint pulse design case, the correct flip angles can be obtained in the temperature rise constrained case as the temperature rise can be progressively increased in a controlled manner. The corresponding temperature rises are shown in FIG. 7, where the hot spots mainly occur in the skin and CSF regions.
  • the TVOP approach described here may require more computing efforts because a different set of TVOPs exists for each pulse sequence duration.
  • the TVOPs can be pre-computed only once.
  • the locations of the TVOPs also indicate that the maximum temperature rises occur in less sensitive regions (e.g., in the skin] than the brain or other internal structures, where the equilibrium temperature (i.e. the steady state obtained without RF source] is less than 37°C. As such, it is possible that a 1.0-2.0°C temperature rise constraint can be implemented while remaining conservative with respect to the absolute temperature guidelines.
  • the proposed pulse design method also presents a natural increase in complexity compared to SAR-constrained optimizations because the sliding-integration window used for SAR calculations can no longer be applied, which is especially problematic when different sequences using different RF pulses follow each other.
  • a pause can be interleaved between different pulse sequences to reset the temperature, or the whole imaging exam can be taken into account.
  • the MRI system 800 includes an operator workstation 802, which will typically include a display 804; one or more input devices 806, such as a keyboard and mouse; and a processor 808.
  • the processor 808 may include a commercially available programmable machine running a commercially available operating system.
  • the operator workstation 802 provides the operator interface that enables scan prescriptions to be entered into the MRI system 800.
  • the operator workstation 802 may be coupled to four servers: a pulse sequence server 810; a data acquisition server 812; a data processing server 814; and a data store server 816.
  • the operator workstation 802 and each server 810, 812, 814, and 816 are connected to communicate with each other.
  • the servers 810, 812, 814, and 816 may be connected via a communication system 840, which may include any suitable network connection, whether wired, wireless, or a combination of both.
  • the communication system 840 may include both proprietary or dedicated networks, as well as open networks, such as the internet.
  • the pulse sequence server 810 functions in response to instructions downloaded from the operator workstation 802 to operate a gradient system 818 and a radiofrequency ("RF" ⁇ system 820. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system 818, which excites
  • the gradient coil assembly 822 forms part of a magnet assembly 824 that includes a polarizing magnet 826 and a whole-body RF coil 828.
  • RF waveforms are applied by the RF system 820 to the RF coil 828, or a separate local coil (not shown in FIG. 8], in order to perform the prescribed magnetic resonance pulse sequence.
  • Responsive magnetic resonance signals detected by the RF coil 828, or a separate local coil (not shown in FIG. 8] are received by the RF system 820, where they are amplified, demodulated, filtered, and digitized under direction of commands produced by the pulse sequence server 810.
  • the RF system 820 includes an RF transmitter for producing a wide variety of RF pulses used in MRI pulse sequences.
  • the RF transmitter is responsive to the scan prescription and direction from the pulse sequence server 810 to produce RF pulses of the desired frequency, phase, and pulse amplitude waveform.
  • the generated RF pulses may be applied to the whole-body RF coil 828 or to one or more local coils or coil arrays (not shown in FIG. 8 ⁇ .
  • the RF system 820 also includes one or more RF receiver channels. Each
  • RF receiver channel includes an RF preamplifier that amplifies the magnetic resonance signal received by the coil 828 to which it is connected, and a detector that detects and digitizes the I and ⁇ quadrature components of the received magnetic resonance signal.
  • the magnitude of the received magnetic resonance signal may, therefore, be determined at any sampled point by the square root of the sum of the squares of the I and ⁇ components: (18);
  • phase of the received magnetic resonance signal may also be determined according to the following relationship:
  • the pulse sequence server 810 also optionally receives patient data from a physiological acquisition controller 830.
  • the physiological acquisition controller 830 may receive signals from a number of different sensors connected to the patient, such as electrocardiograph ("ECG") signals from electrodes, or respiratory signals from a respiratory bellows or other respiratory monitoring device.
  • ECG electrocardiograph
  • Such signals are typically used by the pulse sequence server 810 to synchronize, or "gate,” the performance of the scan with the subject's heart beat or respiration.
  • the pulse sequence server 810 also connects to a scan room interface circuit 832 that receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit 832 that a patient positioning system 834 receives commands to move the patient to desired positions during the scan.
  • the digitized magnetic resonance signal samples produced by the RF system 820 are received by the data acquisition server 812.
  • the data acquisition server 812 operates in response to instructions downloaded from the operator workstation 802 to receive the real-time magnetic resonance data and provide buffer storage, such that no data is lost by data overrun. In some scans, the data acquisition server 812 does little more than pass the acquired magnetic resonance data to the data processor server 814. However, in scans that require information derived from acquired magnetic resonance data to control the further performance of the scan, the data acquisition server 812 is programmed to produce such information and convey it to the pulse sequence server 810. For example, during prescans, magnetic resonance data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 810.
  • navigator signals may be acquired and used to adjust the operating parameters of the RF system 820 or the gradient system 818, or to control the view order in which k-space is sampled.
  • the data acquisition server 812 may also be employed to process magnetic resonance signals used to detect the arrival of a contrast agent in a magnetic resonance angiography ("MRA" ⁇ scan.
  • MRA magnetic resonance angiography
  • the data acquisition server 812 acquires magnetic resonance data and processes it in real-time to produce information that is used to control the scan.
  • the data processing server 814 receives magnetic resonance data from the data acquisition server 812 and processes it in accordance with instructions downloaded from the operator workstation 802. Such processing may, for example, include one or more of the following: reconstructing two-dimensional or three- dimensional images by performing a Fourier transformation of raw k-space data; performing other image reconstruction algorithms, such as iterative or backprojection reconstruction algorithms; applying filters to raw k-space data or to reconstructed images; generating functional magnetic resonance images; calculating motion or flow images; and so on.
  • Images reconstructed by the data processing server 814 are conveyed back to the operator workstation 802 where they are stored.
  • Real-time images are stored in a data base memory cache (not shown in FIG. 8], from which they may be output to operator display 802 or a display 836 that is located near the magnet assembly 824 for use by attending physicians.
  • Batch mode images or selected real time images are stored in a host database on disc storage 838.
  • the data processing server 814 notifies the data store server 816 on the operator workstation 802.
  • the operator workstation 802 may be used by an operator to archive the images, produce films, or send the images via a network to other facilities.
  • the MRI system 800 may also include one or more networked workstations 842.
  • a networked workstation 842 may include a display 844; one or more input devices 846, such as a keyboard and mouse; and a processor 848.
  • the networked workstation 842 may be located within the same facility as the operator workstation 802, or in a different facility, such as a different healthcare institution or clinic.
  • the networked workstation 842 may gain remote access to the data processing server 814 or data store server 816 via the communication system 840. Accordingly, multiple networked workstations 842 may have access to the data processing server 814 and the data store server 816. In this manner, magnetic resonance data, reconstructed images, or other data may be exchanged between the data processing server 814 or the data store server 816 and the networked workstations 842, such that the data or images may be remotely processed by a networked workstation 842. This data may be exchanged in any suitable format, such as in accordance with the transmission control protocol ("TCP" ⁇ , the internet protocol (“IP" ⁇ , or other known or suitable protocols.
  • TCP transmission control protocol
  • IP internet protocol

Abstract

Systems and methods for designing RF pulses using a technique that directly controls temperature rise via a compression model that is based on virtual observation points ("VOPs"} are provided. Thermal pre-simulations are first carried out for a given RF exposure time, coil, and subject model in order to obtain complex temperature matrices, after which the compression scheme follows. As one example, the thermal model employed can be Pennes' bio-heat equation. Focusing design constraints on the temperature rise instead of the absolute temperature allows for uncertain parameters to be dropped from the thermal model, making it more robust and less prone to errors. In some embodiments, the algorithm used for RF pulse design is the active-set ("A-S"} method.

Description

LOCAL TEMPERATURE RISE CONSTRAINED RADIO FREQUENCY PULSE DESIGN IN PARALLEL TRANSMISSION MRI
CROSS-REFERENCE TO RELATED APPLICATIONS
[0001] This application claims the benefit of U.S. Provisional Patent Application
Serial No. 62/053,603, filed on September 22, 2014, and entitled "LOCAL TEMPERATURE RISE CONSTRAINED RADIO FREQUENCY PULSE DESIGN IN PARALLEL TRANSMISSION."
STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH
[0002] This invention was made with government support under EB 015894 awarded by the National Institutes of Health. The government has certain rights in the invention.
BACKGROUND OF THE INVENTION
[0003] The field of the invention is systems and methods for magnetic resonance imaging ("MRI"}. More particularly, the invention relates to systems and methods for designing radio frequency ("RF"} pulses for use in parallel transmission MRI.
[0004] Although there is a general consensus that temperature is the true relevant radio frequency ("RF"} related safety parameter, tracking specific absorption rate ("SAR"} in MR exams and in RF pulse design has remained the gold standard. The fact that SAR has remained the gold standard can probably be explained by the complexity of bio-heat equation models and the lack of fast temperature mapping methods.
[0005] Given the multi-factorial dependence of temperature, SAR computation has been a great simplification for handling safety, which has been supported by experimental evidence. Whereas RF coil technology, static field intensity, and temperature guidelines have evolved throughout the years, SAR thresholds have essentially remained identical.
[0006] In a recent study by A. Massire, et al., "Thermal simulations in the human head for high field MRI using parallel transmission," J Magn Reson Imaging 2012; 35:1312-1321, a series of thermal simulations was conducted on a numerical human head model containing 20 anatomical structures in order to verify the consistency between the SAR and temperature guidelines for parallel transmission (pTx} RF exposures at 7 T. Thousands of different RF scenarios were simulated. In these models, starting from the equilibrium temperature (i.e. the steady-state obtained when there is no RF heat source], the temperature almost never exceeded 38 °C over 30 minutes of RF pulsing as long as 10-g average SAR did not exceed its limit of 10 W/kg.
[0007] These simulations, however, relied on the validity of Pennes' bioheat model whose constant blood temperature assumption is currently debated. Motivated also by the lack of direct correspondence between temperature and SAR, the RF pulse design for a time-of-flight (TOF] sequence under strict absolute temperature constraints was investigated numerically, whereby strict SAR constraints were adjusted in a feedback-based manner to fulfill the temperature constraints. Although the IEC guidelines in general involve upper limits for absolute temperature, the latter actually appears quite dependent on daily fluctuations caused by uncontrolled events occurring in the patient's life (fever, physical activity, stress, weather etc}. Furthermore, in vivo absolute temperature measurement with MR is still extremely challenging.
SUMMARY OF THE INVENTION
[0008] The present invention overcomes the aforementioned drawbacks by providing a method for designing parallel transmission radio frequency ("RF"} pulses for use in magnetic resonance imaging ("MRI"}. The method includes providing a set of pre-computed temperature matrices to a computer system, where the temperature matrices describe temperature rises in response to RF pulses based on a linear thermal model. A set of temperature virtual observation points ("TVOPs"} based on the pre- computed temperature matrices is then selected, and a plurality of RF pulse waveforms for use in MRI are computed by solving an optimization problem that is constrained by a maximum temperature rise evaluated at the set of TVOPs.
[0009] The foregoing and other aspects and advantages of the invention will appear from the following description. In the description, reference is made to the accompanying drawings that form a part hereof, and in which there is shown by way of illustration a preferred embodiment of the invention. Such embodiment does not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims and herein for interpreting the scope of the invention.
BRIEF DESCRIPTION OF THE DRAWINGS [0010] FIG. 1 is a flowchart setting forth the steps of an example method for designing radio frequency ("RF"} pulse waveforms subject to a temperature rise constraint.
[0011] FIG. 2 is a flowchart setting forth the steps of an example method for setting a temperature rise constraint based on determining a set of temperature virtual observation points ("TVOPs"} at which temperature rise should be evaluated.
[0012] FIG. 3A illustrates an example of specific absorption rate ("SAR"} virtual observation points, shown as white squares in two different sagittal slices.
[0013] FIG. 3B illustrates an example of TVOPs, shown as white squares in the same two sagittal slices as depicted in FIG. 3A.
[0014] FIG. 4 illustrates a validation of TVOPs for a six minute sequence. Twenty thousand random RF waveforms were used to compute a SAR distribution, which yielded a maximum temperature rise over a full model, a maximum temperature over the TVOPs, and a maximum temperature over the TVOPs plus a defined tolerance. The results are scaled for 10 W total RF power.
[0015] FIG. 5 depicts histograms of the maximum temperature rises for different sequence durations, calculated over 20,000 random scenarios. The waveforms are scaled to always yield 10 W/kg peak 10-g SAR. The global SAR limit of 3.2 W/kg was never reached.
[0016] FIG. 6 depicts flip angle distributions for different pulse designs, which from top to bottom include SAR-constraint pulse design and 1°C, 1.5°C, 2°C temperature-constraints pulse designs.
[0017] FIG.7 depicts resulting temperature rises for different pulse designs, which from top to bottom include SAR-constraint pulse design and 1°C, 1.5°C, 2°C temperature-constraints pulse designs.
[0018] FIG. 8 is a block diagram of an example of a magnetic resonance imaging
("MRI"} system.
DETAILED DESCRIPTION OF THE INVENTION
[0019] Described here are systems and methods for utilizing parallel transmission in magnetic resonance imaging ("MRI"} to mitigate B1+ inhomogeneities while explicitly constraining the temperature rise during parallel transmission. By way of example, finite difference time domain simulations can be performed on a numerical human head model and for a 16 channel coil at 10.5 Tesla. Based on a set of pre- computations, a virtual observation point compression model for the temperature rise is derived. This compact representation is then used in a nonlinear programming algorithm for pulse design under explicit temperature rise constraints.
[0020] In the example of a time-of-flight sequence, radio frequency ("RF"} pulse performance in some cases is increased by a factor of two compared to specific absorption rate ("SAR"} constrained RF pulses, while temperature rise is directly and efficiently controlled. RF pulse performance can be gained by relaxing the SAR constraints, but at the expense of a loss of direct control on temperature.
[0021] There remains a need for accurate safety control at ultra-high magnetic field strengths, especially given the lack of correspondence between SAR and temperature. It is contemplated that the systems and methods described here will provide for safer and more efficient magnetic resonance imaging at these ultra-high field strengths, as well as at more routine clinical magnetic field strengths.
[0022] As will be described below in more detail, the systems and methods of the present invention design RF pulses using a technique that explicitly controls temperature rise via a compression model that is based on virtual observation points ("VOPs"}. Thermal pre-calculations are first carried out for a given RF exposure time, coil, and subject model in order to obtain complex temperature matrices, T , after which the compression scheme is implemented. As one example, the thermal model employed can be Pennes' bio-heat equation. As will be shown, focusing on the temperature rise instead of the absolute temperature allows for uncertain parameters to be dropped from the thermal model, making it more robust and less prone to errors. A nonlinear programming algorithm can be implemented for designing the RF pulses. In some embodiments, the algorithm used for RF pulse design is based on the active-set ("A-S"} method, which has been previously shown to be fast, robust, and powerful for nonlinear programming problems.
[0023] Referring now to FIG. 1, a flowchart is illustrated as setting forth the steps of an example method for designing one or more RF pulses for use in MRI, in which the RF pulses are optimized to minimize B1+ inhomogeneities while minimizing the temperature rise that will be caused in a subject when the designed RF pulses are applied to that subject. [0024] The method includes selecting at least one region to be imaged within a subject, as indicated at step 102. In some implementations, previously acquired images of the subject can be used for determining or identifying the regions to be imaged, and for which the RF pulses should be designed.
[0025] As indicated at step 104, an objective function is then initialized using the selected regions. In general, this step can include providing an encoding matrix for each region, and a desired flip angle for each region. As one example, the objective function can be based on a ma nitude least squares ("MLS"} problem having the following form:
Figure imgf000006_0001
[0026] where A is a diagonal matrix whose diagonal entries, An , correspond to the encoding matrix for each region in units of flip angle per square root of power, X is a vector whose entries, Xn , correspond to the RF pulses waveforms designed to be applied to the n'h region; and Θ is a vector whose entries, θη , correspond to the flip angle for the RF pulses designed to be applied to the n'h region.
[0027] As a specific example in which two different regions are selected for imaging, Eqn. (1} can have the following form:
Figure imgf000006_0002
[0028] where the matrices Aj and A2 have the units of flip angle per square root of power (W}, and encode the spins' dynamics in the small tip angle approximation for the first and second regions, respectively, of the excitation volume-of-interest. As an example, 0X can equal 20 degrees and θ2 can equal 60 degrees.
[0029] A set of RF waveforms is then determined by optimizing the objective function, as indicated at step 106. The optimization problem can be solved using an appropriate nonlinear programming algorithm. As one example, the nonlinear programming algorithm can include an active-set ("A-S"} algorithm. It will be appreciated, however, that other methods can also be used in the pulse design problem. For example, interior point methods could also be used.
[0030] When using the A-S algorithm, it can be implemented via the optimization toolbox of Matlab (The Mathworks, Natick, MA, USA}. The A-S algorithm is an iterative (Quasi-Newton] procedure where, at each step, the problem is approximated by a quadratic program with the constraints linearized. The number of iterations (i.e., the number of times the quadratic problem was solved] for this implementation of the temperature constraint pulse design can be set to a predetermined number of iterations, such as 300; however, it will be appreciated that other stopping criterion could also be suitably implemented.
[0031] Based on the example objective function in Eqn. (2], as one example for the RF pulse design a kr-spokes technique can be used so that Xj and x2 are, respectively, a 2 kr-spoke pulse and 3 kr-spoke pulse. In this example, each sub-pulse can be a one millisecond sine pulse apodized with a Hanning window and with a time- bandwidth product equal to four. Further, in this example the duty-cycle can be taken as 20 percent.
[0032] Referring still to FIG. 1, and as mentioned above, the optimization problem is constrained using a temperature rise constraint, as indicated at step 108. In general, the temperature rise constraint can be implemented using compressed temperature matrices, thereby reducing the computational complexity of the RF design process. As one example, the compressed temperature matrices can be evaluated at a set of temperature virtual observation points ("TVOPs"}. As one example, the temperature rise constraint can be,
Figure imgf000007_0001
[0033] where NTVOPs is the number of TVOPs obtained after compression and c. (x) are the values over the temperature rises calculated over the TVOPs, with the predefined tolerance. The temperature rise constraint is described below in more detail. In addition, other constraints can be implemented, such as peak power constraint, as indicated at step 110. As one example, a peak power constraint can be defined as,
cpJc (x) = |¾|2 < 2kW,k = 1,...,NcNkT (4).
[0034] The optimization of the objective function can also take into consideration other factors, including RF frequency, average power, number of RF channels, and so on.
[0035] The designed RF pulse waveforms can be provided to an MRI system, as indicated at step 112. The MRI system can then be programmed or otherwise directed to generate RF pulses, preferably using parallel transmission, using the designed RF waveforms in order to achieve a desired effect in the selected regions. Because of the constraint on temperature rise, the applied RF pulses will not cause a significant temperature rise in the subject.
[0036] Referring now to FIG. 2, a flowchart is illustrated as setting forth the steps of an example method for setting a temperature rise constraint based on TVOPs. In general, the method includes pre-computing a number of temperature matrices and then reducing the set of temperature matrices to a compressed set of TVOPs.
[0037] The method thus includes pre-computing a plurality of temperature matrices, as indicated at step 202. In particular, the temperature matrices described here model a temperature rise, but could also be readily adapted to model absolute temperature. One method for pre-computing these temperature matrices is now described.
[0038] A general linear absolute temperature model may be written as, pcp ^ = f (T) + g + pSAR (5);
[0039] where p is the tissue density (kg/m3}, cp is the specific heat capacity
Q/kg/K], T is the local temperature (K], f is a linear function of temperature, and g is an RF source-independent function leading to a nonzero equilibrium temperature in the absence of RF. All terms may be spatially dependent. In the case of Pennes' bio-heat model, the different terms are,
f (T(r)) = V - (k(r)VT(r)) - B (r)T(r) (6); [0040] where g (r) = B (r^Tb (r) + C(r) ; k is the thermal conductivity
(W/K/m}; B is the perfusion constant (W/K/m3}; Tb is the constant blood temperature
(K}; and C is the metabolic rate (W/m3}. As discussed above, modeling based on the absolute temperature introduces uncertainties into the RF pulse design process because both [Τ[τ^ and g (r) must be known. Some uncertainties can be mitigated, however, by modeling based on the temperature rise instead. To restructure the thermal model accordingly, the absolute temperature, T , can first be set to be equal to the equilibrium temperature, T , (i.e. the steady state obtained when there is no RF heat source] plus the temperature rise Δ7" . This modification to the thermal model in Eqn. (5} yields the following: pcp ^ = f (AT) + pSAR (7);
[0041] In making this modification, the variable, g , which contains Tb and C, drops out of the equation. In this instance, the temperature rise does not depend on the equilibrium temperature or on the possibly uncertain function, g .
[0042] SAR can be calculated by using Hermitian local matrices, Q( ) , defined for every voxel as,
Q (r) = ^(E + E + ΕΛ' ) C8);
[0043] where Ex , E , and Ez are normalized complex row vectors containing electric field values coming from each channel, and for the corresponding axis, at the location, r . Given NT -point time waveforms dispatched on Nc channels, IC being the
Nr -by- Ne array containing all these waveforms, for 100 percent duty cycle the SAR at voxel, r , is given by:
Figure imgf000009_0001
[0044] where / and j respectively refer to the row and column indices of the corresponding matrices. Because Eqn. (7} is linear with respect to the RF heat source, it can be numerically solved for each map to obtain, for a given duration, a corresponding Δ7 . (r) complex distribution, which is a solution of the following: dAT. . (r) , .
PCP dt = f (ATu + PQ J (r) (10)-
[0045] Because Q; and Qj ; {f ) are complex conjugates, by construction so are Δ7^. (Γ ) and ΔΓ. . (r) . After Nc (Nc + 1)/2 computations, which may be pre- computations, and by exploiting the linearity of the model, the temperature rise, Δ7" (Γ) , at any point can be calculated using complex matrices, T(r) , which can be defined as: A (r)
T(r) = (ii);
ATN r) - · T,
[0046] so that for any waveform array, / , the following is obtained:
1
(12).
T J
[0047] These matrices, T r) , are the temperature matrices that can be pre- computed. For a static RF shim with the x-vector containing the Nc complex weights, the previous result becomes,
Ar(r) = xT(r)x (13}.
[0048] Using pre-computations as described above, the computation of the temperature rises can be sped up for fast optimization. As one example, pre- computations can be performed using a finite difference time domain ("FDTD"} method.
[0049] In one non-limiting example of implementing the methods described above, FDTD-based electromagnetic simulations were used to perform simulation of the B1 and electric fields. Once computed, the field maps were eventually downsampled to a 4x 4x 4 mm3 grid to make the SAR and thermal calculations more tractable. After computing the Q matrices for each voxel, Q10g matrices were constructed to calculate
10-g SAR values by using a region-growing algorithm. A global SAR QG matrix was obtained by averaging the energy density over the relevant mass. With all these Q matrices available, the T matrices and TVOPs could be computed for each sequence duration.
[0050] While SAR averaging guidelines (10-g or global] dilute the spurious electric field spikes due to stair-casing errors in FDTD simulations, temperature calculations arising from voxel-wise SAR distributions are more sensitive to numerical artifacts and can lead to unrealistic temperature rises, even when taking into account heat diffusion effects. These electric field spikes may lead to special SAR averaging rules at the skin interface, where large discontinuities in conductivity typically occur. Due to the speckle noise nature of these spurious spikes, the electric field maps can be suitably filtered to remove the spikes. As one example, a 3x3x 3 median filter can be applied to the electric field distributions for each channel, due to the efficiency of this filter to remove this particular kind of noise. The filtered electric field distributions can be used both for SAR and temperature calculations.
[0051] Referring again to FIG. 2, the method for setting a temperature rise constraint next selects an initial set of TVOPs based on the computed temperature matrices, as indicated at step 204. For instance, once the temperature matrices, T(r) , have been computed (e.g., built from the AT. . elements], the number of constraints can be compressed based, in part, on the knowledge that the temperature rise for an arbitrary pulse can be compactly calculated. In particular, a TVOP compression scheme can be implemented.
[0052] As one example, assuming a discretized model of the volume exposed to the radio frequency field, there can be as many constraints as there are voxels (i.e., on the order of tens of thousands to hundreds of thousands, depending on the model resolution}. The TVOPs can thus be defined as the matrices among the T matrices obtained via the pre-computations that provide a lower-bound and upper-bound of the maximum temperature rise obtained over the whole model.
[0053] It is noted that the T matrices, and thus the TVOPs, can be computed in a number of different ways. In some embodiments, the TVOPs can be computed based on measurements of and B0 and using model-based E -field values to calculate E- fields. This approach is non-specific to the subject being imaged. In some other embodiments, the TVOPs can be computed based on measurements of and B0 and using a E -field algorithm to generate E -field measurements specific to the subject being imaged. As one example, the E -field algorithm can be the one described in copending US Patent Application Serial No. 14/245,145, which is herein incorporated by reference in its entirety.
[0054] By defining ATG as the average temperature rise (with a corresponding temperature matrix] over the whole volume, the TVOPs can be defined as the set of T matrices, which guarantees that,
Ξ k& TVOPs such that ΔΤ^ < max(AT(r)) < ΔΤ^ + λ ΑΎ0 (14); [0055] where λ is a tolerance factor that controls the level of compression and
ATk denotes the temperature rise calculated for the kth TVOP. In general, larger values of λ will result in more compression (i.e., fewer TVOPs}. Instead of the matrix ATG , the identity matrix can also be used. In this case, the tolerance is expressed in terms of the average RF power.
[0056] Because the T matrices depend on the duration of the RF exposure, the
TVOPs will also depend on the duration of the RF exposure. By construction of the TVOPs, it is thereby ensured that the maximum temperature rise over the full model is always upper-bounded by one of the TVOPs plus the defined tolerance, and lower bounded by the maximum over the TVOPs.
[0057] To select the initial set of TVOPs, the T matrices are first sorted according to their maximum eigenvalues in decreasing order. The set of TVOPs is then initialized as the first T matrix in the sorted list.
[0058] Referring again to FIG. 2, the initial set of TVOPs can be further refined, as indicated at step 206. As one example, the set of TVOPs can be reduced based on a search for a set of N positive coefficients, cn , where N is the current number of TVOPs. In particular, the set of coefficients that is searched for has the sum,
N
∑c„= l (15). n
[0059] A candidate TVOP, T , is added to the initial set of TVOPs unless a set of coefficients that renders the following metric to be nonnegative can be found for that candidate TVOP:
N
∑c TVOP + λΎβ - Ύ (16); n
[0060] where T is the candidate TVOP and TG is the global temperature matrix
(returning the average temperature rise scalar quantity, ATG). If a set of coefficients, cn , is found such that the metric in Eqn. (16) is nonnegative, then the candidate T matrix is not retained as a TVOP; otherwise, the candidate T is selected as an additional TVOP. [0061] After the set of TVOPs has been refined, the temperature constraint is defined, as indicated at step 208. For example, the temperature constraint can be defined based on a selected maximum temperature value, Tmax, such as, ci(x)≤rM for i = l,. . . , iVrraRf (17).
[0062] The maximum temperature, Tmax, corresponds to the maximum allowable rise in temperature, and can be selected as any suitable value that constrains the temperature rise to safe values. For instance, Tmax can be in the range of about 1.0-2.0 degrees Celsius.
Examples
[0063] An example set of SAR virtual observation points ("VOPs"} computed for a subject is illustrated in FIG. 3A, and an example set of TVOPs determined according to the methods described here for the same subject is illustrated in FIG. 3B. In this example, 389 different SAR VOPs and TVOPs were determined over a whole model, and were obtained with a compression factor of = 10 and were determined for a thirty minute acquisition time. The dependence of temperature is multifactorial and does not simply depend on SAR hotspots. Thus, the locations of the TVOPs, shown as white squares in FIG. 3B, are not correlated with the locations of the SAR VOPs, shown as white squares in FIG. 3A. While the SAR VOPs can be located deep inside the brain, the TVOPs are mostly located in low-perfused regions (e.g., skin, muscle, cerebrospinal fluid}. The compressibility of the model improves with the duration of the sequence. For example, keeping λ = 10, the number of TVOPs in the example illustrated in FIG 3B (i.e., 389 TVOPs for a thirty minute sequence] increases to 483 for a ten minute sequence, and to 551 for a six minute sequence.
[0064] FIG. 4 illustrates a validation of an example set of TVOPs determined for a six minute sequence, with the temperature rise calculated over the TVOPs, the same temperature plus the defined tolerance, and the temperature calculated over the whole model. The temperature calculated over the whole model (shown as a dashed line] is upper-bounded by the TVOPs plus the selected tolerance, and is lower-bounded by the TVOPs. The calculation was performed for 20,000 random RF waveforms having a uniform distribution for amplitude and phase. The total power for each waveform was scaled to 10 W. In addition, FIG. 5 depicts histograms of the maximum temperature rises for six, ten, and thirty minute RF exposures, computed with the same 20,000 random RF scenarios. In FIG. 5, however, the results were scaled so that the peak 10-g SAR was always equal to 10 W/kg. The global SAR threshold of 3.2 W/kg was never reached in the results shown in FIG. 5. Despite the peak 10-g SAR being always the same, these examples illustrate that the maximum temperature rise can reach a range of different values, and show that there is no one-to-one mapping between the peak 10-g SAR and the maximum temperature rise. Also, the longer the sequence, the wider the temperature distribution becomes, with a double hump furthermore appearing for the thirty min sequence.
[0065] Table 1 summarizes the results of an example pulse design using the A-S algorithm under various constraints, for six, ten, and thirty minutes time-of-flight sequences.
Table 1
Constraint NRMSE (%) SAR lOg (W/kg) SAR Global (W/kg) Temperature (°C)
6 min
Temperature 10.8 13.4 5.3 1
Temperature 8.5 21 6.9 1.5
Temperature 8.1 27.3 7.6 2
Temperature 11.4 12.6 4.9 0.91
SAR 9.7 13.4 5.3 1.28
SAR 8.3 21 6.9 1.93
SAR 8.1 27.3 7.6 2.38
SAR (IEC) 16.2 10 3.2 0.91
10 min
Temperature 12.6 11.5 4.5 1
Temperature 9.0 17.8 6.3 1.5
Temperature 8.3 24.7 7.3 2
Temperature 11.4 12.4 4.8 1.1
SAR 11.2 11.5 4.5 1.38
SAR 8.5 17.8 6.3 2.07
SAR 8.2 24.7 7.3 2.80
SAR (IEC) 16.2 10 3.2 1.10
30 min
Temperature 15.5 10.7 3.9 1
Temperature 10.6 16.0 5.5 1.5
Temperature 8.7 21.8 6.7 2 Temperature 11.3 14.4 5.1 1.37
SAR 12.9 10.7 3.9 1.60
SAR 9.2 16.0 5.5 2.34
SAR 8.2 21.8 6.7 2.79
SAR (IEC) 16.2 10 3.2 1.37
[0066] The A-S algorithm converged to a local minimum after 300 iterations in about 30 seconds for the SAR and temperature-constrained pulse designs. The peak power constraint (2 kW in this example] was always saturated. For sequence durations up to 30 minutes, the pulse design under explicit temperature rise constraints noticeably outperformed the pulse design constrained under SAR, with a lower Normalized Root Mean Square Error ("NRMSE"}.
[0067] Examples of flip angle maps and temperature rise maps determined for the 6 minute sequence are illustrated in FIGS. 6 and 7, respectively. While the target flip angles in FIG. 6 cannot be achieved in the SAR-constraint pulse design case, the correct flip angles can be obtained in the temperature rise constrained case as the temperature rise can be progressively increased in a controlled manner. The corresponding temperature rises are shown in FIG. 7, where the hot spots mainly occur in the skin and CSF regions.
[0068] An efficient algorithm for directly constraining the temperature rise in RF pulse design using parallel transmission ("pTx"}, which is applicable for high and ultrahigh magnetic field strengths, has been provided. Although the international guidelines generally involve absolute temperature constraints, for which some approaches have been proposed, constraining the temperature rise provides a more conservative and practical approach for designing safe and efficient RF pulses. The systems and methods provided here can allow for a direct control of a more relevant safety metric in RF pulse design than the ubiquitously used SAR.
[0069] While all the considered sequence durations with peak 10-g SAR equal to
10 W/kg could lead to temperature rises higher than 1 °C (see FIG. 5], the results provided in Table 1 show that the pulse design techniques described here can be more powerful and safer by achieving lower NRMSEs with more controlled temperature rises. The same table also shows that SAR constraints can be relaxed to achieve better RF pulse performance, but then with a loss of direct control over the temperature, which increases with sequence duration.
[0070] Compared to SAR-constrained techniques, the TVOP approach described here may require more computing efforts because a different set of TVOPs exists for each pulse sequence duration. However, when generic models are used, the TVOPs can be pre-computed only once.
[0071] The locations of the TVOPs (see FIG. 3B] also indicate that the maximum temperature rises occur in less sensitive regions (e.g., in the skin] than the brain or other internal structures, where the equilibrium temperature (i.e. the steady state obtained without RF source] is less than 37°C. As such, it is possible that a 1.0-2.0°C temperature rise constraint can be implemented while remaining conservative with respect to the absolute temperature guidelines.
[0072] Because the MLS optimization problem is non-convex, it is uncertain whether a given temperature-constrained or SAR-constrained design corresponds to the global minimum (i.e., the most flip-angle homogenizing solution]. The advantages of temperature-regulated RF pulse design numerically demonstrated above appears as a consequence of the multi-factorial dependence of temperature.
[0073] The model described above for the temperature calculations was Pennes' bio-heat model, but the same method can be applied for any linear model. Taking into account a possible dependence of the thermal parameters on temperature would break the linearity; however, disregarding this dependence puts the thermal calculations on the safe side.
[0074] The proposed pulse design method also presents a natural increase in complexity compared to SAR-constrained optimizations because the sliding-integration window used for SAR calculations can no longer be applied, which is especially problematic when different sequences using different RF pulses follow each other. To address this, for temperature-constraint pulse designs, a pause can be interleaved between different pulse sequences to reset the temperature, or the whole imaging exam can be taken into account.
[0075] On a more practical note, because SAR hotspots in pTx spatially depend on the designed pulses, and hence the sequences, the resulting temperature rise is at its most additive in the worst-case scenario, and less so in all other cases. As a result, pulses for two sequences can be designed independently, each with a different temperature rise constraint and, thus, a resulting overall temperature rise being upper- bounded by the sum of the two constraints. This design approach offers more latitude in RF pulse design than if SAR constraints were used.
[0076] Referring particularly now to FIG. 8, an example of a magnetic resonance imaging ("MRI"} system 800 is illustrated. The MRI system 800 includes an operator workstation 802, which will typically include a display 804; one or more input devices 806, such as a keyboard and mouse; and a processor 808. The processor 808 may include a commercially available programmable machine running a commercially available operating system. The operator workstation 802 provides the operator interface that enables scan prescriptions to be entered into the MRI system 800. In general, the operator workstation 802 may be coupled to four servers: a pulse sequence server 810; a data acquisition server 812; a data processing server 814; and a data store server 816. The operator workstation 802 and each server 810, 812, 814, and 816 are connected to communicate with each other. For example, the servers 810, 812, 814, and 816 may be connected via a communication system 840, which may include any suitable network connection, whether wired, wireless, or a combination of both. As an example, the communication system 840 may include both proprietary or dedicated networks, as well as open networks, such as the internet.
[0077] The pulse sequence server 810 functions in response to instructions downloaded from the operator workstation 802 to operate a gradient system 818 and a radiofrequency ("RF"} system 820. Gradient waveforms necessary to perform the prescribed scan are produced and applied to the gradient system 818, which excites
G G
gradient coils in an assembly 822 to produce the magnetic field gradients x , y , and z used for position encoding magnetic resonance signals. The gradient coil assembly 822 forms part of a magnet assembly 824 that includes a polarizing magnet 826 and a whole-body RF coil 828.
[0078] RF waveforms are applied by the RF system 820 to the RF coil 828, or a separate local coil (not shown in FIG. 8], in order to perform the prescribed magnetic resonance pulse sequence. Responsive magnetic resonance signals detected by the RF coil 828, or a separate local coil (not shown in FIG. 8], are received by the RF system 820, where they are amplified, demodulated, filtered, and digitized under direction of commands produced by the pulse sequence server 810. The RF system 820 includes an RF transmitter for producing a wide variety of RF pulses used in MRI pulse sequences. The RF transmitter is responsive to the scan prescription and direction from the pulse sequence server 810 to produce RF pulses of the desired frequency, phase, and pulse amplitude waveform. The generated RF pulses may be applied to the whole-body RF coil 828 or to one or more local coils or coil arrays (not shown in FIG. 8}.
[0079] The RF system 820 also includes one or more RF receiver channels. Each
RF receiver channel includes an RF preamplifier that amplifies the magnetic resonance signal received by the coil 828 to which it is connected, and a detector that detects and digitizes the I and ^ quadrature components of the received magnetic resonance signal. The magnitude of the received magnetic resonance signal may, therefore, be determined at any sampled point by the square root of the sum of the squares of the I and ^ components:
Figure imgf000018_0001
(18);
[0080] and the phase of the received magnetic resonance signal may also be determined according to the following relationship:
Figure imgf000018_0002
[0081] The pulse sequence server 810 also optionally receives patient data from a physiological acquisition controller 830. By way of example, the physiological acquisition controller 830 may receive signals from a number of different sensors connected to the patient, such as electrocardiograph ("ECG") signals from electrodes, or respiratory signals from a respiratory bellows or other respiratory monitoring device. Such signals are typically used by the pulse sequence server 810 to synchronize, or "gate," the performance of the scan with the subject's heart beat or respiration.
[0082] The pulse sequence server 810 also connects to a scan room interface circuit 832 that receives signals from various sensors associated with the condition of the patient and the magnet system. It is also through the scan room interface circuit 832 that a patient positioning system 834 receives commands to move the patient to desired positions during the scan.
[0083] The digitized magnetic resonance signal samples produced by the RF system 820 are received by the data acquisition server 812. The data acquisition server 812 operates in response to instructions downloaded from the operator workstation 802 to receive the real-time magnetic resonance data and provide buffer storage, such that no data is lost by data overrun. In some scans, the data acquisition server 812 does little more than pass the acquired magnetic resonance data to the data processor server 814. However, in scans that require information derived from acquired magnetic resonance data to control the further performance of the scan, the data acquisition server 812 is programmed to produce such information and convey it to the pulse sequence server 810. For example, during prescans, magnetic resonance data is acquired and used to calibrate the pulse sequence performed by the pulse sequence server 810. As another example, navigator signals may be acquired and used to adjust the operating parameters of the RF system 820 or the gradient system 818, or to control the view order in which k-space is sampled. In still another example, the data acquisition server 812 may also be employed to process magnetic resonance signals used to detect the arrival of a contrast agent in a magnetic resonance angiography ("MRA"} scan. By way of example, the data acquisition server 812 acquires magnetic resonance data and processes it in real-time to produce information that is used to control the scan.
[0084] The data processing server 814 receives magnetic resonance data from the data acquisition server 812 and processes it in accordance with instructions downloaded from the operator workstation 802. Such processing may, for example, include one or more of the following: reconstructing two-dimensional or three- dimensional images by performing a Fourier transformation of raw k-space data; performing other image reconstruction algorithms, such as iterative or backprojection reconstruction algorithms; applying filters to raw k-space data or to reconstructed images; generating functional magnetic resonance images; calculating motion or flow images; and so on.
[0085] Images reconstructed by the data processing server 814 are conveyed back to the operator workstation 802 where they are stored. Real-time images are stored in a data base memory cache (not shown in FIG. 8], from which they may be output to operator display 802 or a display 836 that is located near the magnet assembly 824 for use by attending physicians. Batch mode images or selected real time images are stored in a host database on disc storage 838. When such images have been reconstructed and transferred to storage, the data processing server 814 notifies the data store server 816 on the operator workstation 802. The operator workstation 802 may be used by an operator to archive the images, produce films, or send the images via a network to other facilities.
[0086] The MRI system 800 may also include one or more networked workstations 842. By way of example, a networked workstation 842 may include a display 844; one or more input devices 846, such as a keyboard and mouse; and a processor 848. The networked workstation 842 may be located within the same facility as the operator workstation 802, or in a different facility, such as a different healthcare institution or clinic.
[0087] The networked workstation 842, whether within the same facility or in a different facility as the operator workstation 802, may gain remote access to the data processing server 814 or data store server 816 via the communication system 840. Accordingly, multiple networked workstations 842 may have access to the data processing server 814 and the data store server 816. In this manner, magnetic resonance data, reconstructed images, or other data may be exchanged between the data processing server 814 or the data store server 816 and the networked workstations 842, such that the data or images may be remotely processed by a networked workstation 842. This data may be exchanged in any suitable format, such as in accordance with the transmission control protocol ("TCP"}, the internet protocol ("IP"}, or other known or suitable protocols.
[0088] The present invention has been described in terms of one or more preferred embodiments, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention.

Claims

We claim:
1. A method for designing parallel transmission radio frequency (RF] pulses for use in magnetic resonance imaging (MRI], the steps of the method comprising:
(a] providing to a computer system, a set of pre-computed temperature matrices that describe temperature rises in response to RF pulses based on a linear thermal model;
(b] selecting with the computer system, a set of temperature virtual observation points (TVOPs] based on the pre-computed temperature matrices;
(c] computing a plurality of RF pulse waveforms for use in MRI by solving, with the computer system, an optimization problem that is constrained by a maximum temperature rise evaluated at the set of TVOPs.
2. The method as recited in claim 1, wherein the pre-computed temperature matrices are computed based in part on a specific absorption rate (SAR] values computed from simulated electric field distributions.
3. The method as recited in claim 2, wherein the simulated electric field distributions are filtered to remove electric field spikes before computing the SAR values.
4. The method as recited in claim 2, wherein the simulated electric field distributions are computed for a particular subject to be imaged.
5. The method as recited in claim 1, wherein the pre-computed temperature matrices are computed for a particular pulse sequence duration.
6. The method as recited in claim 1, wherein step (b] includes selecting as the TVOPs the pre-computed temperature matrices that satisfy,
AT,≤max(AT(r))≤AT, + - ATG wherein ATk is a temperature rise calculated for a kth TVOP, ΔΤ(Γ) is a temperature rise indicated in a pre-computed temperature matrix, ATG is an average temperature rise over a volume-of-interest, and λ is a tolerance factor that controls a level of compression of the T VOPs.
7. The method as recited in claim 6, wherein ATG is selected to be an identify matrix.
8. The method as recited in claim 1, wherein selecting the set of T VOPs includes determining eigenvalues for the pre-computed temperature matrices and sorting the pre-computed temperature matrices according the maximum eigenvalues in decreasing order.
9. The method as recited in claim 1, wherein the optimization problem is a magnitude least squares problem.
10. The method as recited in claim 1, wherein solving the optimization problem includes computing the plurality of RF pulse waveforms that minimizes an objective function that is based in part on a selected encoding matrix and a selected flip angle.
11. The method as recited in claim 1, wherein step (c] includes selecting at least one imaging region to which RF pulses should be applied and the objective function includes an encoding matrix and a flip angle for each imaging region.
12. The method as recited in claim 1, wherein the optimization problem is further constrained by at least one of a maximum peak power constraint or an average power constraint.
13. The method as recited in claim 1, further comprising providing the plurality of RF pulse waveforms to an MRI system and controlling the MRI system to generate RF pulses using the plurality of RF pulse waveforms, wherein the RF pulses generate a rise in temperature that is below the maximum temperature rise.
PCT/US2015/051425 2014-09-22 2015-09-22 Local temperature rise constrained radio frequency pulse design in parallel transmission mri WO2016049013A1 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US15/513,042 US10545210B2 (en) 2014-09-22 2015-09-22 Local temperature rise constrained radio frequency pulse design in parallel transmission

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US201462053603P 2014-09-22 2014-09-22
US62/053,603 2014-09-22

Publications (1)

Publication Number Publication Date
WO2016049013A1 true WO2016049013A1 (en) 2016-03-31

Family

ID=54251745

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/US2015/051425 WO2016049013A1 (en) 2014-09-22 2015-09-22 Local temperature rise constrained radio frequency pulse design in parallel transmission mri

Country Status (2)

Country Link
US (1) US10545210B2 (en)
WO (1) WO2016049013A1 (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10878639B2 (en) * 2018-02-09 2020-12-29 David Byron Douglas Interactive voxel manipulation in volumetric medical imaging for virtual motion, deformable tissue, and virtual radiological dissection
US10247803B2 (en) * 2014-04-25 2019-04-02 Regents Of The University Of Minnesota Systems and methods for designing magnetic resonance imaging radio frequency pulses that are robust against physiological motion errors
EP3637124A1 (en) 2018-10-11 2020-04-15 Siemens Healthcare GmbH Method for controlling a high-frequency amplifier for a magnetic resonance system

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100134105A1 (en) * 2008-10-15 2010-06-03 Zelinski Adam C Method For Reducing Maximum Local Specific Absorption Rate In Magnetic Resonance Imaging
US20120256626A1 (en) * 2011-04-08 2012-10-11 Siemens Aktiengesellschaft Parallel transmission rf pulse design with local sar constraints
US20140300354A1 (en) * 2013-04-04 2014-10-09 Bin He Systems and methods for spatial gradient-based electrical property properties tomography using magnetic resonance imaging

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10661091B2 (en) * 2014-03-26 2020-05-26 The General Hospital Corporation System and method for hyperthermia treatment using radiofrequency phased arrays

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100134105A1 (en) * 2008-10-15 2010-06-03 Zelinski Adam C Method For Reducing Maximum Local Specific Absorption Rate In Magnetic Resonance Imaging
US20120256626A1 (en) * 2011-04-08 2012-10-11 Siemens Aktiengesellschaft Parallel transmission rf pulse design with local sar constraints
US20140300354A1 (en) * 2013-04-04 2014-10-09 Bin He Systems and methods for spatial gradient-based electrical property properties tomography using magnetic resonance imaging

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
A. MASSIRE ET AL.: "Thermal simulations in the human head for high field MRI using parallel transmission", J MAGN RESON IMAGING, vol. 35, 2012, pages 1312 - 1321
AURÉLIEN MASSIRE ET AL: "Thermal simulations in the human head for high field MRI using parallel transmission", JOURNAL OF MAGNETIC RESONANCE IMAGING, vol. 35, no. 6, 12 January 2012 (2012-01-12), US, pages 1312 - 1321, XP055228455, ISSN: 1053-1807, DOI: 10.1002/jmri.23542 *
CEM MURAT DENIZ ET AL: "Specific absorption rate benefits of including measured electric field interactions in parallel excitation pulse design", MAGNETIC RESONANCE IN MEDICINE, vol. 67, no. 1, 29 August 2011 (2011-08-29), pages 164 - 174, XP055046641, ISSN: 0740-3194, DOI: 10.1002/mrm.23004 *
CLOOS M A ET AL: "Local SAR reduction in parallel excitation based on channel-dependent Tikhonov parameters", JOURNAL OF MAGNETIC RESONANCE IMAGING, SOCIETY FOR MAGNETIC RESONANCE IMAGING, OAK BROOK, IL, US, vol. 32, 28 October 2010 (2010-10-28), pages 1209 - 1216, XP002610147, ISSN: 1053-1807, DOI: 10.1002/JMRI.22346 *
ESRA NEUFELD ET AL: "Analysis of the local worst-case SAR exposure caused by an MRI multi-transmit body coil in anatomical models of the human body;Local worst-case SAR exposure caused by an MRI multi-transmit body coil", PHYSICS IN MEDICINE AND BIOLOGY, INSTITUTE OF PHYSICS PUBLISHING, BRISTOL GB, vol. 56, no. 15, 6 July 2011 (2011-07-06), pages 4649 - 4659, XP020208285, ISSN: 0031-9155, DOI: 10.1088/0031-9155/56/15/002 *
GRAESSLIN I ET AL: "SAR hotspot reduction by temporal averaging in parallel transmission", INTERNATIONAL SOCIETY FOR MAGNETIC RESONANCE IN MEDICINE. SCIENTIFIC MEETING AND EXHIBITION. PROCEEDINGS, INTERNATIONAL SOCIETY FOR MAGNETIC RESONANCE IN MEDICINE, US, vol. 17, 18 April 2009 (2009-04-18), pages 176, XP002564610, ISSN: 1524-6965 *
YUDONG ZHU ET AL: "System and SAR characterization in parallel RF transmission", MAGNETIC RESONANCE IN MEDICINE, vol. 67, no. 5, 2 December 2011 (2011-12-02), pages 1367 - 1378, XP055027275, ISSN: 0740-3194, DOI: 10.1002/mrm.23126 *

Also Published As

Publication number Publication date
US20170307710A1 (en) 2017-10-26
US10545210B2 (en) 2020-01-28

Similar Documents

Publication Publication Date Title
KR102354723B1 (en) tensor field mapping
JP6142073B2 (en) Radiotherapy system with real-time magnetic resonance monitoring
WO2017106469A1 (en) Systems and methods for analyzing perfusion-weighted medical imaging using deep neural networks
US9928596B2 (en) Motion corrected imaging system
US20180238983A1 (en) Systems and methods for multislice magentic resonance fingerprinting
US20170146623A1 (en) Systems and Methods For Segmented Magnetic Resonance Fingerprinting Dictionary Matching
WO2013169368A1 (en) System and method for local sar reduction in multislice parallel transmission magnetic resonance imaging using sar hopping between excitations
US11051711B2 (en) Noninvasive determination of electrical properties of tissues and materials using magnetic resonance measurements
US11875509B2 (en) System and method for determining undersampling errors for a magnetic resonance fingerprinting pulse sequence
EP3151027B1 (en) Mri method for jointly reconstructing an image and estimating the actual k-space trajectory
US10545210B2 (en) Local temperature rise constrained radio frequency pulse design in parallel transmission
CN107427257A (en) The magnetic resonance imaging temperature measuring measured using Proton Resonance Frequency and T1
US9041395B2 (en) MRI method of calculating and generating spatially-tailored parallel radio frequency saturation fields
US20150196223A1 (en) System and method for neutral contrast magnetic resonance imaging of calcifications
US10330752B2 (en) Magnetic resonance imaging apparatus and method for setting RF shimming parameters
US20180292491A1 (en) Systems and methods for calibrated multi-spectral magnetic resonance imaging
CN110709721B (en) System and method for simultaneous multi-slice magnetic resonance fingerprint imaging using parallel transmit radio frequency coil arrays
US10884091B2 (en) Voxelwise spectral profile modeling for use in multispectral magnetic resonance imaging
JP7420834B2 (en) Method of motion artifact detection
US10545209B2 (en) System and method for acquiring both T2*-weighted and T1-weighted data in a single acquisition using a single dose of contrast agent
US20220248973A1 (en) Deep learning of eletrical properties tomography
US20230280428A1 (en) Method for operating a magnetic resonance imaging scanner, magnetic resonance imaging scanner, computer program and storage medium with the computer program
WO2023049524A1 (en) Parallel transmit radio frequency pulse design with deep learning
WO2023205497A1 (en) Multi-spectral susceptibility-weighted magnetic resonance imaging

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 15775050

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 17.07.2017)

122 Ep: pct application non-entry in european phase

Ref document number: 15775050

Country of ref document: EP

Kind code of ref document: A1