Computational Modeling of Deep Tissue Heating by an Automatic Thermal Massage Bed: Predicting the Effects on Circulation

Automatic thermal and mechanical massage beds support self-managed treatment, including reduction of pain and stress, enhanced circulation, and improved mobility. As the devices become more sophisticated (increasing the degrees of freedom), it is essential to identify the settings that best target the desired tissue. To that end, we developed an MRI-derived model of the lower back and simulated the physiological effects of a commercial thermal-mechanical massage bed. Here we specifically estimated the tissue temperature and increased circulation under steady-state conditions for typical thermal actuator settings (i.e., 45–65°C). Energy transfer across nine tissues was simulated with finite element modeling (FEM) and the resulting heating was coupled to blood flow with an empirically-guided model of temperature-dependent circulation. Our findings indicate that thermal massage increases tissue temperature by 3–8°C and 1–3°C at depths of 2 and 3 cm, respectively. Importantly, due to the rapid (non-linear) increase of circulation with local temperature, this is expected to increase blood flow four-fold (4x) at depths occupied by deep tissue and muscle. These predictions are consistent with prior clinical observations of therapeutic benefits derived from spinal thermal massage.

1. Introduction

Direct application of heat to the skin is known to increase blood flow (1, 2), an essential physiological response that also removes excess heat actively produced during metabolism and exercise. The increased circulation produced by passive heating has been exploited in multiple therapeutic settings, most notably in the management of pain (3). Several candidate mechanisms may mediate the therapeutic effect of heat therapy. The vasodilation and increased blood flow generated by thermal stimulation results in an increased supply of nutrients and oxygen to the affected tissue (4–6). Heat transduction is mediated by the TRP vanilloid 1 (TRPV1) receptor, whose activation may modulate antinociceptive pathways (7). Heating also increases the rate of local tissue metabolism (5) and the elasticity of connective tissue (8, 9), both of which may in turn promote pain relief. Understanding the distribution of tissue heating—the specific goal of the present study—is an essential first step in developing quantitative models of bio-thermomechanical responses.

Medical devices that provide thermal stimulation are becoming increasingly sophisticated. For example, automatic thermal-mechanical massage beds are equipped with multi-axis traction and far-infrared thermal projectors (10–12). Such devices allow in-home, self-managed treatment that may potentially leverage the numerous degrees of freedom governing the stimulation “dose.” In theory, stimulation parameters such as temporal dynamics (pulsed vs. continuous), intensity (temperature), duration, and position may be individualized to maximize therapeutic benefit. Understanding the relationship between these numerous parameters and the resulting physiological changes is critical to such optimization efforts.

A central question in thermal therapies pertains to target engagement: given a targeted region (e.g., a specified muscle group), what is the required temperature that must be applied to the surface in order to produce the desired increase in blood flow to the target? Addressing this question requires knowledge of the temperature distribution produced throughout the volume during passive heating, as well as the quantitative relationship between the temperature distribution and resulting change in circulation. Recent work indicates that the increase in circulation during passive heating is largely driven by local temperature gradients (13), implying that the magnitude of circulation increase may be predicted with knowledge of the local temperature. To estimate the temperature distribution produced by a given thermal device, computational modeling that accounts for the heterogeneous anatomical structures of the treated area may be leveraged. Indeed, computer-based simulations are standard tools across medical device design (such as neuromodulation) and increasingly rely on detailed image (e.g., MRI) derived representation of the underlying anatomy (14–17).

Here we performed a realistic simulation of tissue heating and associated circulation changes produced by a commercial thermo-mechanical massage bed. We optimized MRI imaging to resolve key structures of the back, developed a 3D model of the relevant tissues, modeled the action of thermal actuators under static conditions, and employed finite element modeling to predict the intensity and distribution of temperature in the lumbar region. We then leveraged prior work relating local temperature to circulation to construct a model of blood flow changes during thermal massage. Our findings indicate that temperature increases up to 3°C are achievable at deep (3 cm) regions of the lumbar back with conventional actuator settings. Moreover, our findings suggest that blood flow can be increased four-fold at depths of 2–3 cm, corresponding to the lumbar musculature. These findings are consistent with clinical observations performed with the same thermo-mechanical massage bed (10–12, 18).

2. Methods

2.1. MRI Scanning and Tissue Segmentation

