Determining the thermal characteristics of breast cancer based on high-resolution infrared imaging, 3D breast scans, and magnetic resonance imaging

For over the three decades, various researchers have aimed to construct a thermal (or bioheat) model of breast cancer, but these models have mostly lacked clinical data. The present study developed a computational thermal model of breast cancer based on high-resolution infrared (IR) images, real three-dimensional (3D) breast surface geometries, and internal tumor definition of a female subject histologically diagnosed with breast cancer. A state-of-the-art IR camera recorded IR images of the subject’s breasts, a 3D scanner recorded surface geometries, and standard diagnostic imaging procedures provided tumor sizes and spatial locations within the breast. The study estimated the thermal characteristics of the subject’s triple negative breast cancer by calibrating the model to the subject’s clinical data. Constrained by empirical blood perfusion rates, metabolic heat generation rates reached as high as 2.0E04 W/m3 for normal breast tissue and ranged between 1.0E05–1.2E06 W/m3 for cancerous breast tissue. Results were specific to the subject’s unique breast cancer molecular subtype, stage, and lesion size and may be applicable to similar aggressive cases. Prior modeling efforts are briefly surveyed, clinical data collected are presented, and finally thermal modeling results are presented and discussed.

In the context of breast cancer screening and diagnosis, infrared (IR) imaging, also referred to as breast thermography or digital infrared thermal imaging (DITI), is an imaging technique where IR images are taken of a patient's breasts, referred to as "thermograms. " When used adjunctively in a screening or diagnostic environment, IR images are assessed for thermal abnormalities potentially indicating breast cancer. The method is based on the thermographic detection of cancer at the skin resulting from pathophysiologic changes within the breast caused by cancer; these changes include the metabolic and vascular changes associated with cancer development 1 . Lozano and Hassanipour provided a detailed and objective survey of prior clinical studies, both favorable and unfavorable, regarding IR breast thermography 2 . Owing to mixed clinical findings over the past six decades, the United States Food and Drug Administration (FDA) approved IR thermography in 1982 only as an adjunct to other screening modalities (e.g., mammography and ultrasound) 3 . Today, the FDA warns that IR thermography is not a substitute for mammography and should not be used as an isolated diagnostic procedure 4,5 . The primary motivation for using IR thermography has been the fact that IR imaging does not expose the patient to any amount of ionizing radiation, contrary to mammography.
The present study addresses a two-fold, coupled problem: The lack of knowledge of the thermal characteristics of breast cancer and the lack of real clinical data use in the computational thermal modeling of breast cancer, each of which are discussed in detail. First, data regarding the thermal characteristics of cancer, like the metabolic heat generation rate (a critical thermal modeling parameter), are scant and limited in the literature. In fact, the www.nature.com/scientificreports www.nature.com/scientificreports/ metabolic heat generation rate of living tissue (normal or cancerous) has been the most elusive term to quantify in any bioheat formulation 6 . In 1980, Gautherie provided the most extensive (and only) data available to date regarding the metabolic heat generation rate of cancerous breast tissue 7 . Gautherie reported an empirical relationship between the metabolic heat generation rate of a breast carcinoma and its volume-doubling time. Gautherie observed a hyperbolic relationship between these two parameters based on data from 84 breast cancer patients with cancerous growth sizes ranging between 0.9-3.8 cm and doubling times of 49-676 days. Gautherie offered an approximate range for the metabolic heat generation rate of breast cancer: 4,800-65,000 W/m 3 . No other data have been reported in the literature for this parameter. Absent of any other data, prior breast cancer modeling efforts have historically relied on Gautherie's relation when modeling the breast with cancer.
However, of these approximately 30 efforts, only a few prior models have been constructed using experimental data (here, clinical IR images). Ng and Sudharsan matched their steady-state numerical model with real IR images from 3 female subjects; a mannequin of brassiere size 34 C was used for the breast geometry, and the tumor size and location was approximated [10][11][12][13]45 25,46 .
While these efforts in the past have aimed to develop a thermal model of the breast with cancer, the majority have been academic exercises in simulation. These modeling efforts have collectively had the following limitations: 1. The lack of clinical data with which to calibrate the model (here, the thermography of the breast with cancer). 2. The lack of real breast shapes (i.e., surface geometries), which are unique and vary per individual. 3. The lack of real tumor definition (i.e., size and (x, y, z) spatial location within the breast). 4. The use of Pennes' bioheat equation as the governing equation, which inherently has limitations and simplifications. 5. The lack of real internal vasculature and blood flow field (i.e., hemodynamic) information in the breast.
Therefore, the goal of the present study was to determine the thermal characteristics of breast cancer by constructing a computational thermal (or bioheat) model calibrated to real clinical data. Specifically, the study aimed to quantify the thermal characteristics-i.e., blood perfusion rate, ω  b , and metabolic heat generation rate, ′′′ Q m -of normal and cancerous breast tissues that best matched surface temperatures observed on IR images. Clinical data consisted of high-resolution IR images, three-dimensional (3D) breast surface geometries, and internal tumor definition from a female subject histologically diagnosed with breast cancer. A state-of-the-art IR camera recorded IR images, a 3D scanner recorded breast surface geometries, and magnetic resonance (MR) imaging data provided tumor definition such as size and (x,y,z) spatial location. 3D scans and tumor definition served as geometric inputs to the model, whereas IR images were used to calibrate the model. This study is distinguished from prior thermal modeling efforts by its use of clinical data to construct and calibrate the model. The present study aimed to address the first three limitations noted above.

