In silico evaluation of the thermal stress induced by MRI switched gradient fields in patients with metallic hip implant

This work focuses on the in silico evaluation of the energy deposed by MRI switched gradient fields in bulk metallic implants and the consequent temperature increase in the surrounding tissues. An original computational strategy, based on the subdivision of the gradient coil switching sequences into sub-signals and on the time-harmonic electromagnetic field solution, allows to realistically simulate the evolution of the phenomena produced by the gradient coils fed according to any MRI sequence. Then, Pennes’ bioheat equation is solved through a Douglas–Gunn time split scheme to compute the time-dependent temperature increase. The procedure is validated by comparison with laboratory results, using a component of a realistic hip implant embedded within a phantom, obtaining an agreement on the temperature increase better than 5%, lower than the overall measurement uncertainty. The heating generated inside the body of a patient with a unilateral hip implant when undergoing an Echo-Planar Imaging (EPI) MRI sequence is evaluated and the role of the parameters affecting the thermal results (body position, coil performing the frequency encoding, effects of thermoregulation) is discussed. The results show that the gradient coils can generate local increases of temperature up to some kelvin when acting without radiofrequency excitation. Hence, their contribution in general should not be disregarded when evaluating patients’ safety.

put in evidence hot-spots close to the tip of the stem and the screw of the acetabular cup (that was not present in the model used by Mohsin) and, for a local SAR (averaged over 10 g) below 10 W/kg, found a maximum final temperature exceeding the limit prescribed by Standard IEC 60601-2-33 (2010), i.e. 39 °C.
Unlike RF effects, limited effort has been devoted to the thermal effects due to the switching gradient coils (GCs) fields, mainly analyzed for exertion of forces, induced voltages and imaging artefacts in the presence of implants (see for example Erhardt et al (2018) and Kalb et al (2018)). However, under realistic operating conditions, GCs may produce a significant amount of power inside a bulky metallic object, as proved by preliminary theoretical estimates on a model problem (a sphere in Zilberti et al (2017)). Experimental corfirmations are given by Graf et al (2007), which measured a temperature increase up to 2.2 K in an aluminum replica of a hip implant after 210 s of exposure to a 3D true FISP (fast imaging with steady precession) sequence, and by Bruehl et al (2017) on an acetabular cup, as will be discussed in section 3. The role of GCs fields in the increase of tissue temperature in patients with joint arthroplasties is verified also by simulations. In El Bannan et al (2013) the numerical and experimental analysis of twelve rods of different metals inserted in a 1 kHz driven solenoid shows a temperature rise up to 2.45 K. Zilberti et al (2014) consider an anatomical human model with bilateral hip implant placed in a MRI combined with a linear accelerator and compute maximum temperature elevations up to 2.7 K under 1 kHz supply for GCs. The calculation of a maximum temperature rise of about 4 K is described in Zilberti et al (2015), for the case of a patient with hip implant during 30 min of exposure to trapezoidal GC signals produced by a conventional system of coils working with a 20% duty-cycle. These previous analyses put in evidence a potential risk of gradient-induced heating of implants, but are difficult to compare, because of the different features (spatial distribution, waveform, frequency and magnitude of the applied field; shape and material of the test object; duration of the exposure) of each analyzed situation. In addition, they all present some important approximations (for instance, use of a test object with simplified shape instead of a real implant, or use of uniform and sinusoidal magnetic fields, instead of fields with realistic spatial distribution and time waveform) which make the results not completely suitable to describe a real scenario.
Differently from the RF heating, where a significant amount of energy deposition occurs directly in the biological tissues, almost the whole GCs thermal energy is deposed in the metallic implant and then diffuses towards body tissues surrounding it, giving rise to a local phenomenon. The effects of GCs are usually stronger when the implant is placed in the scanner periphery. There, as can be seen in figure 1, the concomitant transversal components B x and B y produced by the GCs, whose presence represents a minor source of distortion for standard MRI, can become stronger than the longitudinal component B z for some GCs (Liu et al 2003). Thus, the implant positioning within the scanner becomes an essential parameter.
Up to now, the computational papers devoted to this subject have approximated the supply conditions by adopting periodic sinusoidal or trapezoidal waveforms, which are quite far from the much more complex waveforms of GC fields nowadays applied in MRI. The rare papers which use actual sequences, e.g. 3D true FISP in El Bannan et al (2013) and Arduino et al (2017), apply a simplified approach, not suitable to take into account the details of many actual sequences (e.g. when the signal level is modified at each repetition). This paper proposes a novel strategy specifically developed to realistically account for the time behavior of GCs fields during any MRI sequence in the electromagnetic and thermal simulations required to estimate the temperature increase. The technique is based on the decomposition of the supply waveforms into sub-signals. It leads to the solution of a limited number of electromagnetic problems in the frequency domain. The power density generated by the GCs fields, increased by the one generated by the RF antenna, becomes the forcing term of a transient thermal problem, which provides the distribution of the temperature increase inside the patient body. In order to focus the discussion on the effects due to the exposure to the GCs switching fields, in the following the power generated by the RF field is discarded.
The computational approach used for the in silico evaluation is described in details in section 2, related to Methods. Then, the proposed procedure is applied to the analysis of the GCs induced heating of a hip implant due to a echo-planar imaging (EPI) sequence, as detailed in section 3. This sequence was chosen for two main reasons: (1) it is quite aggressive from the point of view of energy deposition; and (2) it has characteristics that make it representative of other types of sequences from the viewpoint of the approach here proposed. In addition, the adopted EPI sequence is similar to the sequence that was used by Bruehl et al (2017) to perform an experiment within a real scanner working in 'normal operating mode' (IEC 2010), hence considered safe from the viewpoint of cardiac and peripheral nerve stimulation.
In the analysis, the effects of the body positioning with respect to the scanner isocenter and the choice of slice selection, phase encoding and frequency encoding directions on the GCs induced implant heating are discussed.
It must be remarked that the present work investigates the heating due to the only GCs, disregarding RF. Thus, it cannot give a full response about the admissibilty of a patient to a MRI scan, for which it gives a condicio sine qua non, anyway. In addition, while RF SAR may involve a large portion of the body, producing not only a local heating, but, sometimes, a core temperature increase (van den Brink 2019), gradient-induced heating requires the presence of a metallic implant and keeps confined around it. For this reason, the following analysis focuses on local thermal effects.

