Accelerating implant RF safety assessment using a low‐rank inverse update method

Purpose Patients who have medical metallic implants, e.g. orthopaedic implants and pacemakers, often cannot undergo an MRI exam. One of the largest risks is tissue heating due to the radio frequency (RF) fields. The RF safety assessment of implants is computationally demanding. This is due to the large dimensions of the transmit coil compared to the very detailed geometry of an implant. Methods In this work, we explore a faster computational method for the RF safety assessment of implants that exploits the small geometry. The method requires the RF field without an implant as a basis and calculates the perturbation that the implant induces. The inputs for this method are the incident fields and a library matrix that contains the RF field response of every edge an implant can occupy. Through a low‐rank inverse update, using the Sherman–Woodbury–Morrison matrix identity, the EM response of arbitrary implants can be computed within seconds. We compare the solution from full‐wave simulations with the results from the presented method, for two implant geometries. Results From the comparison, we found that the resulting electric and magnetic fields are numerically equivalent (maximum error of 1.35%). However, the computation was between 171 to 2478 times faster than the corresponding GPU accelerated full‐wave simulation. Conclusions The presented method enables for rapid and efficient evaluation of the RF fields near implants and might enable situation‐specific scanning conditions.


| INTRODUCTION
The group of patients with medical implants that require an MRI scan is constantly growing. However, MRI scanning of a patient with metallic implants bears a potentially severe safety risk. The electromagnetic (EM) fields produced by an MRI scanner can couple to the metallic implant resulting in image degradation and serious health hazards. The largest risk is tissue heating due to the radio frequency (RF) fields. The implant can locally enhance the RF fields causing temperature hotspots 1,2 with potentially severe consequences. 3,4 Therefore, people with an implant are either exempted from MRI scanning or scanned with very conservative RF power limitations degrading the achievable image quality severely.
In order to quantify the risks associated with a particular implant EM simulations are often performed. These EM simulations can compute the RF fields for a given transmit coil, patient model, and implant. The electrical properties, conductivity and permittivity, of both the patient 5 and implant 6,7 affect the resulting RF fields. The geometry and location with respect to each other of the transmit coil, patient model, and implant are relevant for assessing the RF fields. Finally, the drive settings for the transmit coil are required to correctly quantify the RF fields. [8][9][10][11][12][13] The simulated RF fields are often used in thermal simulations to quantify tissue heating.
A popular method for EM simulations is the finite-difference time-domain (FDTD) method. 14 With the FDTD method a single configuration of source, patient and implant can be computed at a time. These FDTD simulations are time-consuming due to the large domain (the whole MRI RF system) that needs to be simulated, even though the implant only affects a small domain within the patient. On top of this, hundredths of thousands of these FDTD simulations, for all the different possible configurations, are required to obtain a conditional label for an implant with the most lenient restrictions on scanning (i.e. the tier 4 approach as specified by ISO/ TS 10974:2018 15 ). For this reason, investigating RF safety for a particular implant in a patient model is a computationally demanding task. This has been demonstrated by B. Guerin et al 16 recently for different deep brain stimulation implants. The full-wave simulation, performed with the finite element method, took up to 6 hours with 13 processors and ∼300 GB of RAM for a single simulation. For FDTD simulations, E. Cabot et al 17 showed that similar types of simulations can take up to 43 hours for a single simulation, even with GPU acceleration.
To alleviate the computational complexity substitute models are used. One such model used for EM simulations including implants is the Huygens' box. 18,19 This method takes a two-step approach to compute the RF fields. First, a simulation without an implant is done where the RF fields are computed and used to construct a box around the implant. Surface currents running on the Huygens' box are computed that create the same incident RF field inside this box. After this, the implant is placed inside the Huygens' box and the surface currents found are applied, on its boundaries, to the simulation which results in the total RF fields. Everything outside of the Huygens' box is ignored in this second simulation resulting in a smaller computational domain. A challenge with this is making the Huygens' box large enough such that the scattered RF fields created due to the implant are not reflected back into the box by something that is outside of the box (i.e. there should be no crosstalk between the two domains).
Another substitute model that is applicable to elongated implants, e.g. pacemaker leads, is the electric field transfer function (TF). This transfer function describes the electric field enhancement at the tip of an elongated implant for a given incident tangential electric field exposure, 20 where the incident tangential electric field is acquired by an FDTD simulation without the implant geometry present. This effectively entails that the scattered RF field due to the implant is superimposed on the incident field, thereby decoupling the concurrent simulations of transmit coil, human model and implant into concurrent simulations with only the transmit coil and human model.
The use of a TF drastically decreases the number of fullwave simulations that need to be performed. However, as mentioned before, the transfer function is only valid for elongated implants, which is a subset of a large number of different possible implants. Furthermore, the TF can only predict the electric field enhancement at the tip of the implant. The idea of the transfer function was extended to a transfer matrix in the work of J. P. Tokaya et al. 21 The transfer matrix can predict the electric field enhancement at any location along the elongated implant, rather than only at the tip.
Although the TF enables quick calculation of the RF field enhancement due to an elongated implant, its use comes at the price of a loss of accuracy compared to a full-wave simulation. This was shown by E. Cabot et al 22 where it was found that there is a difference (up to 48%) between the concurrent full-wave simulation of the implant and the patient model compared to computing the response of the implant by the use of the TF.
Due to the aforementioned problems with the current methods, very long simulation times or sub-optimal accuracy of substitute models, there is a need for a new and more efficient method. Therefore, in this work, we will investigate a fast and accurate generalized methodology to do RF safety assessment for arbitrary implant geometries. This is derived from the work of J. van Gemert et al 23 that describes a method for efficiently computing the scattered RF fields produced by dielectric pads. Here we use the same methodology for medical implants, thus, instead of computing the EM response of an object near the patient we are interested in the EM response of an object (i.e. an implant) within the patient.
In this method, the RF fields are simulated without the implant present and the scattered RF field produced by the implant is computed and afterward superimposed onto the incident field. The computation to include the EM response of the implant is achieved through a small domain inversion, using the Sherman-Morrison-Woodbury formula. 24 We assume that the matrix that needs to be inverted is non-singular, which is normally satisfied. 23 The inversion is computed on a much smaller domain than the initial simulations. Therefore, the effect of the implant can be computed almost instantly. Furthermore, since the simulation with the source and patient is decoupled from the implant, the electrical properties, geometry, and location of the implant can easily be altered without doing another full-wave simulation, making this an efficient method for a tier 4 15 safety assessment.
To compute the scattered RF field due to an implant, a library and the RF fields without the implant are required. The library consists of the EM response for a unitary current density for each location (i.e. voxel edge) that can be occupied by an implant, which can, for example, be simulated using the FDTD method. Computing this library is a one-time effort and once available it facilitates computing the effect of different materials (i.e. electrical properties) and locations of these implant geometries within the patient can be computed in an extremely fast manner. This decreased computational effort allows for efficient evaluation of the RF safety assessment for implants.
In comparison to existing full FDTD simulations, the presented method achieves unprecedented acceleration factors. This may enable RF safety assessment of implants at much lower costs, may benefit the design of implants, and ultimately may even enable online RF safety assessment of implants.

