Efficient prediction of MRI gradient-induced heating for guiding safety testing of conductive implants

Purpose: To propose an efficient numerical method to predict the temperature increase of an implantable medical device induced by any linearly polarized homogeneous magnetic field, according to the ISO 10974 methodology for testing of gradient-induced device heating. Theory and Methods: The concepts of device-specific power and temperature tensors are introduced to mathematically describe the electromagnetic and thermal anisotropic behavior of the device, from which the device heating for an arbitrary exposure direction can be predicted. The proposed method is compared to a brute-force approach based on simulations, and validated by applying it to four reference orthopedic implants with a commercial simulation software. Results: The proposed method requires about 5 % of the time required by the brute-force approach, and 30 % of the memory occupancy. The temperature increase predicted by the proposed method over a range of incident magnetic field exposures deviated from brute-force direct simulations by less than ± 0.3 % . Conclusion: The proposed method allows efficient prediction of the heating of an implantable medical device induced by any linearly polarized homogeneous magnetic field using a small fraction of the simulations required by the brute-force approach. The results can be


INTRODUCTION
MRI is a largely used imaging tool for noninvasive clinical diagnosis. Despite significant efforts in making MRI an increasingly safe technique, many safety issues are still being addressed, especially when it comes to patients carrying implantable medical devices. [1][2][3][4][5] This is of particular concern due to the growing number of implanted orthopedic devices [6][7][8][9] and the related need of MRI examinations aimed at assessing possible perioperative and postoperative complications. 10 One main safety issue in this regard is MRI-induced heating of the tissues surrounding the implant. Radiofrequency (RF) tissue heating has gained a wide attention both from the literature and relevant standards due to its potential to cause significant temperature increases when the incident RF field, implant characteristics, and exposure scenario meet specific criteria. [11][12][13][14][15][16][17] While RF remains the most important source of heating in MRI, recent studies showed that heating induced by gradient fields could also represent a risk in the presence of bulky conductive objects. 5,18,19 The ISO 10974 standard 17 recognizes that the switching gradient fields can induce significant electric currents on an active implantable medical device enclosure. Such currents lead to heating of the device, with a consequent risk for the patient. For this reason, the standard proposes a specific test to experimentally assess the gradient heating risk. During the test, the active implantable medical device is exposed to a linearly polarized harmonic homogeneous magnetic field oriented along the direction which maximizes the heating. The maximum temperature increase on the active implantable medical device enclosure is then measured after a suggested time of 30 min of ongoing exposure. While the worst exposure direction in terms of heating can be identified rather easily for a disk-like active implantable medical device enclosure, it may be difficult to be identified when the same test is performed on objects with a more complex shape, such as many orthopedic implants. 20 In this case, the worst exposure direction has to be investigated through several different experiments where, in each experiment, the orientation of the device under test (DUT) is slightly changed and the temperature increase is recorded and compared with those obtained with the other orientations. While recording the temperature increase for all possible device orientations and all over its surface is not feasible, reducing the number of tested orientations and recording points does not assure that the worst device orientation and maximum temperature increase are properly identified. In this context, demonstrating compliance with MRI safety regulations and obtaining "MR conditional" labeling becomes even more challenging and expensive for implant manufacturers and testing companies. 21 Numerical simulations provide a valuable tool to address this issue whenever the virtual model of the DUT is available. However, approaching the problem through a brute-force strategy is still far from an efficient solution. Indeed, a large number of simulations is required to reliably identify the worst-case direction, costing significant computational time as well as storage of result files. 20,22 In this paper, we propose a strategy to predict the gradient-induced heating of any passive implant under any linearly-polarized homogeneous magnetic field by performing a minimum amount of numerical simulations to generate a number of 3 × 3 tensors which characterize the implant response. We validated the method through a comparison with a set of explicit gradient heating simulations representing a brute-force approach. We performed all simulations in Sim4Life, 23 and developed a Python script that can be used in Sim4Life to automatically compute the worst exposure direction after the virtual model of an implant is imported. The script is available in the authors' GitHub repository 24 along with basic documentation and a relevant example.