Results
Solutions of ( ω ′′′Q , m b ) pairs were found for a range of arterial temperatures, T A , and convective heat transfer coefficients, h (or HTC). Results are presented in the following: Fig. 1 plots solutions for the right (Fig. 1a) and left (Fig. 1b) breast thermal models, which represent thermal characteristics of normal and cancerous breast tissues, respectively; each ( ω ′′′Q , m b ) solution pair is plotted along its corresponding T A line. For the left breast model, solutions are presented for each of the three tissues in the model (lower cancerous, upper cancerous, and non-cancerous tissues). Figure 2 demonstrates surface temperature contours for a representative solution from the right (Fig. 2a) and left (Fig. 2b) breast models, with the IR image shown for reference representing the true temperature distribution. Finally, Table 1 outlines all solutions found for each breast model.
Regarding the right breast model, the accuracy of each solution was verified based on the calculated rate of convective heat loss from the anterior surface (based on the 33.8°C mean surface temperature over the entire breast); solutions were accurate to within 0.003 W of the theoretical value, which varied by HTC. The maximum internal temperature for all cases was 36.0°C. No solutions were found for HTC values greater than 8.0 W/m 2 -K given the ( ω ′′′Q , m b ) ranges selected for either of the two following reasons: First, the maximum internal breast temperature constraint was not satisfied (see Eq. 6), or second, insufficient internal heat was generated to overcome the convective surface cooling caused by the high HTC value. Similarly, no solutions were found for www.nature.com/scientificreports www.nature.com/scientificreports/ T A = 36.0°C that satisfied the maximum internal temperature constraint, which was attributed to the posterior wall boundary condition (see Eq. 3).
Regarding the left breast model, solutions were found for all HTC and T A cases considered, contrary to the right breast model for the reason that no constraints were enforced in this model. Noteworthy were the large magnitudes of ′′′ Q m in lower cancerous breast tissue relative to other tissues required to achieve sufficiently high surface temperatures at the nipple, particularly for lower T A values. This was attributed to the local region of skin thickening at the nipple. Maximum internal breast temperatures for each case varied between 42.2-50.8°C. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
The solutions of normal and cancerous breast tissues presented herein represent a general range of possibilities regarding their thermal characteristics, depending on arterial temperature and environmental conditions (i.e., convective heat transfer coefficient and ambient temperature). For cancerous breast tissue specifically, it is important to note that solutions were subject-specific and cannot necessarily be applied to all breast cancers. However, the subject modeled herein represented a triple negative breast cancer and the thermal characteristics reported may be applicable to this molecular subtype of breast cancer in other subjects. In contrast, for normal breast tissue, solutions found for metabolic heat generation rate were generally less than 20,000 W/m 3 , which were comparable to the maximum value of 12,000 W/m 3 reported by Gautherie 7 .
Mathematically, the Pennes temperature-dependent source term, , can represent either a heat source or heat sink, depending on whether the local tissue temperature is less than T A (heat source) or greater than T A (heat sink). This perfusion term balances the metabolic heat generation term, ′′′ Q m , in order to bring the local tissue temperature closer to T A according to the magnitude of ω  b . Physiologically, this term represents the thermoregulatory mechanism in perfused tissue, and without this term, unrealistically high tissue temperatures may be encountered. Accordingly, for the right breast model, any metabolic heat generation at T A ≥ 36.0°C violated the maximum internal breast temperature constraint (see Eq. 6) and therefore yielded no solution, owing to the fact that a 36.0°C posterior wall temperature boundary condition was assigned based on measurements reported by Gautherie (see Eq. 3) 7 . For the left breast model, the two hotspots caused by the internal cancerous mass matched the IR image. The medial hotspot, which was due to a local blood vessel and not any underlying mass, was not recreated. Lighting shown to demonstrate 3D surfaces.
www.nature.com/scientificreports www.nature.com/scientificreports/ The slope of each T A series in Fig. 1 was unique to the ( ω ′′′Q , m b ) endpoints selected for that individual series (see Table 1). The ranges of ( ω ′′′Q , m b ) were determined based on preliminary trial-and-error modeling that yielded a solution under the most extreme condition modeled (here, h = 8.0 W/m 2 -K, T A = 33.0°C). Various other ( ω ′′′Q , m b ) ranges per T A series were found to yield fewer solutions than the ranges assumed (not reported). Therefore, it was conceivable that an optimal ( ω ′′′Q , m b ) range existed that could yield even more solutions for greater HTC values than those reported herein. However, ascertaining an optimal ( ω ′′′Q , m b ) range would require a numerical parametric sensitivity study, which was beyond the scope of the present modeling study.
In the right (normal) breast thermal model, the following inter-variable relationships were observed: As T A decreased, not enough heat could be generated internally to produce a sufficiently high mean surface temperature of 33.8°C. That is, fewer solutions were found as T A decreased. This was attributed to the fact that the Pennes perfusion term (here, a heat sink), which is a function of T A , counteracted any contribution from the ′′′ Q m source term. The same observation was made as HTC increased for the same mathematical reasoning. Further, increasing ′′′ Q m at a constant ω  b yielded higher surface temperatures, which varied by case. Conversely, increasing ω  b at a constant ′′′ Q m brought tissue temperatures closer to T A according to the magnitude of ω  b . Interestingly, in the left (cancerous) breast thermal model, the elevated perfusion rates characteristic of cancerous tissue were required for all tissues-including non-cancerous breast tissue-in order to find solutions. That is, normal perfusion rates did not suffice in raising the mean surface temperature (excluding hotspots) to the 34.9°C metric. This observation was attributed to the fact that normal perfusion rates, which are approximately ten times lower than cancerous tissue, limited the magnitude of the Pennes perfusion term (here, a heat source). This finding indicated that for invasive and aggressive breast cancers (as observed in the subject), the entire breast may exhibit elevated blood perfusion rates characteristic of cancerous tissue, not only the malignant tumor itself. This finding was reasonable given the 1.2°C mean temperature difference between breasts and generalized hyperthermia on the subject's left breast based on the subject's IR images (see Fig. 3a,b and Table 2).
It was also observed in the left breast thermal model that given the physiologic modeling assumptions made, the ′′′ Q m ranges of cancerous tissue reported by Gautherie did not produce sufficiently high surface temperatures