We performed T2-weighted lumbar spine magnetic resonance imaging (MRI) of a healthy male adult (BMI: 25 kg/m²; age: 42 years) with a 3T Siemens MAGNETOM Prisma scanner equipped with a CP Spine array coil (Siemens Healthineers, PA, USA). A sagittal slice is shown in Figure 1A. The parameters of image acquisition were set according to: TE = 99 ms, TR = 7,040 ms, flip angle = 130°, FOV = 256 mm, SNR = 1, in-plane resolution = 256 x 256, slice thickness = 1 mm; and voxel size: 1 x 1 x 1 mm.

FIGURE 1

Figure 1. Constructing a computational model to predict tissue temperature and blood flow during thermal massage. (A) Sagittal slice of a T2-weighted anatomical MRI used to construct the computational model developed in this study. (B) Based on image contrast, the volume was segmented into nine tissues that were then endowed with physical properties governing the physiological response to an external heat source. (C) The Pennes bio-heat equation was solved with a finite element model (FEM) solver, yielding the temperature during thermal massage throughout the volume. The locations of the four actuators are apparent from the temperature “hot-spots” in the superficial slice. (D) In order to translate model derived tissue temperature to the predicted changes in blood flow, we constructed an empirically-guided sigmoidal model that outputs the expected blood flow based on in situ tissue temperature. Red markers indicate the measurements performed by Chiesa et al., while the gray curve denotes the resulting model fit.

Based on image intensity, the MRI was segmented into nine tissue masks: skin, subcutaneous fat, soft tissue, muscles, intervertebral disc, vertebrae, epidural fat, cerebrospinal fluid, and spinal cord. Manual tissue segmentation was carried out by first thresholding the images, followed by morphological filtering such as flood fill, dilation, and erosion, performed in Simpleware ScanIP (Synopsys Inc., CA, USA). To ensure tissue continuity and improve segmentation accuracy, the data was extensively visualized and compared with the anatomy of the subject's spine while performing manual adjustments to the tissue masks. Based on the results of segmentation, the thicknesses of the tissues comprising the imaged region were measured as: (skin) 1.1 mm, (subcutaneous fat) 13 mm, (muscle) 42 mm, (soft tissue) 2.8 mm, (epidural fat) 2.2 mm, (CSF) 2.9 mm, and (spinal cord) 1.9 mm. Thus, the muscle (i.e., the presumed target of thermal massage) occupied a region approximately 1.5–5.5 cm from the surface.

Note that hair is not resolvable with MRI, and was thus not captured by the computational model constructed in our study. Due to the fact that the stimulation was applied to the lumbar back region, this omission was considered to have minimal impact on the resulting temperature and blood flow distribution.

2.2. Finite Element Modeling

The tissue masks generated by the segmentation process were utilized to construct a finite element model (FEM) of the lumbar back. To that end, we employed the Simpleware Scan IP software running the tetrahedral voxel-based meshing algorithm to generate a FEM consisting of more than 3.38 million tetrahedral elements. The resulting model was then integrated with separate domains that captured the position, orientation, and geometry of a thermal massage device. In particular, we modeled the characteristics of the CGM MB-1901 massage bed (CERAGEM Co. Ltd., Cheonan, Korea), which is equipped with two types of heat sources: (i) a contoured heating “mat” that resides directly under the bed surface and provides a diffuse layer of heat, and (ii) four actuators (“rollers”) that make direct contact with the body to provide more punctate and intense thermal stimulation. The device contains four actuators arranged as the vertices of a rectangle: the lateral distance between actuators is 56 mm, and the vertical distance between top and bottom actuators is 32 mm. Specifically, we modeled cross-sections of the four actuators (width: 65 mm, diameter of the circular ends: 45 mm, spacing between the two circular ends: 30 mm) in Solidworks (Dassault Systems, MA, USA) and imported the components into Simpleware ScanIP for positioning and meshing. The resulting model was then imported into COMSOL Multiphysics 5.5 (COMSOL Inc., MA, USA) to solve the Pennes' bio-heat transfer equation governing heating during thermal massage as:

$\begin{array}{ll}\rho c\frac{\partial T}{\partial t}=\nabla ·(k\nabla T)+{\rho }_b{c}_b{w}_b(T-T_b)+Q_m,& (1)\end{array}$

where ρ is the tissue density (units of kg/m³), c is the specific heat capacity (J kg⁻¹ K⁻¹), T is the temperature (K), t denotes time (s), k is the thermal conductivity (W m⁻¹ K⁻¹), ρ_b is the density of blood (kg m⁻³), c_b is the specific heat capacity of blood (J kg⁻¹ K⁻¹), w_b is the perfusion rate of blood (1/s), T_b is blood temperature (K), and Q_m is the rate of metabolic heat generation (W/m³). The solution was carried out under steady state conditions.