Power tensor
The current density j(x) generated in a spatial point x of the DUT by a linearly polarized harmonic homogeneous magnetic flux density B can be defined as where j u i (x) is a three element complex vector representing the (peak) harmonic current density generated in x by a unitary magnetic flux density oriented along the ith Cartesian direction and with the same frequency of B , B i is the ith component of the magnetic flux density B and J (x) is a 3 × 3 matrix with the j u i (x) vectors as columns. Since B is linearly polarized, we can assume that B ∈ R 3 .
The local power density p(x) dissipated in x by the Joule losses can be computed as where is the electrical conductivity in x , ℜ denotes the real part, and T and H superscripts denote the transpose and Hermitian transpose operators, respectively. The total power P deposited inside the DUT is obtained by integrating (2) over the DUT volume () as Here, Q is a device-specific 3 × 3 positive semidefinite symmetric tensor which relates the power deposited inside the DUT to the linearly polarized harmonic homogeneous magnetic flux density which is responsible for the power deposition. The tensor Q can be computed by means of three electromagnetic simulations aiming at obtaining the current densities j u i (x). The tensor Q describes the anisotropy of the DUT with respect to the electromagnetic problem.

Temperature tensor field
Since the power density distribution inside the DUT and the resulting temperature increase are proportional in a phantom experiment, it is possible to express the latter with a quadratic form similar to that obtained for the power deposition. From (2), it is possible to express the temperature rise in each voxel of the implant as where  T is the linear operator solving the heat equation, Therefore,  T depends on the thermal conductivity , the mass density and the specific heat capacity c p of the DUT and the surrounding environment. The linearity of  T allows introducing M (x), device-specific 3 × 3 positive semidefinite symmetric tensors which relate the temperature increase in any spatial point x of the DUT with the linearly polarized harmonic homogeneous magnetic field responsible for it. While the tensor Q alone represents the whole DUT with regard to power deposition, the temperature is not a global quantity, and the implant is characterized by a tensor field, where each M (x) tensor is relevant to the temperature increase at a specific position in the DUT. Considering that the tensors M (x) are symmetric, the whole tensor field can be computed with six thermal simulations. Indeed, the m, n component of M (x) is equal to the temperature rise distribution following the application of a fictitious power density The tensor M (x) describes the anisotropy of the DUT in x with respect to the complete (electromagnetic and thermal) problem.

Worst exposure orientation
With reference to the quadratic form in (3), the maximum power deposition inside the DUT following the application of a linearly polarized unitary harmonic homogeneous magnetic flux density is given by the largest eigenvalue of Q . The maximum power deposition is obtained with a magnetic flux density oriented as an eigenvector associated with the largest eigenvalue (see Appendix S1). Similarly, the maximum temperature increase in a DUT point x is given by the largest eigenvalue of M(x) and it is induced by a magnetic flux density oriented as one of the corresponding eigenvectors. The overall maximum temperature increase in the DUT and the relevant magnetic flux density direction can be retrieved by identifying the maximum among the maximum eigenvalues of the M(x) tensors.