Methods
Since the expected temperature increase, in the order of a few kelvin, is small enough to not alter sensibly the electric properties of the implants' materials and the biological tissues (according to Trujillo et al (2013), the relative variation is limited to few percents), the proposed numerical procedure can be divided into two successive separate steps. A set of EMF solutions provides the instantaneous power deposed in the orthopedic implant; then, a thermal problem, starting from the previous result, describes the consequent heat diffusion and temperature increase in the surrounding tissues.

Electromagnetic problem
Due to the extremely low electric conductivity of the human body with respect to metal, the electromagnetic problem is developed only inside the medical implant, under the reasonable assumption that, for the involved frequencies, the currents induced into the body tissues neither generate a significant thermal power, nor modify the magnetic field produced by GCs. The latter assumption would introduce a relative error on the induced electric field within the body tissues lower than 10 −4 (estimated using the analytical solution reported in Zilberti et al (2017) at the frequency of 100 kHz). As an original alternative to a step-by-step procedure, the proposed approach is based on time-harmonic EMF solutions and is structured in successive phases. The related operations are illustrated in the following, making reference to an EPI sequence, whose waveforms, reported in figure 2, are rather complex and challenging to be simulated. The EPI gradient waveforms have been recorded from nominal-current monitor of the gradient amplifier of a Siemens Verio 3 T scanner at PTB. The vendor supplied EPI sequence was used and adjusted to achieve maximum heating. Shim currents and noise have been removed.
The proposed procedure can be divided into three steps: signal pre-processing, electromagnetic simulations, power signal synthesis by post-processing.
The former step consists in subdividing the waveform supplied to each GC during a given time interval Δ (e.g. the time frame of the EPI sequence, or the repetition time TR for non-single-shot sequences) into sub-signals, which can be either periodic or aperiodic. Each sub-signal is then represented through a truncated Fourier Figure 1. Spatial distribution of the GC field components (mT) along the mid-plane (plane x-z) of a tubular scanner for a rated gradient of 20 mT/m. From left to right, the magnetic flux density components (B x , B y and B z ) generated by GC for gradient-X, gradient-Y and gradient-Z (from upper to lower, respectively). The B z distribution generated by the gradient-Y coil in the plane y -z is identical to the one generated by the gradient-X coil in the plane x-z and, in the same way, the B x and B y components are identical, but exchanged. Vice versa, for the gradient-X coil in plane y -z. expansion via fast Fourier transform. To determine the level of the truncation, an estimate of the relative error made with respect to the deposed energy is provided by the error index ε n defined as being Bñ the sub-signal waveform approximated by the Fourier series truncated at the n-th harmonic order, B the original sub-signal waveform, and Δ s the duration of the considered sub-signal. Precisely, the Fourier series is truncated when the error index ε n is less than 5%. The integrals that appear in (1) are proportional, with the same proportionality factor, to the energy that would be deposed by the EMF in the radiated object if the skin effect (i.e. the confinment of the induced currents at the periphery of the object) were negligible; thus, they can be used for a preliminary comparison of the thermal effects produced by the actual and the truncated waveforms. Once the sub-signals are identified, the second step of the proposed methodology requires, for each harmonic of each sub-signal, the solution inside the implant of an EMF problem with unitary source. It is important to underline that, when more sub-signals in the same coil have the same fundamental frequency, the field solution is performed just once for that frequency and the related harmonics. The eddy currents problem is formulated through an electric vector potential and a magnetic scalar potential (T − φ formulation), where the magnetic field distribution produced by the considered GC is assumed as the driving term . A hybrid Finite Element/Boundary Element method, solved by a generalized minimal residual (GMRES) algorithm running on graphics processing units (GPUs), provides the EMF distributions (Bottauscio et al 2015a). The adopted homemade implementation of the electromagnetic solver has already been applied and validated for the analysis of RF problems (Bottauscio et al 2015b). Details of the numerical implementation are reported in the appendix. The value of the computed electric field E in each implant voxel is stored.
Finally, in the last step of the procedure, the time interval Δ is sampled. For each voxel, the stored electric field harmonics are scaled by the complex Fourier coefficients of the related sub-signal. It must be noted that the change of one or more signal amplitudes in successive repetition times (e.g. in the phase encode of Gradient Echo sequences) involves the update of the scale factors only, without requiring additional field computations. The electric fields are then moved back in the time domain, located in the correct time interval and superimposed coil by coil to reconstruct the instantaneous distribution of the total electric field in all time samples. Finally, the instantaneous Joule power density and the energy density deposed in each voxel during each time step are computed. At the end of the process, in each implant voxel the power is averaged on the time step of the marching method to solve the thermal problem, inserting possible idle times between two successive intervals Δ. At this point, the power contribution due to the radiofrequency (RF) coil, provided by another numerical code, could be added in all body voxels in case of a complete simulation of the MRI session.
As an example, the process based on the harmonic decomposition, applied to the gradient-Z coil waveform, is presented in figure 3, where sub-signals 1, 2 and 3 are colored in red, blue and green, respectively. It is worth noting that 'idle' intervals must be inserted before and after the aperiodic sub-signals, in order to make the Fourier approximation reasonable for both magnetic and induced electric fields (i.e. to take into account possible transient evolutions when superposing the effects of different sub-signals). Moreover, by virtue of the arbitrariness in the choice of the time window covering the aperiodic signals, the fundamental frequency can be suitably 'tuned' to reduce the number of required electromagnetic simulations. Indeed, on the one hand, the same duration can be adopted for different sub-signals in the same coil, reducing the number of time-harmonic problems to be simulated numerically. On the other hand, in some cases idle times can be adjusted to obtain a waveform whose symmetry rules out the even harmonics (apart from the DC component, which is irrelevant in the process of electromagnetic induction, anyway). In figure 3, the Fourier expansion is truncated to the 63rd harmonic for the two aperiodic sub-signals (using the same fundamental frequency) and to the 15th harmonic (with only odd components) for the periodic sub-signal. Thus, only 71 simulations are required to evaluate the energetic effect of the gradient-Z coil effectively. The corresponding values of the error index ε n are reported in table 1.
The reliability of this approch for aperiodic signals has been tested on the model problem consisting of a sphere (diameter 40 mm, comparable with the femoral head of a hip implant) made of the alloys usually adopted for hip implants with the higher and the lower electric conductivity (stainless steel and Ti6Al4V, respectively, see table 2). A trapezoidal waveform for the magnetic field is imposed to simulate one of the aperiodic sub-signals (e.g. the last sub-signal in the gradient-Z coil of the EPI sequence). The waveform of the power density p (t) generated in an inner point close to the sphere surface is computed using the proposed time-harmonic approach with a decomposition truncated at the 31st component. The result is compared with the one given by a time-domain 2D homemade solver, suitable for general purpose electromagnetic simulations. The comparison, presented in figure 4, shows a good agreement between the two solutions. In particular, the decay time of the power density p (t), during the plateaus of the magnetic field waveform, are accurately reconstructed. About this aspect, it is worth noting that the reported power density is proportional to the square of the local electric field, which does not simply reflect the time derivative of the applied magnetic field, but exhibits the delay typical of a ohmic circuit with non-negligible self-inductance. This occurs because the eddy currents in the metallic objects are able to perturb the local distribution of the magnetic field itself.

