Implanted electrically-powered medical devices are commonly used to provide a wide range of therapeutic functions. Examples include cardiac pacing, continuous microinjections, and neurostimulation. The electrical power used to energize these devices ultimately ends up as waste heat and is absorbed into the adjoining tissue. For devices equipped with batteries of small storage capacity, recharging must be performed frequently. The recharging may be accomplished by means of a skin-surface-mounted antenna which produces a magnetic field whose lines of force intersect the wires of a coil contained within the implant. The antenna is positioned in close proximity to the implant and aligned with it. Heating occurs during the recharging period, both within the implant and the antenna.
Predictions of tissue temperature levels and potential tissue damage due to the recharging process for neuromodulation implants are set forth in  -  . For the most part, these studies were focused on comparisons among available neuromodulation devices with respect to their tissue injury potential. Investigated operating conditions included implant placement in the tissue bed, recharging rate, heating rate control algorithm, implant-antenna misalignment, and heat loss from the skin surface. It was found that the conditions that increase the likelihood of tissue damage include 1) suppression of skin surface heat loss (e.g., patient seated on an upholstered surface), 2) high rates of power dissipation in the implant and/or the antenna, and 3) absence of a control circuit which modulates power dissipation in response to elevated temperatures.
Failsafe approaches to the avoidance of tissue damage can be designed by making use of suitable control algorithms. For instance, in  , it was shown that modulations of the recharge rate based on temperature measurements can effectively limit the potential for tissue damage. On the other hand, it is also possible to modulate recharge rates in a preprogrammed manner. High rates of charging can be performed when tissue temperatures are low, and the recharge rate can be reduced as tissue temperatures rise. Such an approach has been adopted by Medtronic in a newly designed neuromodulation device (Intellis). That device operates with appropriately programmed low- and high-recharging modes, in combination with sensed tissue measurements, to shorten the total recharge duration while insuring safe tissue temperatures.
The goal of the present investigation is to make use of advanced numerical simulation techniques to predict the impact of preprogrammed charging on tissue temperatures in the near neighborhood of both the implant and the antenna. Aside from the issue of preprogramming, the present approach differs fundamentally from those used in the preceding studies. There, the rate of heat transfer at the surface of the implant was determined experimentally by the immersion of the implant in a calorimetric fluid. Here, in sharp contrast, detailed simulation-based analysis of heat generation at the component level was implemented internal to both the implant and the antenna. Another important factor taken into account in the present simulation model is the heat capacity effect associated with the mass of the implant. A final enhancement of the present approach is the use of two independent methodologies to assess the potential of thermal injury. These numerical methods and injury assessment tools have a long history of successful use such as     .
The approach taken here, while applied to a specific implantable device can be extended to other products provided appropriate geometrical and power generation changes are made in the analysis.
2. Experimental Determination of the Implant Heat Generation
The objective of the experimentation was to determine the rate of heat generation that occurs within the implant. The experimental setup included a Pyrex beaker that was used to house the implant and a calorimetric fluid. The implant was supported in the fluid bath by miniature 2-mm high insulating pedestals. The antenna was situated 20 mm below the implant, and the respective coils in the antenna and the implant were precisely aligned. This separation distance was assured by means of a precisely machined slab of extruded polystyrene foam. These components were swaddled all around, on top, and on bottom by insulation media which included 1) vacuum insulation panels (thermal conductivity equal to 20% of that of still air), 2) silica aerogel powder (conductivity 75% of air), and 3) extruded polystyrene foam (conductivity 5% greater than air). The overall thickness of the insulation assembly was about 100 mm.
The test setup was thoroughly instrumented to measure temperatures to high accuracy. A thermocouple was attached to each of the two principle faces of the implant, and three thermocouples were distributed throughout the calorimetric fluid. Corroborating measurements of the temperature of the fluid were made with a 0.05˚C thermometer traceable to NIST. The thermometer performed dual functions as a temperature sensor and as a stirrer. Thermocouples and heat flux gauges were affixed to the bottom and sides of the glass beaker. The selected thermocouples (Type E) were of a very fine diameter (0.075 mm) and of low conductivity alloys to ensure that heat conduction along their leads would not affect their readings. The rationale for measuring temperatures at the implant proper, in the calorimetric fluid, and on the bounding walls of the containment beaker is that the heat generated in the implant distributes itself to those media.
Other relevant instruments included an electronic balance (0.01 g resolution) for determining the mass of the implant and a finely resolved graduated cylinder for the volumetric measurement of the calorimetric fluid.
With regard to the experimental protocol, the implant was initially situated outside of the test setup and its battery was charged for ten minutes. After that, a period of two hours was allowed for the implant to equilibrate to room temperature. Following that, the implant was placed in the calorimetric bath for a half hour prior to the initiation of the experiment. The experiment was operated for a period of 25 minutes during which time the thermocouples were read automatically every 10 seconds. The thermometer was utilized at both the beginning and the end of the experimental period.
At all times during the experiment, the calorimetric-bath thermocouples displayed temperatures uniform to within 0.03˚C. Agreement between the in-bath thermometer and the thermocouples was within the 0.05˚C resolution of the thermometer. The temperatures of the two principle faces of the implant differed by 0.04˚C, reflecting a slight but not unexpected asymmetry. All surfaces of the glass beaker had temperatures about 0.15˚C below that of the calorimetric fluid.
These temperature data were employed to evaluate the thermal energy absorbed at the three measurement zones throughout the duration of the experiment. For this purpose, the first law of thermodynamics was separately invoked for each of energy absorbing media. The appropriate form of the first law is
where Q represents the absorbed energy, m is the mass, c the specific heat, and Tfinal and Tinitial are the respective temperatures at the end and at the start of the experiment. Equation (1) was separately evaluated for the implant proper, the calorimetric fluid, and the glass wall of the beaker. The sum of the three contributions was divided by the duration of the experiment to obtain the rate of heat generation of the implant. For the conditions of the experiment, the heat generation rate was found to be 0.614 W, of which 82%, 10%, and 8% were respectively absorbed by the calorimetric fluid, the glass, and the implant proper.
Although the experimental results for the implant heat generation rate were obtained for an antenna-to-implant separation distance of 20 mm, the results continue to apply without change for other separation distances.
3. The Numerical Model
The targeted results of the numerical simulations were temperature and tissue-damage profiles in the tissue due to heat generated in both the implant and the antenna. Prior to the initiation of the simulations, preparatory work was necessary to establish the respective spatial distributions of the generated heat within the two components (Antenna and Implant shown in Figure 1). The simulation model was also employed for several sensitivity studies to support the accuracy of the results. In particular, systematic variations of both the material properties of the tissue and the rates of heat generation were made to implement the sensitivity studies.
Admittedly, the thicknesses and properties of tissue layers depend on the involved organs and on the age, somatotype, and gender of specific patients. For the sake of generality, a typical set of tissue characteristics was selected for the simulation work  , as listed in Table 1.
Plan views of the antenna and the implant are respectively shown in parts (a) and (b) of Figure 2. In the recommended recharge procedure, the coils of the two components are to be aligned. The images (a) and (b) do not represent the respective sizes of the components. The implant (a) is, in truth, notably smaller than the antenna (b).
In order to facilitate the numerical simulations, it is necessary to create geometric facsimiles of the devices. The facsimiles maintain all important geometric features but omit small and/or thermally insignificant features. Plan views of the facsimile geometries
Figure 1. Illustration of the positioning of the antenna and implant in the surrounding tissue.
Figure 2. Plan views of the (a) implant and (b) antenna (not to scale).
Table 1. Properties and dimensions of tissue layers  .
are shown in Figure 3. There, the external casings have been made semi-transparent so as to reveal the innards of the devices. The components are identified by annotations and accompanying leader lines. For the antenna, the heat generation is confined to a single component (the coil). In contrast, for the implant, the heating is distributed among multiple components.
The physical properties of each of the internal components are required for the simulation. A listing of all the components and their associated properties is provided in Table 2 for the antenna and Table 3 for the implant.
Figure 3. Semi-transparent views of the (a) implant and (b) antenna with identified components.
Table 2. Thermophysical properties of the components of the antenna.
Table 3. Thermophysical properties of the components of the implant.
Heat Generation Protocols
The neuromodulation device in question is designed to operate on a predetermined recharging schedule. The protocol is first to recharge at a low rate for the first minute, followed by a boost (higher) recharge rate for nine minutes and then returning to the initial low rate for the remainder of the recharging process. The heating rates used for the simulations were experimentally determined by electrical measurements and confirmed by the calorimetric experiments discussed earlier.
Table 4 lists the heating rates for the heat-generating components during normal operating conditions. Two entries are listed for each component. The first entry relates to the higher heat dissipation during the boost phase, while the second corresponds to the lower-power of the initial and terminal recharge phases. Values are listed both as absolute heating rates (W). There is only a single heat-generating component in the antenna, whereas five implant components generate heat. The combined heat-genera- tion rate for all components in the implant was 1.102 W during the boost phase and 0.601 W for the low-power regulated phase.
Table 4. Rates of heat generation for electrically active components during normal operation (both high and low heating rates).
As described earlier, careful calorimetric experiments were performed to independently validate the electrical measurements and were within a few %. This outcome provides strong support for the heating rates listed in Table 4.
While the aforementioned discussion relates to normal operating conditions, it is possible to operate the device in a manner that elevates the heat generation rate. For example, higher heating rates can result from misalignment of the coils in the implant and the antenna or from tilting of the antenna away from the plane of the skin surface. The time wise heating pattern for these cases of elevated heating is similar to the normal protocol. A listing of heat generation rates for the elevated heating case is conveyed by Table 5 from which it can be seen that the combined heating rate for the implant components was 1.56 W during the boost phase and 0.604 W during the low-rate phase. Results for both the normal and elevated heating rates will be provided in this report.
4. Mathematical Model
Governing Equations: The model that was employed to determine the time wise evolution of the temperatures within the implant, the antenna, and the surrounding tissue is based on the first law of thermodynamics, the Fourier law of heat conduction, and the Pennes description  of perfusion-based advection of thermal energy in tissue. The corresponding mathematical description of these processes is
This equation applies selectively to each of the three participating region: the tissue, the implant, and the antenna so that there are, in effect, multiple such equations to be solved. These equations are coupled and must be solved simultaneously. The symbols ρ, c, k, and T represent, respectively, the density, specific heat, thermal conductivity, and temperature of the medium in question. The quantity ω is the volumetric rate of blood perfusion per unit volume of tissue, and the subscript b refers to blood. The Sgen term represents the various volumetric heat generation rates that were detailed in the preceding section of the paper.
The term on the left-hand side of Equation (1) represents the thermal inertia of the
Table 5. Rates of heat generation for electrically active components during elevated heating operation operation (both high and low heating rates).
medium, whereas the successive terms on the right are heat conduction, energy transport by blood perfusion, and volumetric heat generation.
Since both the implant and the antenna are inanimate objects, the blood perfusion term in Equation (1) is deleted for those media. For regions without heat generation, Sgen = 0. Living tissue may have a small metabolic heating term which has been ignored in this analysis on the basis of an order of magnitude argument.
The use of the Pennes perfusion model is widely accepted based on its simplicity and record of accurate predictions of biological temperature evolutions. It has been particularly effective at predicting skin temperatures (for example,       -  ).
4.1. Boundary and Initial Conditions
To complete the mathematical description of the problem, thermal boundary conditions were applied to the exposed surfaces of the skin and the antenna as well as at the base of the deep muscle. The selected surface boundary conditions are believed to encompass the expected situations that may occur when recharging is accomplished in practice. For an upper bound on tissue temperatures, an adiabatic condition was enforced at the exposed surfaces which correspond to the user seated on an upholstered surface. The other bounding case was based on heat transfer through a thin layer of fabric and the combined effects of natural convection and thermal radiation to the environment    . For that case, an effective heat transfer coefficient of 6.5 W/m2-˚C was used along with an environment temperature of 20˚C. At the base of the deep muscle, a uniform temperature of 37˚C was applied, while at the outboard lateral surfaces of the solution domain, insulated conditions were imposed.
Attention will now be turned to the initial condition of the temperature prior to the onset of recharging. While the body-core normothermic condition is 37˚C, temperatures significantly lower than this occur near skin surfaces. For this reason, great care was given to the temperature distribution just prior to the onset of heating. To this end, a steady-state temperature distribution was determined by numerical simulation wherein heat generation in both the implant and the antenna were omitted, and the aforementioned thermal boundary conditions were fulfilled. The thus-obtained tissue temperature distribution served as the initial condition for the recharging operation. For all cases, the initial antenna temperature and the ambient temperatures were taken to be 20˚C.
4.2. Issues of Numerical Accuracy
All methods of numerical simulation, including that used in the present study, make algebraic approximations to the governing partial differential equations. In order to be assured that the discretization does not affect the accuracy of the results, a careful study of the spatial and temporal discretization scales was performed. Successive refinements of the mesh dimension scales and the time step size were performed until there was no appreciable change in the resulting temperatures.
For the time-step study, values of 2, 5, 10, 20, and 40 seconds were used. It was found that time-step independence was achieved for step sizes of 20 seconds and less.
With regard to spatial discretization, attention was focused on the maximum tissue- temperature results corresponding to meshes respectively composed of 500,000, 750,000, and 10,500,000 elements. The successive deviations among these maximum values were 0.45˚C and 0.12˚C. All of the final results presented in this paper correspond to a grid having 10,500,000 elements.
In addition to the issue of numerical accuracy, uncertainty with respect to material properties and boundary conditions were evaluated. With respect to properties, simulations were performed with all thermal diffusivity values (all tissues and device components) increased/decreased by 10% through alteration of the thermal conductivity. Assessment of the impacts of property uncertainties is necessary since they vary from person to person and also differ by body region.
4.3. Summary of Investigated Cases
A summary of the investigated cases is provided in Table 6 to aid the later discussion of the results.
4.4. Tissue Damage Evaluation
The temporal and spatial temperature distributions provided by the numerical simulations enable the prediction of thermal tissue damage. There are two broad categories of thermal damage prediction methods. The first is based on an Arrhenius integral model which yields an index of the tissue damage at a selected location as a function of time. That approach, often the termed the Henriques-Moritz method, was developed in a series of studies     . It has found high favor in the literature.
Mathematically, the injury index Ω is calculated by
Here, C(0) and C(t) are the viable cell concentrations at time t = 0 and at a later time t. The symbols A and Ea are experimentally determined tissue-specific constants, R is the universal gas constant, and T is the time- and spatial-dependent tissue temperature
Table 6. Summary of the investigated numerical simulations.
in K. It is widely accepted that tissue necrosis is defined by Ω = 1, which corresponds to 63% of cells injured during a heating process. This criterion will be used in the present study.
An alternative approach, used by regulatory communities to compare the effects of exposures to different temperatures, is the Equivalent Dosimetry Method. Commonly referred to as CEM43, this method relates exposures to an arbitrary temperature to an equivalent exposure at 43˚C    .
Tissue damage evaluation using the CEM43 method is performed with
for isothermal exposures, with T expressed in degrees Celsius and
for exposures during which temperatures may vary continually or in discrete steps. Values of RCEM depend on tissue type, and a large collection of tissue-specific RCEM values can be found in   . Both methods will be employed to estimate the potential of thermal injury during recharge of the Intellis device.
5. The Results
5.1. Spatial Temperatures Variations
The numerical simulations provide spatial temperature distributions at all instances of time throughout the duration of the recharging. Two useful ways to present this information are contour plots whose color tones identify spatial temperature distributions and graphs showing time wise temperature variations at selected locations. The contour diagrams are displayed first.
The temperature contour diagrams were constructed by passing a plane which bisects the antenna, implant, and tissue layers and then plotting the distribution of temperature on the plane. The contour diagrams are conveyed in Figures 4-7, with each figure consisting of (a) and (b) parts. The (a) part corresponds to the moment of termination of the boost phase (10 minutes of recharging), and the (b) part represents the condition at the completion of recharging. A legend connecting color tones and temperatures is displayed in each figure. Tissue temperatures are shown in these figures,
Figure 4. Temperature contours for Case 1 at (a) the cessation of the boost power phase and (b) the termination of the recharge process.
Figure 5. Temperature contours for Case 2 at (a) the cessation of the boost power phase and (b) the termination of the recharge process.
Figure 6. Temperature contours for Case 3 at (a) the cessation of the boost power phase and (b) the termination of the recharge process.
but not implant and antenna temperatures.
Figure 4 conveys results for Case 1 of Table 6, which corresponds to normal heating and properties and to an adiabatic condition at the surfaces exposed to the environment. The temperature levels seen in Part (b) are substantially higher than those of Part (a), reflecting the longer durations of heating for the former. At all heating durations, the temperatures in the tissue that borders the implant are the highest in the entire tissue bed. At the termination of recharging, the maximum tissue temperature is less than 42˚C.
Figure 7. Temperature contours for Case 4 at (a) the cessation of the boost power phase and (b) the termination of the recharge process.
The impact of external heat loss can be seen by comparing Figure 4 (Case 1, adiabatic) and Figure 5 (Case 2, external convection and radiation). Up to the end of the boost phase (10 minutes of heating), there is very little difference in the tissue temperatures corresponding to adiabatic and non-adiabatic conditions (Parts (a)). However, at the termination of recharging, the external heat transfer has brought about significant lowering of the tissue temperatures (Parts (b)).
Figure 6 and Figure 7 both correspond to elevated heating conditions (Table 6), respectively for adiabatic (Case 3) and non-adiabatic (Case 4) external conditions. In general, the elevated heating rates do give rise to somewhat higher tissue temperatures. With regard to the effects of external thermal conditions, they are the same as those identified for Figure 4 and Figure 5.
5.2. Temporal Temperature Variations
A major safety concern with the use of externally recharged implants is the possibility of attaining excessively high tissue temperatures. The temperature results already presented clearly demonstrate that tissue temperatures will not exceed 42˚C. This temperature is widely regarded as benign. Support of this assessment will shortly be conveyed when tissue damage metrics will be presented.
The importance of the maximum tissue temperature has motivated Figure 8, where the time wise variation of the maximum tissue temperature is displayed for Cases 1, 2, 4, and 5 of Table 6. Comparisons of the results for Cases 1 and 2 and Cases 3 and 4 illustrate the effect of adiabatic vs. non-adiabatic external boundary conditions. The difference in boundary condition does not assert itself during the boost heating phase, but subsequently brings about substantial deviations between the respective maximum temperatures. Elevated vs. normal heating is highlighted by comparing the curves for Cases 1 and 3 and Cases 2 and 4. Clearly, elevated heating has an important impact on the maximum temperature for all times during the recharging period.
Perhaps the most important message of Figure 8 is that the highest value of the temperature maxima does not exceed 42˚C, thereby corroborating the observations noted in Figures 4-7.
An evaluation of properties variations is presented in Figure 9, which includes results
Figure 8. Time wise variations of the maximum tissue temperature for various operating conditions.
Figure 9. Effect of material property variations on the maximum tissue temperature.
for all three property-related cases (Case 2 of Table 6) and a refined temperature scale. It is seen from the figure that the largest differences in the temperatures for the three cases are approximately 0.7˚C.
In order to show a more complete view of the time wise temperature evolution, it is relevant to provide information at important locations within the solution domain. The selected locations, displayed in Figure 10, are deployed along a line extending from the antenna to the implant. Point 1 lies in the plane of the outward facing surface of the implant, point 3 is at the interface between the skin and subcutaneous fat layer, and point 2 is equidistant between locations 1 and 3.
The time wise temperature variations at locations 1, 2, and 3 are displayed in Figure 11 for the Case 2 operating conditions of Table 6. For perspective, a curve is shown that depicts the maximum temperatures for Case 2. It is seen from the figure that among the three points in the space between the implant and the antenna, the highest temperatures occur at the surface of the implant. These temperatures are lower than the maximum temperatures for the given operating condition, which occurs just beneath the implant. There is a substantial drop off in temperature in the direction toward the skin.
5.3. Thermal Injury Evaluation
The spatial and temporal temperature results provided by the numerical simulations
Figure 10. Identification of three in-tissue locations between the implant and the antenna chosen for display of time wise temperature variations.
Figure 11. Time wise temperature variations for Case 2 operating conditions at in-tissue locations identified in Figure 10. The uppermost curve corresponds to the maximum temperatures for Case 2 operating conditions.
enable the evaluation of the potential for thermal injury. As discussed earlier, two independent approaches are used here for the evaluation. To facilitate the Arrhenius model, values of A and Ea are needed. For a conservative approach (overprediction of thermal injury), the values selected are A = 1.30 × 1095 s−1, and Ea = 6.04 × 105 J/mole  .
The execution of the Arrhenius integral, Equation (2), for the highest temperature operating conditions which correspond to Case 3, the injury metric Ω was found to be approximately 0.02. This value is well below the tissue necrosis value of 1.0.
Reinforcement of this outcome is obtained by independent calculation of the CEM43 from Equation (4) with RCEM = 0.233. This calculation yields a CEM43 value of approximately four minutes, which is also well below the potential level for thermal injury.
From these analyses, regardless of the tissue injury prediction methodology and utilizing conservative estimates for the highest temperature operating conditions considered here, it can be concluded regardless of the tissue injury prediction method used, no evidence for a potential for thermal injury was identified when the Intellis device undergoes recharge.
6. Concluding Remarks
Neuromodulation implants generally contain relatively short-lived batteries that must be frequently recharged. The recharging is accomplished by magnetic means whereby a coil housed in the antenna produces lines of force which cut across an implant coil. Both coils generate heat, providing the possibility of overheating tissue in the neighborhood of these devices. Present neuromodulation systems provide off-on recharging dependent upon sensed values of the tissue temperature. A new conceptual approach embodied in the Medtronic Intellis is to combine the sensed tissue measurement approach with a preprogramed recharging protocol. The goal of the research reported here is to determine whether the preprogrammed time-varying recharge protocol can give rise to tissue temperatures that are either absolutely safe or may otherwise cause tissue damage.
In sharp contrast to previous studies of dated recharging protocols, the present approach is based on a highly accurate model of all the thermal processes that are activated by the recharging. The earlier approaches bypassed these processes and depended solely on a global experiment. The present detailed model was implemented by numerical simulations performed for several realistic operating conditions and provided in-depth spatial and temporal tissue temperature distributions. Those temperature distributions were employed to estimate possible tissue damage by making use of two independent methodologies. Independent calorimeter-based experiments were performed to provide validation for the calculated rates of heat generation in the coils of the implant.
None of the aforementioned spatial and temporal tissue temperature distributions gave rise to temperatures above 42˚C, regardless of the operating conditions. It is broadly accepted that tissue temperatures less than 43˚C are not injurious to tissue. Numerical values of thermal tissue damage metrics were determined and compared with accepted values which correspond to the absence and the presence of tissue damage. Both the tissue temperature results and the thermal damage metrics strongly support the total absence of tissue injury when time-varying preprogrammed protocols are used in the recharging of neuromodulation implant-encased batteries. The methodology carried out here has been validated in multiple settings including animal model experiments  . That work, and other work that is currently in submission have used lamb and pig surrogates in experimentation. They have found that measured temperatures and temperatures calculated using the numerical method described in this study are in excellent agreement with maxima values within approximately 0.1C.
Support for this research project was provided by Medtronic.
 Lovik, R., Abraham, J. and Sparrow, E. (2009) Potential Tissue Damage from Transcutaneous Recharge of Neuromodulation Implants. International Journal of Heat and Mass Transfer, 52, 3518-3524.
 Smith, D., Lovik, R., Sparrow, E. and Abraham, J. (2010) Human Tissue Temperatures Achieved during Recharging of New-Generation Neuromodulation Devices. International Journal of Heat and Mass Transfer, 53, 3292-3299.
 Lovik, R., Abraham, J. and Sparrow, E. (2011) Surrogate Human Tissue Temperatures Results from Misalignment of Antenna and Implant during Recharging of a Neuromodulation Device. Journal of Neuromodulation, 14, 501-511.
 Chen, S., Li, Q., Wang, W., Ma, B., Hau, H. and Li, L. (2013) In Vivo Experimental Study of Thermal Problems for Rechargeable Neurostimulators. Journal of Neuromodulation, 16, 436-441.
 Lovik, R.D., Sparrow, E.M., Abraham, J.P., Zelmer, C., Oh, S., Friend, K., et al. (2011) Effect of Component Misalignment on Human Tissue Temperatures Associated with Recharged Neuromodulation Devices. Journal of Medical Devices, 5, Paper No. 0207516.
 Abraham, J.P., Plourde, B.D., Vallez, L.J., Stark, J.R. and Diller, K.R. (2015) Estimating the Time and Temperature Relationship for Causation of Deep-Partial Thickness Skin Burns. Burns, 41, 1741-1747.
 Abraham, J.P., Plourde, B.D., Vallez, L.J., Nelson-Cheeseman, B.B., Stark, J.R., Gorman, J.M. and Sparrow, E.M. (In Press) Skin Burn. In: Shrivastava, D., Ed., Theory and Application of Heat Transfer in Cells and Organs, Wiley.
 Johnson, N.N., Abraham, J.P., Helgeson, Z.I., Minkowycz, W.J. and Sparrow, E.M. (2011) An Archive of Skin-Layer Thicknesses and Properties and Calculations of Scald Burns with Comparisons to Experimental Observations. Journal of Thermal Science and Engineering Applications, 3, Paper No. 011003.
 Dewhirst, M.W., Abraham, J.P. and Viglianti, B.L. (2015) Evolution of Thermal Dosimetry for Applications of Hyperthermia Treatment to Cancer. Advances in Heat Transfer, 47, 397-421.
 Abraham, J.P. and Sparrow, E.M. (2007) A Thermal Ablation Model including Liquid-to-Vapor Phase Change, Necrosis-Dependent Perfusion, and Moisture-Dependent Properties. International Journal of Heat and Mass Transfer, 50, 2537-2544.
 Diller, K., Hayes, L. and Blake, G. (1991) Analysis of Alternate Models for Simulation Thermal Burns. Journal of Burn Care & Rehabilitation, 12, 177-189.
 Reinders, D., Baldwin, S. and Bert, J. (2003) Endometrial Thermal Balloon Ablation Using a High Temperature, Pulsed System: A Mathematical Model. Journal of Biomechanical Engineering, 125, 841-851.
 Ng, E. and Chua, L. (2002) Prediction of Skin Burn Injury. Part 2: Parametric and Sensitivity Analysis. Proceedings of the Institution of Mechanical Engineers, 216, 171-183.
 He, X., McGee, S., Coads, J., Schmidlin, F., Iaizzo, P., Swanlund, D., Klugey, S., Rudie, E. and Bischoff, J. (2004) Investigation of the Thermal and Tissue Injury Behavior in Microwave Thermal Therapy Using a Porrcine Kidney Model. International Journal of Hyperthermia, 20, 567-593.
 Sun, J., Zhang, A. and Xu, L. (2008) Evaluation of Alternate Cooling and Heating for Tumor Treatment. International Journal of Heat and Mass Transfer, 51, 5478-5485.
 Dai, W., Wang, H., Jordan, P., Mickens, R. and Bejan, A. (2008) A Mathematical Model for Skin Burn Injury Induced by Radiation Heating. International Journal of Heat and Mass Transfer, 51, 5497-5510.
 Khanafer, K. and Vafai, K. (2009) Synthesis of Mathematical Models Representing Bioheat Transfer. In: Minkowycz, W., Sparrow, E. and Abraham, J., Eds., Advances in Numerical Heat Transfer, Vol. 3, CRC Press, London.
 Henriques, F.C. (1947) Studies of Thermal Injury V. The Predictability and Significance of Thermally Induced Rate Processes Leading to Irreversible Epidermal Injury. Archives of Pathology, 43, 489-502.
 Henriques, F.C. and Moritz, A.R. (1947) Studies of Thermal Injury: I. The Conduction of Heat to and Through Skin and the Temperatures Attained Therein. A Theoretical and an Experimental Investigation. American Journal of Pathology, 23, 530-549.
 Moritz, A.R. and Henriques, F.C. (1947) Studies of Thermal Injury II. The Relative Importance of Time and Surface Temperature in the Causation of Cutaneous Burns. American Journal of Pathology, 23, 695-720.
 Sapareto, S. and Dewey, W. (1984) Thermal Dose Determination in Cancer Therapy. International Journal of Radiation Oncology, Biology, Physics, 10, 787-800.
 Dewhirst, M., Viglianti, B., Lora-Michiels, M., Hanson, M. and Hoopes, P. (2003) Basic Principles of Thermal Dosimetry and Thermal Thresholds for Tissue Damage from Hyperthermia. International Journal of Hyperthermia, 19, 267-294.
 Yarmolenko, P., Moon, E., Landon, C., Manzoor, A., Hochman, D., Viglianti, B. and Dewhirst, M. (2011) Thresholds for Thermal Damage to Normal Tissues: An Update. International Journal of Hyperthermia, 27, 320-343.
 Abraham, J.P. and Plourde, B.D. (2016) Validation of Numerically Simulated Tissue Temperatures during Transcutaneous Recharge of Neurostimulation Systems. Journal of Nueromodulation, 19, 161-170.