METHODS
The above procedure was implemented in Sim4Life 23 and applied on four orthopedic implants, described in Table 1, in order to identify the exposure conditions that maximize the power deposition and the temperature increase. To validate the method, the deposited power and the resulting maximum temperature increase of the shoulder implant were computed in Sim4Life, by exposing the implant to a homogeneous magnetic flux density with amplitude equal to 35 mT and frequency equal to 270 Hz (i.e., the values proposed in ISO 10974 17 ) while varying the magnetic flux density direction to cover the surface of a hemisphere with angular steps of 22.5 • . After that, the Q and M (x) tensors were computed and, for each magnetic flux density direction, the power deposited into the shoulder implant and the maximum temperature increase previously obtained were compared to the same quantities computed according to (3) and (4). Eventually, the worst exposure orientation deduced by the proposed method, and the corresponding result, were compared to those obtained by the brute-force approach.
To comply with the requirements of ISO 10974, 17 the thermal simulations, required to compute the M (x) tensor field, were performed with the implants centrally positioned inside a phantom. The dimensions of the phantom (130 mm × 200 mm × 400 mm) guaranteed that the phantom-air boundary conditions did not play a significant role in the thermal results for the 30min simulated thermal evolution (see Figure S1). The thermal properties of the phantom (specific heat capacity: 4150 J kg −1 K −1 , thermal conductivity: 0.54 W m −1 K −1 ) complied with those reported in the ASTM standard 16 for the polyacrylic acid-based solution. The phantom was not included in the electromagnetic simulations performed to compute the Q tensor, since its low electrical conductivity (with respect to that of the conductive implants) would not significantly affect the deposited power. 25 All simulations were performed using an isotropic 1 mm 3 mesh on an Intel(R) Xeon(R) E5-2650 CPU with 64 GB of RAM.

RESULTS
The validation of the proposed method is illustrated in Figure 1,  The same difference is obtained looking at the maximum temperature increases after 30-min exposure. Figure 2 represents, for each implant, the worst exposure orientations, estimated with the proposed approach, which maximize the power deposition inside the implant (blue arrow) and the peak temperature increase after 30-min of exposure (red arrow).
The first row of Figure 3 shows the maximum temperature increase that each voxel of the implants can reach when exposed to a homogeneous magnetic flux density of 35 mT at 270 Hz. The second row of Figure 3 shows the distributions of the temperature increase induced in each implant exposed to the homogeneous magnetic flux density of 35 mT at 270 Hz oriented along the direction which maximizes the peak temperature increase within the implant.
The orientation of the magnetic flux density which results in maximum temperature increase at each point is illustrated in Figure 4. The arrow colors correspond to those in the second row of Figure 3.

DISCUSSION
The negligible difference in power and peak temperature estimates between the brute-force approach and the proposed method demonstrates the reliability of the proposed procedure. In order to compute the M (x) tensor field, the computations required about 4 h to perform the three electromagnetic and six thermal simulations, yielding 4.5 GB of data. The brute-force approach required almost 80 h of computations in order to perform the 60 electromagnetic and thermal simulations, producing 13.4 GB of data (noting that both times could have been further shortened by enabling GPU acceleration). Finally, the proposed method determines the actual maximum power deposition and peak temperature increase. With a 22.5 • angular step, the brute-force approach, used for validating the proposed method, estimates a maximum power deposition of 0.577 W and a maximum temperature increase of 2.53 K. In both cases, these values are obtained for magnetic flux density oriented according to = 45 • and = −22.5 • ). The proposed procedure computes a maximum power deposition of 0.582W (for = 47.41 • and = −27.46 • ) and a maximum temperature increase of 2.54 K (for = 47.28 • and = −30.46 • ). In this case, the brute-force approach revealed to be reliable but this may not be the case in general. In fact, the discretization of the magnetic flux density direction sweep could lead to an underestimation of the above quantities if an inadequate discretization step is selected.
The magnetic flux density directions which maximize the power deposition and the peak temperature increase are slightly different (Figure 2), especially for the knee implant ( Figure 2C) and the ankle plate ( Figure 2D). Looking at the angle formed by the two magnetic flux density directions, we obtain a very good direction coherence for the shoulder and hip implants (3.89 • and 8.39 • , respectively) and a slightly worse coherence for the knee implant and the ankle plate (15.55 • and 15.5 • , respectively). The relative difference between the maximum peak temperature increases and those resulting from the exposition of the implant to a homogeneous magnetic flux density oriented according the direction which maximize the power deposition are equal to 0.8 % for the shoulder implant, 1.27 % for the hip implant, 3.7 % for the knee implant and 5.85 % for the ankle plate. Figure 3 shows that, while there is a direction of the magnetic flux density which maximizes the peak temperature increase within the implant, peak temperatures close