| THEORY
We follow similar steps as 23 and start with the Maxwell equations, given by here H is the magnetic field, E is the electric field, σ is the conductivity, ɛ is the permittivity, ω is the angular frequency at the Larmor frequency, μ is the magnetic permeability, and J ext is the external current density, i.e. the current running through the RF coil. In an MRI environment, all materials have a magnetic permeability of μ = μ 0 . Equations 1a and 1b are defined on a continuous domain. For numerical analysis the domain is typically discretized into a voxelized grid. The discretization of Equations 1a and 1b can be written in matrix vector notation as where D contains the curl operators and C bg , the electromagnetic properties matrix, contains the electrical properties and is defined as C bg = diag(c bg ), where the vector c bg is written as where bg is used as shorthand notation for background, indicating that there is no implant present. The subscripts k and l indicate the number of edges and faces of the discretized domain respectively. The vector f bg contains the E and H fields and q contains the external current densities. Equation 3 can be written more compact as where A = (D + C bg ). Solving the field distributions for a given external current density is performed through the inversion of A It should be noted that A encompasses the complete simulation domain, which can be dozens of millions of voxel edges for realistic situations. Therefore, this inversion is not feasible and the fields can only be computed using numerical methods (e.g. FDTD or FEM).
However, we are now interested in a small perturbation in this A matrix created by a change in the dielectric properties, for example, due to an implant.
If we were to add this implant to the simulation domain and keep the discretization the same, we would need to solve where C imp contains the change in electrical properties between the newly added implant and the background. Similar to where the superscript imp is used as shorthand notation for implant. This operation is equivalent to deleting the electrical properties of the background and adding those of the implant. Since the magnetic permeability of objects in an MRI should be, approximately, equal to μ 0 the bottom half of the vector in Equation 8 is equal to zero. Furthermore, at locations where the implant is not present the change in conductivity and permittivity is zero too. Therefore, the change in the medium property matrix is confined to a very small (low-rank) domain within the matrix A. This small domain consists of M edges whereas the entire domain on which A is defined has N edges.
To map quantities from this large domain to the small domain the support matrix S is introduced. The . and solving for the field distributions results in Solving this still requires an inverse operation on the large domain. This is not possible for realistic simulation domains due to the number of edges in the simulation domain. However, there is a matrix identity that allows us to rewrite this equation to our advantage. This is called the Sherman-Morrison-Woodbury matrix identity 24