Thermal problem
The time evolution of the temperature T inside the biological matter can be modelled through Pennes' equation (Pennes 1948): where ρc p is the volumetric heat capacity, λ is the thermal conductivity and h b is the blood perfusion coefficient in the human tissues, T b is the temperature of blood, P met is the volume power density associated with the metabolic process and P em is the volume power density deposed by the radiation and derived by the previous electromagnetic solution. The same equation can be used to describe heat transfer also within the implant materials simply by putting h b and P met equal to zero. On the body surface, the thermal field T satisfies the Robin condition: being h amb the heat exchange coefficient, n the outward normal direction and T air the unperturbed air temperature of the external environment. The thermal problem can be solved starting from the knowledge of the initial spatial distribution of temperature within the body before the exposure to the EMF (T 0 ). The thermoregulation affects both the blood perfusion coefficient h b and metabolic heating P met . According to the model proposed in (Laakso et al 2011), coefficient h b is modified by a local temperature-dependent multiplier L B = 2 (T -T0)/ΔB , assuming ΔB = 1.6 K. Thus, the perfusion term in equation (2) becomes: where h b0 is the blood perfusion at T 0 . Factor L B is saturated to 32 for skin and to 15 for all other tissues. Similarly, the metabolic heat production is assumed to be dependent on the local tissue temperature through a multiplier L M = 1.1 (T−T0) , as in Bernardi et al (2003), so that the metabolic heat term modifies as: where P met,0 is the metabolic power density at T 0 . Other non-linear contributions due to thermoregulatory processes, like sweating, that would affect the heat exchange coefficient h amb (Bernardi et al 2003), are neglected.
Since the temperature elevation produced by the EMFs is the quantity of interest, the bioheat equation can be conveniently simplified by introducing the temperature elevation ϑ with respect to the local temperature T 0 before the exposure (ϑ = T − T 0 ), following the approach proposed in Arduino et al (2017). The thermal equation before the exposure (when P em = 0): Figure 3. Sub-signals (from 1 to 3, starting from the top) deriving from the G z signal (reported on top) and related harmonic decomposition. Fourier series expansion truncated to the 63rd harmonic for the aperiodic sub-signals and to the 15th harmonic for the periodic one. In the plots of the harmonic decomposition for the aperiodic signals (1 and 3), only the first 31 harmonics are shown (but 63 harmonics are needed to obtain a satisfactory error index).
is subtracted from equation (2), leading to: In the inner tissues, the temperature at rest T 0 can be approximated with the blood temperature T b . By using the same approximation in the whole body, equation (7) can be rewritten in terms of the temperature elevation θ, as: with Robin boundary conditions λ ∂ϑ/∂n| ∂V = −h amb ϑ. The adopted approximation is quite large for the skin, which is colder than the rest of the body and plays a significant role in the whole body thermoregulation. However, in the considered simulations the heat source is the metallic implant located deeply within the body, and the obtained temperature maps are completely controlled by local parameters (mainly the blood perfusion, that dissipates the heat before it reaches the body surface). This fact makes the presented results trustworthy, despite the approximations in the mathematical modelling and sweating has been neglected. The form of Density kg/m 3 8445 4420 7900 Figure 4. Example of the reconstruction of the induced power density in a point in close proximity to surface of the metallic sphere when an aperiodic behaviour of magnetic field is assumed. The upper plot shows the aperiodic trapezoidal waveform of the applied magnetic field. The middle and lower plots report the waveforms of p (t) for a sphere made of Stainless steel or Ti6Al4V alloy, respectively. In both diagrams the waveform of p (t) computed with the harmonic decomposition is compared with the reference one given by a time-domain electromagnetic solver. The power density values are normalized to the maximum peak value (in the Stainless steel sphere) reached during the considered time interval.
Pennes' equation (8) avoids computing the initial temperature distribution within the body, determined by the equilibrium between metabolic heat, diffusion and perfusion phenomena, and heat exchange with the air. Homemade code is employed for the thermal computations. To solve equation (8), the thermal problem is developed in a domain involving the implant and the surrounding biological tissues. To reduce the computational burden, instead of simulating the whole body it is possible to select a portion of it, around the implant (see figures 7 and 9). The suitability of the adopted domain can be verified a posteriori, checking that the tempearture elevation does not extend on the artificial cuts that have been introduced to bound the domain itself. The adoption of a structured Cartesian mesh to discretize the domain (as in voxelized human models) allows the application of the finite difference method (FDM) with the Douglas-Gunn (DG) time split implemented to work on GPUs, whose efficiency has been previously verified .
In the time marching algorithm, the correcting factor L B is evaluated according to the temperature increase value at the previous time instant, in order to deal with a linear system at each step. By keeping sufficiently small time-steps, this linearization introduces a negligible error in the solution. It is worth noting that thermal phenomena evolve much more slowly than the electromagnetic ones, acting like a low-band filter for the rapidly varying electromagnetic power density. Without thermoregulation, problem (8) is solved assuming all thermal parameters to be invariant. This last assumption provides the upper limit for the temperature elevation computed through Pennes' equation (i.e. the worst-case scenario), which corresponds, for instance, to a completely impaired thermoregulation.