F I G U R E 3
First row: Chromatic maps representing, in each voxel, the temperature increase following the exposure to the homogeneous magnetic flux density whose direction maximizes the temperature increase in the same voxel. Second row: Temperature increase distributions resulting from the exposure to the homogeneous magnetic flux density whose direction maximizes the peak temperature increase in the implant. Results are shown for the shoulder implant (A), hip implant (B), knee implant (C), and ankle plate (D). All temperature increases are after 30 min of exposure to a 35 mT magnetic flux density with a 270 Hz frequency.

F I G U R E 4
Vector field plots representing the direction of the magnetic flux density which maximizes the temperature increase at the arrow position. The arrow colors correspond to those in the second row of to the maximum value can be attained by significantly different directions (see also Figure 4). This is particularly true for the knee implant ( Figure 3C) and the ankle plate ( Figure 3D). Dealing with the knee implant exposed to the magnetic flux density oriented along the worst direction (bottom row), the femoral component shows a temperature rise much lower than that encountered on the tibial component. However, the top row of Figure 3C shows that the femoral component can heat up quite as much as the tibial component for certain magnetic flux density directions. The ankle plate is seen to behave similarly when comparing the top and bottom rows of Figure 3D.
This information can help to define the position of the temperature probes in the ISO 10974 experiment and, combining the information contained in Figure 4, to identify further implant orientations to be tested for backup experiments. Indeed, if the digital model of the implant differs from the physical model, the point-wise worst possible temperature increase (second row of Figure 3) indicates locations which could be sensitive to variation.
Regarding the skin effect: the Sim4Life magneto quasi-static solver, used to compute the Q matrix, solves the Maxwell equations neglecting the reaction due to the eddy-currents induced into the conductive components of the implant. By comparison with a homemade code, 26,27 we found that the skin effect does not play a significant role for the considered implants at 270 Hz. However, this would have to be rechecked for higher test frequencies, such as the 1750 Hz proposed in the latest draft edition of ISO 10974. 20 Another direct consequence of the low frequency value is that the power deposited into the phantom is negligible with respect to that into the conductive implant. This behavior is typical for gradient coil-induced heating which sensibly differs from RF induced heating where the majority of the power is deposited into the tissues surrounding the implant. 5,25 This allowed us to save some disk space by removing the phantom in the electromagnetic simulations after verifying, for a selected case, that the results were not affected.
Finally, the proposed method allows to easily calculate the implant-related anisotropy coefficients proposed by Arduino et al., 22 used to avoid over-conservative safety limits when dealing with gradient-induced heating of implantable medical devices.

CONCLUSION
We have presented a numerical method to efficiently predict the power deposition and temperature increase inside a conductive medical implant experiencing a linearly polarized harmonic homogeneous magnetic flux density. The proposed strategy reduces the efforts and time required to perform gradient-heating measurements of passive medical implants following the ISO 10974 test methodology, compared to a brute-force approach. The proposed method has been implemented in Sim4Life, successfully validated, and applied to four orthopedic implants. A Sim4Life Python script that implements the method and automatically computes the orientations that maximize the deposited power and the induced temperature increase is provided in the authors' GitHub repository. 24 The script also provides functions for visualizing the results, and it is provided together with basic documentation and a usage example.

ACKNOWLEDGMENTS
The project (21NRM05 STASIS) has received funding from the European Partnership on Metrology, co-financed from the European Union's Horizon Europe Research and Innovation Programme and by the Participating States. Hip, knee and shoulder implants and related digital models have been kindly provided by Laboratorio di Tecnologia Medica, IRCCS Istituto Ortopedico Rizzoli, Bologna, Italy. The ankle plate and related digital model has been kindly provided by Medartis AG. The authors sincerely thank Dr. Tolga Goren and Dr. Lena Kranold for the valuable discussions and the contributions provided for the language editing of the paper.