and is given for Equation 12 by
where I M is an M by M identity matrix. We will now introduce a new matrix called the library matrix Z This matrix is an N by M matrix where every column is the field response for a unitary current density of the corresponding edge in the support matrix S. This matrix needs to be simulated before computing the response of any implant. Building the library matrix poses an extensive simulation effort, M numerical simulations need to be performed. However, each separate simulation converges quickly because there is only a single edge source present. The library only needs to be computed once for a given dielectric background environment (e.g. for a specific patient model). After this, the response of any arbitrary implant can be calculated almost instantly. Now substituting Equations 6 and 14 into Equation 13, we obtain Note that the inverse in this equation only needs to be computed on the small domain. This allows the computation of the E and H to be possible with this methodology. From Equation 15, two key points can be observed. The first is that the total field is the sum of the two different RF fields, the incident fields and a field that is dependent on the scatterer (i.e. the implant). The second is that a generalized form of the transfer matrix 21 is defined within this equation. This can be seen when realizing that the library matrix, Z, has to be multiplied by the scattered current density within the implant. Furthermore, we can observe from Equation 1a that the conductivity and permittivity, the quantities we are changing, are only multiplied with the electric field within f bg . Therefore, the terms in between Z and f bg must be equivalent to the transfer matrix. This point is explained in more detail in Appendix A. The generalized transfer matrix can help provide insight into what implant characteristics significantly influence the scattered RF field. The above-mentioned key points show another way of looking at how and why this lowrank inverse computation works and more specifically which electromagnetic quantities affect the total RF field.