The physical and thermal properties of the model tissues were assigned according to the values listed in Table 1, and were obtained from a survey of prior literature (19–25). In order to solve the bio-heat equation (1), the temperature at the model's external boundaries was fixed to core body temperature (37°C). We assumed the absence of convective gradients across all internal tissue boundaries. Convective heat loss from model boundaries to the environment was modeled as:

$\begin{array}{l}{q}_0=h(T_{\text{amb}}-T)\end{array}$

where q₀ (units of W m⁻²) is the convective heat flux, h = 5 W m⁻² K is the heat transfer coefficient for the tissues comprising the model, and T_amb is the ambient temperature (i.e., 25 °C). The initial temperature of the tissues was set to 37°C. The temperature of the actuators was fixed within a simulation to a value ranging from 45 to 65 °C. The external boundary of the skin making contact with the bed surface, corresponding to the heating mat, was assigned a temperature of 40°C.

TABLE 1

Table 1. Thermal properties of the tissues and sources comprising the computational model developed here to estimate temperature gradients during thermal massage.

2.3. Model of Temperature-Dependent Circulation

In order to estimate the magnitude and distribution of blood flow produced during thermal massage, we developed a non-linear model for mapping the in situ tissue temperature to the resulting blood flow. The model was guided by simultaneous empirical measurements of temperature and blood flow in the human leg during heat stress as conducted by Chiesa and colleagues (13). The general form of the model developed here was given by:

$\begin{array}{ll}F(\Delta T)=F_o+(F_{\text{max}}-F_o)e^{-(a\Delta T+b)},& (2)\end{array}$

where F is the blood flow (L/min), ΔT is the change in temperature at the site of the vessel (°C), F_o is the blood flow in the absence of exogenous heating, F_max is the maximum achievable blood flow due to the physical limitations of the cardiovascular system, and the free model parameters a and b were estimated by fitting the model to the empirical measurements reported by Chiesa and colleagues. Namely, the following (T, F) pairs of measurements were reported, where T was taken here to be the temperature in the muscle of the leg (i.e., the deepest measurement taken and the one closest to the site of the vessel) and F was taken to be the empirical blood flow in the central femoral artery (CFA): (34.9°C, 0.31 L/min), (36.2°C, 0.53 L/min), (37.0°C, 0.86 L/min), (37.6°C, 1.24 L/min), and (38.3°C, 1.22 L/min).

From these empirical measurements, the values F_o = 0.31 L/min and F_max = 1.22 L/min were assigned to the baseline and maximum blood flow, respectively. We then performed conventional least-squares fitting employing the Matlab (Mathworks, MA, USA) function lsqcurvefit. The minimization procedure yielded values of a = −2.67 and b = 4.88, producing excellent fits to the empirical measurements while capturing a sigmoidal transition from baseline to maximal flow (Figure 1D). The resulting model was then employed to transform local temperature to the corresponding blood flow.

3. Results

We simulated the action of a commercial automatic thermal massage device, where contact heating was delivered through (i) a contoured plane (“mat”) at 40°C covering the full extent of the lower back, and (ii) four actuators (“rollers”) with a diameter of 4.5 cm and a height of 1 cm. The actuators were arranged on the back as the vertices of a rectangle (Figures 1B,C). The actuator temperature was varied from 45 to 65 °C in increments of 5°C. For each intensity, the temperature distribution throughout the lumbar region was numerically computed. The resulting solution was then coupled to the predicted change in blood flow. Our objective was to determine the magnitude of the circulation increase as a function of tissue depth and actuator temperature.

3.1. Spatial Distribution of Lower Back Temperature During Thermal Massage

The spatial distribution of temperature at a sagittal slice immediately anterior to a pair of actuators is shown in Figure 2A (shown here for an actuator temperature of 65°C), while the contours of the temperature field are depicted in Figure 2B. The spatial dynamics of the temperature increase are apparent, and it is clear that the largest changes occur near the stimulation site (note the higher density of contour lines). To quantify the effect of tissue depth on achieved temperature, we analyzed temperature immediately anterior to the centroid of an actuator, as a function of depth. The tissue temperature exhibited an exponential decay with depth (Figure 2E, colors indicate the input temperature, markers indicate depths of 1, 2, 3, and 4 cm). With an actuator temperature of 45 °C, the temperature at depths of 1, 2, and 3 cm into the tissue was 42.2, 39.6, and 38.0 °C, respectively (Figure 2E, blue). With an increased roller temperature of 65°C, the corresponding temperatures were 54.9, 45.5, and 40.6°C (Figure 2E, green). Thus, the temperature gradient at a depth of 2 cm ranged from 3 to 8°C, while the gradient at 3 cm was 1–3°C.