Ranges Solutions
Breast model www.nature.com/scientificreports www.nature.com/scientificreports/ to meet the hotspot target metrics for the subject modeled in this study. In fact, based on the 10 cm tumor size exhibited by the subject (with an unknown doubling time), the appropriate metabolic heat generation rate per Gautherie's relation was unclear 7 . Also, it is impractical to clinically ascertain tumor doubling times with routine . Infrared images (a,b) and 3D breast surface scans (c) of Subject 03. IR images in (a) were exported directly from FLIR ResearchIR Max; IR images in (b) were processed using MATLAB to modify background color to black. Clear thermovascular asymmetry was observed between breasts due to generalized hyperthermia in subject's left breast (see Table 2 for summary of temperatures). Physical deformity in subject's left breast due to malignancy was evident at presentation during research procedures. High-resolution IR images available for download online (see Supplementary Figs. S1-S6).  46 .
In this study, the subject's medial hotspot was not modeled for the reason that it was attributed to a local blood vessel and not to the underlying malignant tumor. Accurately recreating the subject's medial hotspot would require incorporating the internal vasculature (i.e., blood vessel branching structure) into the thermal model. Therefore, these results highlight the inherent limitation to implementing Pennes' bioheat equation. Modeling the computational fluid dynamics (CFD) of the detailed 3D vasculature reconstructed from MR imaging-with accurate hemodynamic velocities and arterial temperatures-has recently made progress and is needed to achieve a higher degree of fidelity in the thermal modeling of any tissue 47,48 .
The following were three limitations to the modeling study presented herein. First, as previously mentioned, modeling results were subject-specific and are not necessarily applicable to all breast cancers. Second, the metabolic heat generation rate of skin was unknown. No empirical data were found in the literature, and in an effort to avoid modeling arbitrary values, it was consequently assumed to be zero. Skin tissue likely has a degree of metabolic heat generation, even if it were to have a minimal impact on surface temperatures. Accordingly, the thermal characteristics of breast cancer reported herein may be interpreted as upper-limit values, pending empirical data on skin's metabolic heat generation rate. Third, high maximum internal temperatures were observed for lower cancerous tissue in the left breast model, despite the reasonable matching of surface temperatures on the breast. In fact, internal temperatures in cancerous tissue exceeding 40°C approach hyperthermia treatment ranges [49][50][51][52] . It is expected that incorporating the detailed internal vasculature would resolve this limitation and eliminate the need for Pennes' bioheat equation altogether.
In summary, the following are key takeaways from the thermal modeling results presented herein: Potential clinical applications for the computational modeling of breast cancer include: Monitoring tumor response to cancer treatment over other imaging modalities; tracking tumor growth over time, including from presentation to surgical removal 53 ; and hyperthermia cancer treatment [49][50][51][52] .
In conclusion, the present study developed a computational thermal model of breast cancer based on high-resolution IR images, real 3D breast geometries, and real tumor definition from a breast cancer subject. The novelty of the study was its use of clinical data to construct and calibrate the model. The goal of the study was to quantify the thermal characteristics of breast cancer. The results reported herein represent a first step toward the accurate computational thermal modeling of breast cancer. Clearly, the female breast is a complex, inhomogeneous tissue. Thus, results provide estimates of the ranges expected regarding the metabolic heat generation rate of breast cancer, particularly for the aggressive triple negative subtype exhibited by the subject. Future directions of this study include validating the calibrated model with future female subjects presenting with similar breast cancer characteristics (here, molecular subtype, stage, and tumor size) and incorporating the real 3D breast vasculature into the model as discussed earlier 47,48 .