| METHODS
In order to compute the scattered RF field created by an implant using the presented method, a simulation with the transmit coil and patient model is required (i.e. the implant is not present). The RF field computed in this simulation represents the incident/background RF field, f bg . Further, the library matrix, Z, needs to be computed. The columns of the library matrix represent the RF fields on the edges in the large domain for a unitary current density, J = 1A∕m 2 . All the edges that can be occupied by the implant need to be simulated. Therefore, constructing the library matrix requires a considerably large set of simulations. The described inputs have been computed using an FDTD software package (Sim4Life, ZMT, Zurich, Switzerland).
To validate the method, a separate simulation is performed with the transmit coil, patient model, and implant present. This simulation will produce the total electric and magnetic fields, f, which are compared to the total fields produced by the computation performed with the presented method.
Three different implant structures are used to test the method. The first represents an orthopaedic surgical implant: a metallic screw. The geometry of the orthopaedic screw is shown in Figure 2A, while the location with respect to Duke and the birdcage coil in the transverse and sagittal plane are shown in Figure 2B,C, respectively. The second implant resembles a deep brain stimulator (DBS) lead structure and the third is a DBS lead structure that is tilted with respect to the FDTD grid axes. Both of these DBS lead structures are shown in Figure 3. The FDTD simulations for the passive implant are simulated at 128 MHz (3T) and for the DBS implants the simulations were done at 298 MHz (7T). For all implant types, the convergence level of the simulation with and without implant was set at −50 dB, while the library matrix simulation had a convergence level of −30 dB.
The simulations for all three implants were calculated on a 1mm isotropic grid. The orthopaedic screw was simulated with a conductivity of 2.38 ⋅ 10 6 S/m, the conductivity of titanium, and a relative permittivity of 1. The DBS electrode consists of two different materials, a conductive lead and an insulation layer around the lead. For the lead 2.38 ⋅ 10 6 S/m, and r = 1 was chosen, while for the insulation material the electrical properties were chosen to be σ = 0 S/m and r = 4.
The computations for the FDTD simulations were calculated using a GPU, NVIDIA TITAN X. The update was performed with the Julia programming language 25 using a CPU, Intel Xeon E3-1270 v5 (@ 3.60 GHz), and 64GB of RAM available. To solve the inverse in Equation 15, the generalized minimal residual method (GMRES) was used.
As a sanity check that GMRES finds the correct solution, we look at the physical interpretation of the solution of the F I G U R E 4 Comparison of the electric field components obtained by FDTD and the proposed inverse computation method from a surgical screw. The three rows show the magnitude of the E x, y, z components respectively. The first column shows the magnitude of the electric field if there is no implant present. The second column shows the electric fields with the implant present computed by the FDTD method. For the same implant, the third column shows the output of the computations performed with the presented methodology. The last column shows the error percentage as computed by Equation 18 system Ax = b in our case. As shown in Appendix A, the solution of our system, x, is the scattered current density given by However, we can also write the scattered current density, per definition, as when written in this form, the scattered current density can be computed using quantities from the FDTD simulations for the incident and total RF fields.

| RESULTS
To validate the presented method, we compare the results from Equation 15 with the simulation from the FDTD solver when the implant is present. For the orthopaedic screw, the results are shown in Figures 4 and 5 for the E and H fields respectively. In both figures, the magnitude of the x, y, z components of the fields is shown. Furthermore, error plots are shown where we defined the error between FDTD fields and the fields as computed by Equation 15 as