Experimental validation
The homemade numerical codes which implement the procedure described in the previous section, have been validated by comparison with the experimental results published in Bruehl et al (2017). In that paper, the heating of an acetabular cup (outer diameter equal to 54 mm) of a real hip implant made of Ti6Al4V alloy was studied. The cup (see picture in figure 5(a)) was placed in a 3 T clinical scanner, programming an EPI-like sequence with continuous, trapezoidal gradient-Z. The experiment was performed with the implant thermally 'insulated' (using a polystyrene cover which surrounds the cup leaving a layer of air around it) or 'embedded' (in gelatin gel). Details of the experiment can be found in Bruehl et al (2017).
Both testing conditions have been here simulated and the comparison has been carried on making reference to the temperature increase of the sensor measuring the highest temperature increase.
In the 'embedded' implant, the simulated results are obtained assuming for the gel the following thermal properties: λ = 0.5 W/(m K), c p = 3700 J/(kg K) and ρ = 1100 kg/m 3 . For the thermally 'insulated' case, simulations have been performed assuming Robin boundary conditions at the implant surface, with heat exchange coefficient h amb between 2.5 W/(m 2 K) and 3 W/(m 2 K).
The relative uncertainty of the measured temperature increase is estimated to be around 5%, which includes the instrumental relative uncertainty and the imperfect coupling of sensor and cup. When comparing the experimental data with computations, we have to include additional uncertainties caused by the incomplete knowledge of the overall experimental parameters, such as the position of the cup within the gradient field of the MRI scanner, the physical properties of the cup (electrical and thermal conductivity, specific heat capacity) and its thickness. A global budget of all these uncertainties leads to an estimate of a standard relative uncertainty of 20%.
The maximum discrepancy observed with respect to the computed temperature increase occurs for the implant embedded in gel and it amounts to about 3.2% of the measured value; hence, significantly below the uncertainty level. The obtained good agreement between measured and computed temperature increases is supported by the trends depicted in figure 5(b).
It is worth noting that the comparison reported here validates the proposed method only partially, since the phantom that surrounds the acetabular cup cannot account for the blood perfusion and thermoregulation that affect the thermal process in a living body. However, the main novelty of the presented approach is given by the treatment of the electromagnetic problem, which is completely supported by the described experimental validation, whereas the thermal computations rely on a standard approach widely adopted in the field of electromagnetic dosimetry.