FIGURE 2

Figure 2. Quantifying temperature increases during thermal massage. (A) Temperature of a sagittal slice of the volume during thermal massage with the actuators set to 65°C. The selected slice was located directly anterior to the actuators. The resulting temperature decayed exponentially with distance from the heating elements. (B) The contours of the temperature distribution shown in (A). (C) Temperature of a sagittal slice located between actuators. The resulting temperature increase was markedly lower. (D) Contours corresponding to the temperature distribution in (C). (E) The temperature directly anterior to an actuator as a function of depth, shown for five different actuator temperatures. At a depth of 2 cm, temperature increased by 3–8°C. The temperature increase at 3 cm depth was between 1 and 3°C. (F) Between actuators, the achieved temperature increase was greatly dampened, with an increase of only 1°C at a depth of 2 cm. (G) The tissue temperature distribution from left to right, shown for a fixed depth of 2 cm. The proximity of the tissue to the actuator placement determines the level of heating (actuators are positioned at ±6 cm). (H) Same as (G) but now for a depth of 3 cm, showing the increment in deep tissue temperature for different actuator settings.

Next we considered the effect of passive heating on tissue that is laterally displaced from the contact site. Figure 2C depicts the spatial distribution at a sagittal slice positioned at the midline (i.e., between pairs of actuators, 6 cm from the actuator centroid and 3.75 cm from the most lateral edge). At an actuator temperature of 45°C, tissue temperatures were computed as 39.1, 37.9, and 37.5°C at depths of 1, 2, and 3 cm, respectively (Figure 2F, blue). The corresponding temperatures with an actuator temperature of 65°C were only marginally higher: 39.3, 38.3, and 37.8°C (Figure 2F, green), indicating that increasing source intensity does not compensate for a lateral displacement of the target.

To better characterize the effect of thermal massage on tissue temperature, we also analyzed the temperature across the full lateral extent of the back, but at fixed depths of 2 cm (Figure 2G) and 3 cm (Figure 2H). We considered a vertical location matching one (horizontal) pair of actuators. The resulting temperature profile confirmed the preferentially local action of the actuators, as the temperature fell to near mean body temperatures at the midline (37.9°C at a depth of 2 cm, and 37.5°C at a depth of 3 cm, Figures 2G,H). With an actuator temperature of 45°C evaluated at a depth of 2 cm, the temperature fell from a peak of 39.3–38.9 °C with a lateral displacement of 1 cm toward the midline, to 37.9 °C with a displacement of 2 cm, and to 36.4 °C with a displacement of 3 cm (Figure 2G, solid blue). With an elevated actuator temperature of 65°C, the corresponding values were 44.6°C (peak), 43.8°C (1 cm displacement), 41.6°C (2 cm), and 38.3°C (3 cm; Figure 2G, solid green). Evaluated at a depth of 3 cm, the peak temperature with actuators set to 45°C was 38.0°C, which fell to 37.7, 37.1, and 35.2°C when laterally displaced by 1, 2, and 3 cm, respectively. When simulating an actuator temperature of 65°C, lateral displacements of 1, 2, and 3 cm reduced a peak temperature of 40.5–40.2°C, 39.2 and 36.7°C, respectively (Figure 2H).

3.2. Predicted Blood Flow Increases in the Lower Back During Thermal Massage

In order to predict the magnitude of increased blood flow during thermal massage, we developed a simple non-linear model of the relationship between in situ tissue temperature and blood flow. The model was based on concurrent empirical measurements of temperature and blood flow in the human leg during heat stress, where flow was found to non-linearly increase from a baseline value of 0.31 L/min to a maximum value of 1.22 L/min (Figure 1D; see Methods for details of the empirical data and model fitting procedure). We employed the resulting sigmoidal model to transform the computationally derived tissue temperature values (shown in Figure 2) to estimates of blood flow in the lower back during thermal massage.