F I G U R E 6
Comparison between the RF fields computed with the FDTD and the presented method for the straight deep brain stimulator lead (aligned with grid axes). On the left, the location of the computed domain within the model is indicated with a red contour. The top row of figures shows the magnitude of the electric field for the FDTD simulation, the inverse computation and the error percentage as computed by Equation 18. Equivalent plots are shown for the magnetic field in the bottom row where f FDTD and f inv are substituted for the x, y, z components of the E and H fields, f FDTD are the fields obtained from the FDTD solver, whereas f inv denotes the RF fields obtained from the inverse computation. The error is scaled by the maximum value in the field, rather than the local field value. This was performed to suppress errors in regions where the magnitude of the fields are very small (e.g. inside the implant). Otherwise, these small deviations inside the implant would result in large error values although they are of minor concern because the peak values in the electric field contribute significantly more to the heating of the tissue. The ratio between the peak value of the electric field and the electric field inside the implant is a few orders of magnitude. In Table 1, the maximum errors as  computed by Equation 18 are shown. For the DBS electrodes, the magnitude of the E and H fields are shown in Figures 6 and 7. Again the maximum errors, as defined by Equation 18, are shown in Table 1. Between the three different implants shown, we find that the range of the maximum errors is given by 0.04% to 1.35%.
In Table 2 domain for the FDTD simulation is given. Furthermore, the number of edges that the implants consist of is shown. This determines both the dimensions of the square matrix that needs to be inverted according to Equation 15 and the number of columns of the library matrix. The computation time, i.e. on the GPU, per column for the library matrix and the total computation time are also given. Finally, the computation times for both methodologies are given together with the acceleration factor. The latter is defined as with Acc as the acceleration factor, t FDTD as the computation time for the FDTD simulation (using either the CPU or GPU) and t inv for the proposed inverse updating method (CPU based), without the calculation of the library matrix and incident field included. The acceleration factor that is found between the two methods is between 2478 and 171 times faster for the proposed method. This acceleration in simulation time entails that the break even point (BEP) of simulations, meaning that using the proposed method with its corresponding precomputation step is as fast as FDTD, when 22 and 55 simulations are done for the case of the first and second implant, respectively. When more implant geometries/locations with varies incident field exposures are required, which for implant safety assessment standards is certainly the case, the proposed method is faster than FDTD. The BEP is calculated as, Finally, to show that the minimization process finds the correct solution of the system the scattered current density is computed for both implant geometries according to Equations 16 and 17. The current density is summed for all the transverse slices (xy-plane) of the implant to make the plots readable. The result is shown in Figure 8, where it is clearly seen that the minimization process finds the correct solution, i.e. the blue and black line are directly on top of each other and the difference between them is two orders of magnitude smaller than the actual magnitude of the current density.