Methods
The present study was comprised of a clinical study and computational modeling effort. Clinical data were used to develop the thermal model. Specifically, 3D surface scans defined the breast geometry; diagnostic imaging procedures defined the tumor geometry; and IR images were used to calibrate the model. Procedures and equipment. Subjects entered and sat in the examination room and disrobed above the waist. Subjects had previously been thermally acclimated to the environment. High-resolution IR images were recorded www.nature.com/scientificreports www.nature.com/scientificreports/ with a FLIR A655sc infrared camera (FLIR Systems Inc., Wilsonville, OR, USA). The IR camera had a 640 × 480 resolution and 30 mK (0.030°C) sensitivity. The standard lens on the IR camera was used (25° field of view, 24.6 mm focal length) and subjects were positioned approximately 1 m away from the IR camera. IR images were processed using FLIR ResearchIR Max 4.40.8.28, a vendor-provided proprietary software. The emissivity was set to a constant value of 0.98 as has been previously determined empirically for human skin [54][55][56] . Additionally, IR images were post-processed using MATLAB 2019b (The MathWorks Inc., Natick, MA, USA).
Then, 3D breast surface scans were recorded with an Artec Eva Lite (Artec Europe SARL, Luxembourg), a portable, handheld 3D scanner with an accuracy of 100 μm and resolution of 500 μm. The 3D scanner was held approximately 50-60 cm away from subjects' breasts and gently swept in an airbrush motion for approximately 30-60 seconds to fully record the 3D topography of both breasts. 3D scans were processed using Artec Studio 13 Professional, a vendor-provided proprietary software.
Finally, the ambient temperature in the clinical examination room ranged between 22.0-24.0°C. Ambient temperature was measured with an Omega HH506RA digital thermometer (Omega Engineering Inc., Norwalk, CT). The thermometer had a 0.1°C resolution and ±0.3°C accuracy. A T-type thermocouple with 36 AWG wire (0.13 mm diameter) with an exposed junction was used with the digital thermometer (Omega 5SRTC-TT-T- . The thermocouple had a ±0.5°C accuracy ("special limits of error" grade wire). Manufacturer datasheets for all equipment are publicly available online.
Clinical data for representative subject. Clinical data and modeling results are herein presented for a representative subject, Subject 03. Subject 03 was selected from among the study population due to the subject's clear thermal asymmetry between breasts, visible physical deformity, and aggressive breast cancer subtype. Subject 03 (herein referred to as "the subject") was histologically diagnosed with Stage 4 breast cancer in the left breast (invasive ductal carcinoma of no special type, Grade 2, triple negative receptor status).
The subject exhibited a mean thermal asymmetry between breasts of 1.2°C. In addition to the generalized hyperthermia, the subject's left breast exhibited three localized (or focal) regions of increased temperatures: (1) the nipple region (herein, "nipple hotspot"); (2) approximately 8.0 cm superior to the nipple at 11 o' clock (herein, "superior hotspot"); and (3) approximately 5.0 cm medial to the nipple at 10 o' clock (herein, "medial hotspot"). The nipple and superior hotspots were caused by the underlying malignant lesion (i.e., tumor or mass); that is, the cancer was located directly behind these two regions. The medial hotspot, in contrast, was not due to any underlying mass, but rather was attributed to a local blood vessel. The subject's IR images are presented in Fig. 3a,b and are summarized in Table 2. The subject's two hotspots that were due to the underlying cancerous mass were the focus of the modeling effort.
The subject's left breast also exhibited an overall skin thickening with an average thickness of approximately 5.0 mm based on MR images (in contrast to the right breast with an average skin thickness of 1.  www.nature.com/scientificreports www.nature.com/scientificreports/ Meshing and ANSYS Fluent. ANSYS Meshing was used to generate the mesh for the geometry and ANSYS Fluent was used to define the governing equation, assign thermal properties, apply boundary conditions, and numerically solve the governing equation. ANSYS Fluent implements the finite volume method (FVM) to integrate the governing conservation laws over each control volume cell in the geometric domain and numerically solve the resulting discretized equations.
Governing equation and boundary conditions. Absent of the subject's real internal vasculature definition, Pennes' bioheat equation was solved as the governing equation in order to obtain the surface temperature distribution of the breast (see Eq. 1). Here, T represents the temperature of the tissue, ρ represents its density, c represents its specific heat, k represents its thermal conductivity, ω  b represents its blood perfusion rate, T A represents the temperature of its arterial blood supply, and ′′′ Q m represents its metabolic heat generation rate. Pennes' bioheat equation is a modified form of the heat equation that includes an additional temperature-dependent source term, The justification for using Pennes' bioheat equation as an analysis tool, despite its physiologic assumptions and limitations, has been previously established [57][58][59][60] .
In order to solve this governing equation using ANSYS Fluent, a user-defined function (UDF) that enforced the Pennes temperature-dependent perfusion source term in the governing equation was required. The UDF is a short code in the C++ programming language using the appropriate Fluent macros; in this case, the "DEFINE_ SOURCE" macro was used. The UDF was written and loaded into the model.
The numerical validity of the custom UDF was confirmed by solving a simple, ad hoc hemispherical model using both ANSYS Fluent 2019 R2 and COMSOL Multiphysics 5.4 and comparing results. COMSOL Multiphysics contains a built-in bioheat transfer module that applies the temperature-dependent source term to the numerical solution. By applying the same mesh parameters, material properties, boundary conditions, thermal characteristics, and numerical solution settings in both software, the steady-state models were observed to be identical to within 0.01°C with regard to the maximum and minimum temperature values and internal and surface temperature distributions; results were omitted for brevity. This exercise validated the accuracy of the UDF written for ANSYS Fluent.
Because the steady-state model did not simulate any fluid flow, the ANSYS Fluent flow solver was turned off. The energy convergence criterion was set to 1E-12, which was sufficient for a numerically accurate solution. The energy under-relaxation factor, which controls the change in temperature over each iteration, was set to 0.95 in order to avoid divergence. Solutions generally required 20 iterations to satisfy this convergence criterion using the parallel processor solver setting on a quad-core processor computer. Results were verified to vary by less than 0.001°C with a convergence criterion of 1E-15.
Finally, the boundary conditions selected were consistent with prior modeling studies. First, a variable convective heat transfer coefficient was assigned to the anterior surface area of the breast (n normal direction at any point on the surface; see Eq. 2). A reference ambient temperature ( ∞ T ) of 23.0°C was assigned, based on the room temperature measured during the IR imaging procedure. Second, the temperature of the surface posterior to breast tissue ("posterior wall") was held at a constant and uniform temperature of T P = 36.0°C (i.e, isothermal), based on in vivo temperature measurements of deep breast tissue reported by Gautherie (see Eq. 3) 7 . This posterior wall represented the anterior surface of the pectoral muscle on the thoracic wall and was constructed based on the subject's MR imaging. Additionally, the connecting surfaces between the anterior and posterior surfaces were assumed adiabatic (see Eq. 4). Geometry. Real 3D breast surface geometries as recorded by the 3D scanner were imported into the thermal model (see Fig. 3c). This allowed for an accurate representation of the subject's unique breast shape (i.e., the natural ptosis of the breast) and convective surface area. With regard to modeling workflow, Artec Studio 13 Professional exported triangular mesh files (STL format) from raw 3D scan data. STL mesh files were then imported into Solidworks 2018 (Dassault Systemes, Velizy-Villacoublay Cedex, France), a computer-aided design (CAD) software, to generate a smooth, continuous surface using the "ScanTo3D" add-in. Then, using standard functions in Solidworks, a solid body was generated from the surface. Finally, the solid body CAD model was imported into ANSYS Meshing for discretization and later into ANSYS Fluent for numerical solution.
A separate CAD model was constructed for each of the subject's breasts (see Fig. 5a-c). Each breast CAD model consisted of a skin layer on the anterior surface, breast tissue, and cancerous tissue (if applicable) www.nature.com/scientificreports www.nature.com/scientificreports/ representing the tumor. The posterior surface boundary in each model was created based on the natural curvature of the pectoral muscle's anterior surface according to corresponding MR images.
Skin layers were modeled based on MR images. The subject's right breast was modeled with a uniform skin thickness of 1.5 mm. The subject's left breast, however, was modeled with a uniform skin thickness of 5.0 mm that locally thickened to 10.0 mm (nipple region) and 6.5 mm (superior to nipple) as measured on MR images. Similarly, the location and irregular geometry of the internal tumor were approximated based on multiple MR image cross-sections (see Fig. 5b).
Radii from the posterior surface to the anterior surface of each CAD model were approximately 5.5 cm (left breast) and 5.3 cm (right breast). Bounding box dimensions (i.e., geometric envelope) of each model, in transverse × axial × depth format, were approximately 16.3 × 15.4 × 9.5 cm (left breast) and 16.0 × 14.1 × 8.7 cm (right breast). Solid volumes for each breast model were approximately 876.1 cm 3 (left breast) and 588.2 cm 3 (right breast); the discrepancy in volumes was explained by the deformity in the left breast caused by malignancy (see Fig. 3c). Exterior convective surface areas were approximately 312.0 cm 2 (left breast) and 267.5 cm 2 (right breast), derived from 3D surface scans.
Due to the irregular geometry of each breast, 1.0 mm tetrahedral mesh elements (maximum element size) were applied in ANSYS Meshing, which were later converted into polyhedral elements in ANSYS Fluent for faster computational processing. Results between tetrahedral and polyhedral mesh types were verified to be consistent to less than 0.005°C, which was approximately 10 times less than the IR camera's sensitivity. Final polyhedral mesh sizes ranged between 5.1-7.7 million nodes, 6.1-9.1 million faces, and 0.9-1.3 million cells, depending on breast model (see Fig. 5d,e). Solutions reported herein were verified to be independent of mesh size; temperatures varied by less than 0.004°C when mesh size was reduced to 0.7 mm and 0.5 mm tetrahedrons.
Thermophysical properties. Thermophysical properties reported in the literature based on empirical data were assumed in the model as outlined in Table 3. Each tissue was modeled with homogeneous properties, with the exception of cancerous tissue as discussed in the proceeding section. Due to the subject's heterogeneously dense classification of breast tissue based on mammograms (not shown), glandular breast tissue properties were assigned. Blood perfusion rates of cancerous breast tissue reported in the literature were approximately ten times greater than of normal breast tissue 61,62 . www.nature.com/scientificreports www.nature.com/scientificreports/ The arterial blood temperatures of the skin, breast, and cancerous tissues in the model (T A S , , T A B , , and T A C , , respectively) were considered variable and assumed to be equal (see Eq. 5) based on the physiologic assumption that the same blood temperature supplied all three tissues. Arterial blood temperatures ranging between 33.0-35.5°C were considered, based on temperature measurements of the internal mammary artery reported in the literature 7,63,64 . Arterial temperatures, along with the thermophysical properties of blood outlined in Table 3, were incorporated into the UDF.
Finally, the perfusion rate of skin was considered constant based on values reported in the literature (see Table 3). The metabolic heat generation of skin was neglected due to its low volume relative to the underlying breast tissue and the lack of data available in the literature.

Parametric formulation and target metrics for model calibration.
In this study, ′′′ Q m and ω  b of normal and cancerous breast tissues were modeled as adjustable variables in order to determine the ( ω ′′′Q , m b ) solution pairs that recreated the surface temperatures observed on the subject's IR images. The subject's right breast was used to determine the thermal characteristics of normal breast tissue in the right normal breast model (see Fig. 5a), whereas the subject's left breast was used to determine the thermal characteristics of cancerous and non-cancerous breast tissues in the left cancerous breast model (see Fig. 5b). The ranges of adjustable variables assumed per tissue in each model are outlined in Table 1. In contrast, the properties of skin tissue in each model were kept constant for all solutions (see Table 3).
For each breast model, the ranges of blood perfusion rate were constrained to values reported by Delille et al. for normal and cancerous breast tissues as outlined in Table 3 62 . Each ω  b range was held constant for all T A and HTC cases simulated. Next, a corresponding range of metabolic heat generation rate was assigned to each breast tissue, which varied by T A case. A linearly proportional relationship was assumed between ′′′ Q m and ω  b ; that is, a low ′′′ Q m corresponded to a low ω  b and a high ′′′ Q m corresponded to a high ω  b . This physiologic modeling assumption was supported by recent blood oxygen level-dependent (BOLD) resting-state functional magnetic resonance imaging (rs-fMRI) observations regarding the increased metabolic and vascular activity in brain tumors 65,66 . Solutions for each tissue were found within the ( ω ′′′Q , m b ) ranges selected for each T A series (i.e., solutions lied along the solid T A lines illustrated in Fig. 1). The ′′′ Q m relation reported by Gautherie was not used in order to ensure the necessary degree of freedom in finding solutions.
Moreover, because breast tissue under normal physiologic conditions cannot physically exhibit higher temperatures than the posterior wall, a constraint on the maximum internal breast temperature was enforced for the right breast model (see Eq. 6). This constraint ensured the exclusion of extraneous solutions yielding artificially high internal temperatures for the right normal breast. This constraint was not enforced in the cancerous breast model due to the underlying thermophysiologic assumption that cancerous tissue (i.e., tumor) represented an internal heat source for the subject. Finally, owing to the non-homogeneous vascularization possible within a tumor 67,68 , the tumor was modeled with two regions (upper and lower portions) of varying metabolic heat generation (and proportional blood perfusion) in order to best match clinical data. Figure 5c demonstrates the upper and lower portions of the cancerous tissue, which was partitioned at the approximate axial midpoint between the two regions of skin retraction outlined earlier. Thermal characteristics of the lower cancerous tissue governed the subject's nipple hotspot, whereas the upper cancerous tissue governed the subject's superior hotspot. This heterogeneity in tumor metabolic heat generation and perfusion provided an additional degree of freedom in finding solutions.
Therefore, in all, solution pairs of ( ω ′′′Q , m b ) were determined for a total of four tissues: Normal breast tissue (right breast model), non-cancerous breast tissue (left breast model), upper cancerous breast tissue (left breast model), and lower cancerous breast tissue (left breast model). Solution pairs were found for various HTC cases by marching along each T A line, solving the model with the specified parameters, analyzing resulting surface temperatures, and comparing with the subject's IR images. For the right normal breast model, the mean surface temperature over the entire breast (33.8°C) was the sole metric for finding solutions, given its constraint (see Eq. 6). Here, ( ω ′′′Q , m b ) solutions of normal breast tissue governed the mean surface temperature. In contrast, for the left malig-