The spatial distribution of predicted blood flow is shown in Figures 3A,B for a sagittal slice immediately anterior to a (vertical) pair of actuators, and in Figures 3C,D for a mid-sagittal slice (both shown for an actuator temperature of 65°C). Upon comparing Figure 3A to Figure 2A, it is evident that blood flow exhibits a more favorable decay with depth.

FIGURE 3

Figure 3. Predicting blood flow increases during thermal massage. (A) The predicted blood flow at a sagittal slice of the back, immediately anterior to the actuator. The actuator temperature was set to 65°C. Compared to temperature (Figure 2), pronounced increases in blood flow are observed at greater depths. (B) Contours of the blood flow distribution in (A), where the concentration of the contours indicates a sharp transition region from high to low blood flow. (C) The predicted blood flow at a sagittal slice located between actuators. The increase in blood flow is visibly dampened. (D) Contours of the blood flow distribution in (C), where the transition region is observed at a relatively shallow depth compared to (B). (E) The blood flow achieved immediately anterior to the actuators, shown for five different actuator intensities. Due to the sigmoidal nature of the temperature-circulation relationship, the falloff with depth occurs at a larger depth (between 2-4 cm, depending on actuator setting). Four-fold increases in circulation are predicted at depths of 2–3 cm. (F) Same as (E) but now shown for a sagittal slice located between the actuators. The increase in blood flow is no longer apparent beyond a depth of 3–4 cm. (G) The distribution of blood flow as a function of left-to-right position, shown for a fixed 2 cm depth. The predicted blood flow saturates (4x increase relative to baseline) across large sections of the back. (H) Same as (G) but now shown for a depth of 3 cm. When located directly anterior to the actuator, the predicted blood flow reaches a three-fold increase at an actuator temperature of 55°C.

At an input temperature of 45°C, the predicted blood flow was 1.22 L/min, 1.11 L/min, 0.41 L/min, and 0.33 L/min at depths of 1, 2, 3, and 4 cm, respectively (Figure 3E, blue line). With an elevated input temperature of 65°C, the predicted flow saturated at depths of up to 3 cm (1.22 L/min), before sharply falling to 0.56 L/min at a depth of 4 cm (Figure 3E, blue line). It is evident that the nonlinearity of the relationship between temperature and blood flow leads to abrupt transitions between low, baseline flow and large, saturated flow. Indeed, the sharp sigmoidal transition (Figure 1D) predicts that large increases in blood flow are achievable at depths of 3 cm, where a flow of 1.03 L/min is predicted with an actuator temperature of 55 °C (Figure 3E, yellow line).

When considering blood flow as a function of lateral displacement from the center of the actuator, we observed that, given a sufficiently high input temperature (i.e., 55 °C), the blood flow at a depth of 2 cm saturated across a large lateral extent of the lower back (Figure 3G). Even at a depth of 3 cm, the predicted flow produced by an actuator temperature of 55°C reached a three-fold increase at locations immediately anterior of the actuators (Figure 3G). Overall, the effective sensitivity of blood flow to distance from the heat source is less than that observed with temperature.

4. Discussion

We employed finite element modeling to estimate the physiological effects of an automatic thermal massage bed on the lumbar back region. The model predicts that it is feasible to produce large (e.g., four-fold) increases in circulation at relatively deep locations with high but practical actuator temperatures. The basis of this finding is that blood flow increases quickly and saturates at in situ temperature gradients of approximately 3° (13). Specifically, our FEM model predicts that such temperature increases are attainable at depths of 2 to 3 cm with actuator temperatures between 55° and 65°, values that overlap with those tested in earlier studies (11, 12). Thus, despite the exponential decay of temperature with depth (Figure 2), there is a sufficient gradient produced at the depths occupied by muscle tissue (1.5–5.5 cm) to modulate blood flow to the muscles of the lumbar region, which represent a common target of therapeutic massage beds.

To our knowledge, this study represents the first computational model of the spatial distribution of tissue temperature changes induced by automatic heat therapy devices. Detailed segmentation of tissue compartments resolved by MRI facilitated the construction of the model, which captured the idiosyncrasies of the lumbar region. By then solving the biophysical equations governing heat transfer in tissue in response to an exogenous source, we were able to quantify the depth of heat penetration, as well as the associated modulation of blood flow. Notwithstanding the experimental data relating local temperature to blood flow (13), and the general recognition that heating can increase perfusion, there remains an absence of empirical data describing the spatial distribution of temperature and blood flow produced by automatic massage beds such as the one modeled here. Consequently, empirical validation of the predictions of the computational model developed here represents a key focus for future studies. For the device parameters simulated here, heating was applied at comfortable levels. However the present methodology can be similarly deployed to assess the likelihood of tissue injury in the case of much higher heat source temperatures.