| DISCUSSION
This work has demonstrated an alternative approach to calculate the RF field response of a medical implant in an MRI. As an input, the method requires the incident RF field distribution (RF field without an implant present) and a library consisting of the unitary current density response of every voxel edge on the discretized implant geometry. To demonstrate the validity of the method, the method is tested for a screw and a deep brain stimulator lead where the input fields are determined by FDTD simulations.
From the maximum errors shown in Table 1, it is clear that this methodology is very accurate. The accuracy is only subject to the numerical precision of the supplied incident and library fields. This is further substantiated by the results shown in Figures 4-8, where the RF fields and scattered current densities computed by the presented method are shown to be equivalent to those computed by the FDTD method.
One major difference between the presented method and the Huygens' box is that the reduced domain in the presented method is only as large as the implant itself, whereas with the Huygens' box the reduced domain should be large enough that there is no crosstalk between the full and reduced domain. Due to the nature of the inverse computational complexity, (M 3 ), the acceleration factor for this method is dependent on the number of edges that the implant occupies. Therefore, the larger the number of edges the implant occupies the longer the simulation time becomes. This occurs when either the implant size is increased or if the discretization is performed on a finer grid. This is also shown in 23 and can be observed by Table 2. The computation time of the inverse, however, is independent of the frequency of the RF fields and the voxel size, i.e. geometric resolution, while FDTD simulations are dependent on these properties. This means that very small implants on a very fine grid would require a precomputation step, i.e. computing the library matrix and incident fields, that is slower while computing the EM response of the implant will be equivalently fast for a similar number of edges that need to be updated.
Another, potentially more restricting factor is the memory requirement. The library matrix grows linearly with the number of edges. For the orthopaedic implant given here, the library matrix is already 4GB for the electric fields (and another 4GB for the magnetic fields, however only the electric fields are needed for the computation). On top of this, the memory requirements for the inversion that needs to be computed grows with the square of the number of edges the F I G U R E 8 Comparison between the true solution, as computed by Equation 17, and the solution found by the inverse computation, as defined by Equation 16. The current density is summed for the transverse (xy-plane) slices. The top row shows the result for the orthopaedic implant and the second row shows the result for the straight DBS implant and the bottom row shows the result for the tilted DBS lead implant occupies, i.e. 0.1 GB for the orthopaedic implant in this work and 0.6 GB for the DBS electrode. Therefore, while theoretically possible, in the current state of the presented method, it would be difficult to compute the response of highly realistic lead structures. Both due to the large structure of such an implant and the high resolution required to capture all the details, i.e. the helical lead structure. This would increase the simulation time for the incident field and the library matrix. The resulting matrices required for our method would become too large, both M and N grow cubed with the factor increase in resolution. The memory requirements of the library matrix and the inverse scale with N by M and M by M, respectively, i.e. with the 6 th power of the factor increase in resolution. We are currently investigating ways to decrease the memory requirements for the presented method. Some of the ideas are discussed below. The current setup and implementation of the presented method would serve well for orthopaedic implants which usually are not tested for RF safety and are either smaller in size or can be modelled on a coarser grid.
To tackle the previously mentioned memory problems, we could approximate the library matrix, Z, by exploiting two properties to introduce sparsity into the library matrix. First, the presented method involves the simulation of a full library matrix, while simulations of current density sources that are spatially located near each other have very similar EM responses due to the equivalent dielectric surrounding. Second, the magnitude of the RF fields decays very rapidly for increasing distances away from the source location. This implicates that the value of the current density at any edge of the implant is dominated by the edges that are located close to it. By either interpolating between columns of the library matrix or truncation of the data if the magnitude becomes too small, sparsity can be introduced into the library matrix at the cost of the accuracy of the computation. These and other alterations for improved performance will be investigated in subsequent studies.
Assuming that the limitations described above can be addressed sufficiently, the presented method bears strong potential for applications in RF safety assessment of implants in MRI, since the calculation time of the RF fields is now in the order of seconds. One example is the safety assessment of implants according to the ISO/TS 10974 technical specification. The output of this technical specification is a conditional label for the implant that specifies the maximum B +, rms 1 and/or other RF power-related settings a patient with the particular implant can safely undergo an MRI examination. For the most rigorous RF safety assessment level (tier 4) of this technical specification concurrent simulations of the implant, patient and transmit coil are required for a wide variety of potential implant locations and trajectories. Although this method will result in the least restrictive scanning constraints, it is often considered too demanding. With the presented method, the field response of every voxel edge in the domain only needs to be calculated once after which the RF field distribution for any potential lead wire trajectory can be computed almost instantly. This may greatly reduce the workload for tier 4 safety evaluations of implants, given that we have access to a library of different RF field exposures and the libraries of different human models.
Another application could be to predict the local RF field enhancement prior to MRI examination of the patient. The implant structure and location could be revealed by the help of previously acquired X-ray photos of the patient. After the implant is localized a quick RF field calculation could be performed based on pre-calculated RF field distributions, both for the incident fields and the library matrix, using generic body models. This calculation would result in a situationspecific power threshold by which the overestimation is reduced to a minimum. This could possibly be achieved for implants without a conditional label, the RF safety assessments could be performed beforehand to verify if a patient with such an implant can safely undergo an MRI exam.

| CONCLUSION
In this work, we have shown a new methodology for RF safety assessment of implants in an MRI setting without assumptions on the implant geometry or composition. With appropriate simulations done beforehand, the presented method can perform the RF safety assessment in a greatly accelerated fashion compared to full-wave simulations, e.g. FDTD.
The incident fields when no implant is present and a library matrix, containing the EM response of every edge the implant can possibly occupy, need to be precomputed. Afterward, the effect of any arbitrarily shaped and positioned implant, with arbitrary material properties, can be calculated within seconds. The result of the computation is numerically equivalent to the solution of a full-wave simulation. For the implants shown in this work, the maximum error was 1.35%. However, using this method a significant acceleration is obtained (a factor 171 to 2478 compared to GPU accelerated FDTD simulations). This is excluding the calculation of the library matrix and the incident RF field.