In silico simulations
The proposed procedure is applied to the evaluation of the GC thermal effects in a patient with a unilateral right hip implant, during a MRI session in a tubular scanner. The implant has a total length of 23 cm (including the screw), that is ~20 cm from the top of the femoral head to the lower tip of the stem, and includes an acetabular shell and a liner. Three different materials for the metallic parts of the implant are analyzed: CoCrMo alloy, Ti6Al4V alloy and austenitic stainless steel. They are representative of the materials usually adopted for this type of implant. The electrical and thermal properties of the metallic parts are detailed in table 2. The liner is made of polyethylene and has a negligible electric conductivity.
The hip implant has been inserted in the 'Duke' anatomical model, belonging to the virtual population (ViP) (Gosselin et al 2014), segmented into 77 biological tissues, whose electric and thermal properties are deduced from the database developed by the IT'IS Foundation (IT'IS Database 2016). The heat exchange coefficient h amb has been always assumed equal to 7 W/(m 2 K). Due to the essentially local heating, as evident in the results shown in the following, this coefficient has a weak effect on the global results.
In the simulations, the scan performed following the EPI sequence presented in figure 2, continuously applied for about 12 min in order to simulate the repeated acquisition of multiple body slices. According to the division into sub-signals shown in figure 3, eight elementary waveforms are identified: six of them are aperiodic (two for each coil) and two are periodic (one in the gradient-Y coil, driven by a fundamental frequency of ~1922 Hz and reconstructed by 63 harmonics, and one in the gradient-Z coil at ~961 Hz with 8 non-null odd harmonics). By choosing time windows with the same duration (4 ms) for all the aperiodic sub-signals, the same 250 Hz fundamental frequency allows the reconstruction of all of them with the stated error index using 63 harmonic components. In this way, the total number of EMF solutions required for an accurate simulation is limited to 260.
The same result could be obtained if the whole EPI sequence were expanded in truncated Fourier series that is assuming the time interval Δ as a period. However, in this case, the waveform is more complicate and a satisfactory reconstruction requires much more harmonic components (at least 1023 to ensure an index error lower than 5%). Consequently, the total number of EMF problems would increase up to 3069, which puts in evidence the efficiency of the proposed strategy based on sub-signals.
The computations are performed on a server with Intel Xeon CPU E5-2680 v2, 128 GB RAM, and a NVIDIA K80 GPU card. The computational time required for each EMF solution (i.e. each harmonic) is ~30 s, when the hip implant is discretized with 2 mm voxels (1250 7 voxels, ~16 000 nodal unknowns, ~32 000 edge unknowns) and a limit error of 10 −3 is set for the GMRES residual (reached after about 200 GMRES iterations). Using the same hardware, the computational time for the thermal problem (on a portion of the body involving 840 742 cubic voxels with 2 mm side and ~872 000 nodal unknowns) is ~300 s having used a 0.5 s time-step.
In the analysis we have considered the three possible imaging planes (sagittal, coronal and transversal), associating the role of slice selection to the corresponding GC and the role of phase and frequency encoding alternatively to the other two GC (the identification of the six cases is summarized in table 3).
The computations have been performed considering different realistic positions of the hip implant, where the body axis is aligned with the scanner axis. Making reference to the femoral head, the implant has been located in the 'natural' position assumed by the body on the tomograph table, that is x-coordinate ≈85 mm and y -coordinate ≈0, with respect to the MR isocenter. The z-coordinate is varied in the range from −450 mm to +300 mm.
Due to the spatial variation of the GC magnetic field (see figure 1), the heating significantly varies with position. As an example, for the simulated case #1 (sagittal imaging with frequency encoding associated to the gradient-Z coil), the maximum heating is found for a longitudinal shift equal to −300 mm, which approximately corresponds to a thoracic MRI scan ( figure 6). Predictably, this is the position where the GC magnetic field shows the highest intensity, including a significant concomitant component along the x-axis produced by the gradient-X coil. For the same GC configuration, the minimum heating is found when the femoral head is located at the level of the isocenter. The results reported in figure 6 refer to an implant made of CoCrMo alloy. A qualitatively similar behavior can be observed for the other materials. Figure 7 shows the volume power density within the CoCrMo implant (computed as the ratio between the energy density induced by the GC and the sequence duration) and the produced increase of temperature in the surrounding tissues for case #1 with the femoral head at −300 mm. For both quantities, the highest values are localized in the right side of the acetabular shell, with a peak temperature increase of about 3.2 K (without thermoregulation). The related temperature increase trend and the spatial temperature evolution are reported in figure 8. This diagram underlines that the maximum value in the whole domain is reached inside the metallic implant. Therefore, in the following, the maximum temperature increase always refers to a point in the metal. Figure 8 also puts in evidence that the heating of the implant and the surrounding tissues develops on a relatively long time scale, which, on the whole, makes the heating curve smooth (masking minor oscillations due to the waveform and the switching of the gradient signals that drive the process).
When the implant is moved to a position with zero shift (i.e. the target of the MRI exam is the pelvis, which is located at the isocenter) the peak of temperature elevation reduces to 0.08 K (figure 9), as a result of the reduction in the field amplitude and the different distribution of the power density within the implant.
The change of the imaging plane, and the consequent role of the three GC, affects the thermal stress sensibly. By comparing the different possible configurations detailed in table 3, the maximum temperature elevation ranges from 1.9 K to 3.2 K with the femoral head shifted longitudinally of −300 mm. The values corresponding to each simulated case are summarized in table 4. It is worth noting that the six cases could be grouped in three couples having the same frequency encoding direction and the same predicted maximum temperature increase. These results highlight that the choice of the direction of the frequency encoding plays a main role on the temperature increase induced by the considered EPI sequence.
As a final example, in the most severe conditions (i.e. cases #4 and #6, with frequency encoding along the x-direction and implant in Stainless steel), the maximum temperature increase reduces to around 12% moving from the position of maximum heating (thoracic MRI) to the position of minimum heating (abdominal/ pelvis MRI). The minimum heating (occuring when scanning the pelvis) for cases #4 and #6 is significantly higher than the one found for cases #1 and #5. This happens because, for small values of the longitudinal shift, all Cartesian components of the magnetic flux density produced by the gradient-Z coil (performing the frequency encoding in cases #1 and #5) are relatively weak, whereas the z-component of the field produced by the gradient-X coil (playing the role of the most 'energetic' signal in cases #4 and #6) is non-negligible at the place of the implant (which is naturally displaced laterally). Figure 10 compares the distribution of power density and z-component of the GC field, on the implant surface, for cases #1 and #4, making reference to the implant made of CoCrMo alloy. A higher value of the B-field in the acetabular region, which determines higher induced currents and deposed power, is evident for case #4.

Discussion
The results of the research highlight a potential risk in the application of the gradient fields of an EPI sequence on a bulky metallic implant, due to a possible temperature elevation which compares to that occurring, under some circumstances, owing to the RF fields and considered as a danger. About this point, it is interesting to draw a parallel between the two fields. The maximum volume power density reported in figure 7 inside the implant, due to GC fields, amounts to about 1.5 · 10 5 W/m 3 , which corresponds to about 17.8 w/kg (the density of the CoCrMo alloy is 8445 kg/m 3 ). If we imagine that such an 'equivalent specific absorption rate' drives an adiabatic heating (as happens at the first instants of the process, when diffusion has not taken place yet), given the specific heat capacity of the alloy (450 J/(kg K)), we obtain a rate of heating equal to 0.04 K/s. The SAR required to obtain the same rate in bones and muscle (where the specific heat capacity is about 1300 J/(kg K) and 3500 J/(kg K), respectively), would be 52 W/kg and 140 W/kg, which would exceed the limits on local SAR for both normal and first level controlled operating mode recommended by IEC 60601-2-33 (2010). Of course, this estimate applies in the first stage of the heating process, when diffusion has not taken off yet, and cannot provide a complete term of comparison with the SAR generated at RF because in the metallic parts of the implant, where thermal conduction is much larger than that of biological tissues, the heat diffuses faster. The computational results reported in the previous section have also put in evidence the importance of simulating the thermal stress in metallic implants by taking into account both the actual spatial distribution of the magnetic field gradient (also outside the imaging region) and the time waveforms of realistic clinical sequences. The combined effects of implant position and signal shape determines the thermal stress experienced by the patient.
The dependence of the heating on the position within the MR bore due to the GC field is more complex with respect to the largely studied effect of positioning in case of RF heating.
In particular, while RF heating is maximized when the metallic implant is located within the RF coil (imaging region), for the GC heating the most severe conditions are always found when the implant is relatively far from the MR isocenter. Such a result is justified by the intrinsic features of the GC fields, which inside the imaging region are requested to linearly increase moving away from the isocenter (along the axial z-direction for the field of gradient-Z coil or toward the radial directions for the other coils, as well illustrated by figure 1). Outside the imaging region, the spatial distribution of the field is less regular for gradient-X and -Y coils, that also show significant spurious components along x and y directions. In particular, the B x distribution generated by the gradient-X coil in figure 1 exhibits the absolute flux density peak in proximity of the side wall of the x-z plane (the same for B y generated by the gradient-Y coil in the y -z plane). Such a field concentration implies a strong power deposition by the gradient-X coil when the implant is placed in this region. The same energetic effect could be provided by the gradient-Y coil (in the y -z plane), but the limited extent of the human body in the y -direction avoids that the metallic implant and the biological tissues are subjected to the B y highest values.
Ultimately, two basic considerations can be deduced from this analysis. The first one, as already mentioned, is that the thermal stress in the implant and the surrounding tissues are higher when the implant is far from the isocenter of the MR scanner. This effect is found for all the simulated situations, independently from the imaging plane that is considered, even if the entity of the variation depends on the supply conditions. The other  conclusion is that, with this kind of GC, the gradient-X coil (because of the B x spurious component) and gradient-Z coil (because of the B z component) are responsible of the strongest thermal stress, while the contribution of the gradient-Y coil is less relevant, although non-negligible.
The other essential factor affecting the thermal behaviour is the waveform of the sequence signals flowing in the GC. The analysis of the EPI sequence signals evidences that the prevailing contribution to the thermal stress is due to the frequency encoding waveform and, in particular, to its periodic sub-signal. A simulation performed under the same conditions as in case #1 with the CoCrMo alloy, but considering only this sub-signal (all other sub-signals in the three coils are suppressed), leads to a temperature elevation of 3.1 K, to be compared with 3.2 K when the actual EPI sequence is considered. It follows that, on the basis of the previous considerations, the maximum heating is reached when the frequency encoding is imposed on the gradient-X coil or on the gradient-Z coil. On the contrary, the temperature elevation significantly reduces if this signal is carried by the gradient-Y coil (see table 4).
More in general, the waveform of the signals, and in particular the presence of high-order harmonics in the Fourier expansion, is crucial for the thermal stress. If case #1 with the CoCrMo implant is simulated assuming a pure sinusoidal waveform for the periodic sub-signal of the frequency encoding signal (having the same amplitude as the actual signal), the temperature elevation decreases from 3.2 K to 2 K. In this context, the slope of the trapezoidal waveform of the periodic sub-signal of the frequency encoding plays an important role. If, with the CoCrMo implant, the slope is increased from 167 (T/m)/s (value adopted in case #1) to 218 (T/m)/s and to 311 (T/m)/s, keeping constant both the plateau of field gradient and the signal period, the maximum temperature elevation increases from 3.2 K to 3.85 K and to 4.80 K, respectively.
The physical properties of the metallic components, and in particular the electrical conductivity, affect the heating of the implant. At the increase of the electrical conductivity, the power density dissipated in the metal increases, as well as the heating. However, the skin effect in the metallic component (which increases with the  harmonic frequencies) makes this behavior non-linear with the conductivity value. By changing the electrical conductivity from 0.58 × 10 6 S/m (Ti6Al4V alloy) to 1.25 × 10 6 S/m (Stainless steel) the maximum temper ature increase computed in case #1 varies from 2.26 K to 3.30 K, while a linear behavior would lead to 4.96 K. Finally, we point out that the effect of the thermoregulation is found to be almost negligible (the variation of the maximum temperature is less than 0.15 K for the considered cases). This relatively low impact of thermoregulation complies with results reported by other authors (Bernardi et al 2003, Kodera et al 2018 that put in evidence a notable effect only at temperature elevations higher than 2 K. In the specific case here analyzed, such a weak effect can be further explained noting that the heating produced by GC is local (around the implant) and the temperature elevation involves tissues with relatively low basal blood perfusion coefficient.

Conclusion
The paper presents a procedure for the evaluation of the thermal effects produced by the gradient coils in the body of patients with orthopedic implants during the execution of an EPI sequence without RF excitation. A specific strategy, based on the subdivision of the actual gradient waveforms into sub-signals and the EMF solution in the frequency domain, has been developed. This approach allows investigating the exposure scenario in a very realistic way, and, at the same time, limiting the computational burden with respect to approaches directly applied to the periodic gradient pattern (or the whole sequence). The large variability in the dosimetric results obtained for different body positions and operating sequences confirms the importance of taking into account the specific features of the exposure situation.
This aspect, together with the fact that, in some of the proposed examples, the maximum temperature elevation exceeds 3 K, indicates that the GC heating effects must be carefully estimated, to check the existence of conditions that could represent a safety concern. If needed, safety measures can be taken, for instance, by adopting 'less aggressive' sequences when scanning patients with implants, by changing the imaging plane, or by introducing waiting times between the application of a sequence and the successive, to allow thermal diffusion within tissues and consequent redistribution of the thermal energy. All these aspects need to be studied in details in a future work.