Increasing the temperature of the blood is expected to reduce its viscosity (26, 27). Interestingly, the reduced viscosity may in turn promote an increase of circulation. In this study, such a mediating effect of blood viscosity was accounted for by the employment of experimental data that aggregated all potential mechanisms.

Thermotherapy (3, 28) and massage (29, 30) have an extensive history in wellness and medicine but recent advances in automatic massage beds suggest new potency for self-managed treatment. Evidently, the therapeutic outcomes of thermal-mechanical massage will depend on the technical features of the mechanical and thermal actuators, as well as (individualized) anatomy and physiology. Notice, for example, that the position of tissues relative to the actuators had a strong influence on the resulting temperature and blood flow changes (Figures 2, 3). Knowledge of the distance and orientation of the target relative to the actuator(s), which will generally vary with the body composition of the individual, may be leveraged by future devices to personalize treatment. For example, the actuator temperature and treatment duration required to achieve a certain benchmark circulation increase may be determined automatically, potentially systematizing treatment.

Computational models of medical devices relate device features (exogenous and set by the operator) to the resulting physiological changes, which in turn govern therapeutic actions. Therefore, the effects on heating and circulation predicted here have direct implications for understanding (and further optimizing) results from clinical trials using the same device. As both temperature and blood flow are known to impact immune function, antioxidant and anti-inflammatory processes, and autonomic function (31–36), these predictions help to explain the clinical biomarkers observed during application of the device modeled here (10–12). In turn, these biomarkers provides a mechanistic substrate for thermal massage therapy in pain (10, 12).

As with any model, our predictions are subject to assumptions, including the transfer function employed here to relate temperature increases to blood flow. Valuable next steps will include (1) directly verifying model predictions by physiological measurements conducted either during or after thermal massage; and (2) designing and validating improvements in automatic massage bed operation. The present study assumed a static (stationary) regime when solving for the temperature during thermal massage. The rationale for this assumption is that the typical duration of thermal massage (e.g., tens of minutes) is expected to be larger than the time constants defining the transients of the bio-heat equation (37). In any case, the predictions made here correspond to the physiological effects produced when transients have passed and steady-state values are achieved. Nevertheless, the addition of a dynamic component to the simulation may allow one to identify potential advantages of pulsed stimulation, if for example, the vascular response is reinforced by a time-varying heat source marked by frequent onsets. Moreover, the mechanical deformations that occur in the lumbar region during thermal-mechanical massage may bring target tissues closer to the heat source during certain periods of the treatment, potentially increasing the improvement in circulation that was measured in a static configuration here.

Future modeling efforts should therefore integrate the mechanical effects of massage with the simultaneous temperature changes produced by thermotherapy in a multi-physics simulation that captures the potentially synergistic property of thermal massage. This may reveal that thermal stimulation modulates mechanical tissue properties and promotes the beneficial effects of tissue manipulation. Such symbiotic outcomes of thermal-mechanical massage are expected to translate to improved clinical outcomes, and will rely on computational modeling to design therapies that optimize the interplay between heating and mechanical forces.

Data Availability Statement

The datasets presented in this article are not readily available because of the proprietary nature of the data. Requests to access the datasets should be directed to jdmochowski@ccny.cuny.edu.

Ethics Statement

This study was not required to obtain ethical review and approval, as the study utilized only de-identified imaging data that was previously collected for another study.

Author Contributions

JD and NK designed the research, performed the research, and wrote the manuscript. LC designed the research. EM performed the research. KL, SK, YJ, and MB designed the research and wrote the manuscript. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The City University of New York holds patents on brain stimulation with MB as inventor. MB has equity in Soterix Medical Inc. MB consults, received grants, assigned inventions, and/or serves on the SAB of SafeToddles, Boston Scientific, GlaxoSmithKline, Biovisics, Mecta, Lumenis, Halo Neuroscience, Google-X, i-Lumen, Humm, Allergan (Abbvie), and Apple. KL, SK, and YJ are employed by Ceragem Clinical Inc. This study received funding from Ceragem Clinical Inc. in the form of a grant to MB, LC, and JD. The funder had the following involvement with the study: providing information about the parameters of the medical device being simulated; providing feedback on simulation results.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Rowell LB, Brengelmann GL, Murray JA. Cardiovascular responses to sustained high skin temperature in resting man. J Appl Physiol. (1969) 27:673–80. doi: 10.1152/jappl.1969.27.5.673

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Charkoudian N. Skin blood flow in adult human thermoregulation: how it works, when it does not, and why. In: Mayo Clinic Proceedings. vol. 78. Rochester, MN: Elsevier (2003). p. 603–12.

PubMed Abstract | Google Scholar

3. Malanga GA, Yan N, Stark J. Mechanisms and efficacy of heat and cold therapies for musculoskeletal injury. Postgrad Med. (2015) 127:57–65. doi: 10.1080/00325481.2015.992719

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Rabkin JM, Hunt TK. Local heat increases blood flow and oxygen tension in wounds. Arch Surg. (1987) 122:221–5. doi: 10.1001/archsurg.1987.01400140103014

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Nadler SF, Weingand K, Kruse RJ. The physiologic basis and clinical applications of cryotherapy and thermotherapy for the pain practitioner. Pain Physician. (2004) 7:395–400. doi: 10.36076/ppj.2004/7/395

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Petrofsky JS, Lawson D, Suh HJ, Rossi C, Zapata K, Broadwell E, et al. The influence of local versus global heat on the healing of chronic wounds in patients with diabetes. Diabet Technol Therapeut. (2007) 9:535–44. doi: 10.1089/dia.2007.0231

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Palazzo E, Rossi F, Maione S. Role of TRPV1 receptors in descending modulation of pain. Mol Cell Endocrinol. (2008) 286:S79-S83. doi: 10.1016/j.mce.2008.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hardy M, Woodall W. Therapeutic effects of heat, cold, and stretch on connective tissue. Journal of Hand Therapy. (1998) 11:148–56. doi: 10.1016/S0894-1130(98)80013-6

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Bleakley CM, Costello JT. Do thermal agents affect range of movement and mechanical properties in soft tissues? A systematic review. Arch Phys Med Rehabil. (2013) 94:149–63. doi: 10.1016/j.apmr.2012.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kim KE, Shin NR, Park SH, Nam SY, Yoon YS, Park SK, et al. Modulation of the human immune status by spinal thermal massage: a non-randomized controlled study. Signa Vitae. (2021) 1:10. doi: 10.22514/sv.2021.144

CrossRef Full Text | Google Scholar

11. Lee YH, Park BNR, Kim SH. The effects of heat and massage application on autonomic nervous system. Yonsei Med J. (2011) 52:982–9. doi: 10.3349/ymj.2011.52.6.982

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Kim KE, Park JS, Cho IY, Yoon YS, Park SK, Nam SY. Use of a spinal thermal massage device for anti-oxidative function and pain alleviation. Front Public Health. (2020) 8:493. doi: 10.3389/fpubh.2020.00493

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Chiesa ST, Trangmar SJ, Kalsi KK, Rakobowchuk M, Banker DS, Lotlikar MD, et al. Local temperature-sensitive mechanisms are important mediators of limb tissue hyperemia in the heat-stressed human at rest and during small muscle mass exercise. Am J Physiol Heart Circ Physiol. (2015) 15:H369–80. doi: 10.1152/ajpheart.00078.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Datta A, Bansal V, Diaz J, Patel J, Reato D, Bikson M. Gyri-precise head model of transcranial direct current stimulation: improved spatial focality using a ring electrode versus conventional rectangular pad. Brain Stimul. (2009) 2:201–7. doi: 10.1016/j.brs.2009.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

15. De Paolis A, Watanabe H, Nelson JT, Bikson M, Packer M, Cardoso L. Human cochlear hydrodynamics: a high-resolution μCT-based finite element study. J Biomech. (2017) 50:209–16. doi: 10.1016/j.jbiomech.2016.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kreisberg E, Esmaeilpour Z, Adair D, Khadka N, Datta A, Badran BW, et al. High-resolution computational modeling of the current flow in the outer ear during transcutaneous auricular Vagus Nerve Stimulation (taVNS). Brain Stimul. (2021) 14:1419–30. doi: 10.1016/j.brs.2021.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Myoung HS, Kim DH, Kim HS, Lee KJ. Design of a stimulation protocol to predict temperature distribution in subcutaneous tissue using the finite element model. Biomed Eng Lett. (2017) 7:261–6. doi: 10.1007/s13534-017-0029-0

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Choi JH, Kim ES, Yoon YS, Kim KE, Lee MH, Jang HY. Self-management techniques and subsequent changes in pain and function in patients with chronic low back pain. J Digit Convergence. (2020) 18:547–55. doi: 10.14400/JDC.2020.18.10.547

CrossRef Full Text | Google Scholar

19. Wilson SB, Spence VA. A tissue heat transfer model for relating dynamic skin temperature changes to physiological parameters. Phys Med Biol. (1988) 33:895. doi: 10.1088/0031-9155/33/8/001

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Hodson D, Barbenel J, Eason G. Modelling transient heat transfer through the skin and a contact material. Phys Med Biol. (1989) 34:1493. doi: 10.1088/0031-9155/34/10/011

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Çetingül MP, Herman C. A heat transfer model of skin tissue for the detection of lesions: sensitivity analysis. Phys Med Biol. (2010) 55:5933. doi: 10.1088/0031-9155/55/19/020

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Mcintosh RL, Anderson V. A comprehensive tissue properties database provided for the thermal assessment of a human at rest. Biophys Rev Lett. (2010) 5:129–151. doi: 10.1142/S1793048010001184

CrossRef Full Text | Google Scholar

23. Collins CM, Smith MB, Turner R. Model of local temperature changes in brain upon functional activation. J Appl Physiol. (2004) 97:2051–5. doi: 10.1152/japplphysiol.00626.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zannou AL, Khadka N, Truong DQ, Zhang T, Esteller R, Hershey B, et al. Temperature increases by kilohertz frequency spinal cord stimulation. Brain Stimul. (2019) 12:62–72. doi: 10.1016/j.brs.2018.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zannou AL, Khadka N, FallahRad M, Truong DQ, Kopell BH, Bikson M. Tissue temperature increases by a 10 kHz spinal cord stimulation system: phantom and bioheat model. Neuromodulation. (2021) 24:1327–35. doi: 10.1111/ner.12980

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Virgilio RW, Long DM, Mundth ED, McClenathan JE. The effect of temperature and hematocrit on the viscosity of blood. Surgery. (1964) 55:825–30.

PubMed Abstract | Google Scholar

27. Snyder GK. Influence of temperature and hematocrit on blood viscosity. Am J Physiol Legacy Content. (1971) 220:1667–72. doi: 10.1152/ajplegacy.1971.220.6.1667

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Brosseau L, Yonge K, Welch V, Marchand S, Judd M, Wells GA, et al. Thermotherapy for treatment of osteoarthritis. Cochrane Database Syst Rev. (2003) 2003:CD004522. doi: 10.1002/14651858.CD004522

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kumar S, Beaton K, Hughes T. The effectiveness of massage therapy for the treatment of nonspecific low back pain: a systematic review of systematic reviews. Int J Gen Med. (2013) 6:733. doi: 10.2147/IJGM.S50243

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Field T. Massage therapy research review. Complement Ther Clin Pract. (2014) 20:224–9. doi: 10.1016/j.ctcp.2014.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Johnson RH. Disorders of the autonomic nervous system. Anat Auton Nervous Syst. (1974) 1–22.

32. Bergholm R, Mäkimattila S, Valkonen M, Liu Ml, Lahdenperä S, Taskinen MR, et al. Intense physical training decreases circulating antioxidants and endothelium-dependent vasodilatation in vivo. Atherosclerosis. (1999) 145:341–9. doi: 10.1016/S0021-9150(99)00089-1

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Mestre-Alfaro A, Ferrer MD, Banquells M, Riera J, Drobnic F, Sureda A, et al. Body temperature modulates the antioxidant and acute immune responses to exercise. Free Radic Res. (2012) 46:799–808. doi: 10.3109/10715762.2012.680193

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Evans SS, Repasky EA, Fisher DT. Fever and the thermal regulation of immunity: the immune system feels the heat. Nat Rev Immunol. (2015) 15:335–49. doi: 10.1038/nri3843

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Pober JS, Sessa WC. Inflammation and the blood microvascular system. Cold Spring Harb Perspect Biol. (2015) 7:a016345. doi: 10.1101/cshperspect.a016345

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Sheng Y, Zhu L. The crosstalk between autonomic nervous system and blood vessels. Int J Physiol Pathophysiol Pharmacol. (2018) 10:17.

PubMed Abstract | Google Scholar

37. Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. (1948) 1:93–122. doi: 10.1152/jappl.1948.1.2.93

PubMed Abstract | CrossRef Full Text | Google Scholar