A Concise Model and Analysis for Heat-Induced Withdrawal Reflex Caused by Millimeter Wave Radiation

In this study, we consider the heat-induced withdrawal reflex caused by exposure to an electromagnetic beam. We propose a concise dose-response relation for predicting the occurrence of withdrawal reflex from a given spatial temperature profile. Our model is distilled from sub-step components in the ADT CHEETEH-E model developed at the Institute for Defense Analyses. Our model has only two parameters: the activation temperature of nociceptors and the critical threshold on the activated volume. When the spatial temperature profile is measurable, the two parameters can be determined from test data. We connect this dose-response relation to a temperature evolution model for electromagnetic heating. The resulting composite model governs the process from the electromagnetic beam deposited on the skin to the binary outcome of subject’s reflex response. We carry out non-dimensionalization in the time evolution model. The temperature solution of the non-dimensional system is the product of the applied power density and a parameter-free function. The effects of physical parameters are contained in non-dimensional time and depth. Scaling the physical temperature distribution into a parameter-free function greatly simplifies the analytical solution, and helps to pinpoint the effects of beam spot area and applied power density. With this formulation, we study the theoretical behaviors of the system, including the time of reflex, effect of heat conduction, biological latency in observed reflex, energy consumption by the time of reflex, and the strategy of selecting test conditions in experiments for the purpose of inferring model parameters from test data.


Introduction
Millimeter wave (MMW) is a subset of radio frequency (RF) in the 30 -300 gigahertz (GHz) range. Human body exposure to MMW radiation at sufficiently high intensities increases the skin temperature and induces thermal pain. Since the energy penetration depth of MMW irradiation in tissue is very shallow (millimeter or less), the primary effect of MMW exposure is temperature increase near the skin surface. Due to the rapid growth in using MMW in common applications including wireless communications systems, automobile collision avoidance systems, airport security screening, non-lethal crowd control weapons, medical imaging, and detection of vital signs such as respiration and heartbeat rates [1] [2] [3], it is important to understand the thermal responses of humans to MMW irradiation and to evaluate biological safety with respect to thermal hazards.
Considerable work has been done in MMW interactions with the human body [4]. We briefly review a few papers that motivate us for this study. In [5], cutaneous thresholds for thermal pain were measured in 10 human subjects (3 female and 7 male Caucasian volunteers, 31 -70 years old with an average age of 43.7, military or DoD civilians or contractors). Each subject was tested with 3-second exposures to 94 GHz electromagnetic wave of intensity up to 1.8 W/cm 2 .
During each exposure, the temperature increase at the skin's surface was measured by infrared thermography. The mean baseline temperature of the skin was 34˚C. The irradiated area of the skin has a diameter of 4 cm. The threshold for pricking pain was 43.9˚C (at skin surface). The measured surface temperature was in good agreement with a simple one-dimensional thermal model that accounted for heat conduction and for the penetration depth of the electromagnetic energy into tissue. One important observation in [5] was that the skin surface temperature increased roughly by 10˚C after a 3-second exposure to 1.8 W/cm 2 at 94 GHz.
Heating of tissues by electromagnetic waves has also been studied in [6] using Pennes' bioheat equation, which includes the effect of the blood convection in dissipating the heat. The model was used to estimate the threshold for perception of pain as a function of frequency and exposure duration in the case of plane-wave irradiation. The numerical results in [6] agree well with the experimental thresholds for warmth perception evoked by electromagnetic waves of 2.45 -94 GHz applied to the back of the test subject for 10-second intervals.
In [7] [8], a one-dimensional multi-layer tissue model was used and solved with the finite difference method to predict temperature variations in the skin exposed to electromagnetic waves. This approach allows estimation of tempera-H. Y. Wang et al.
ture distribution and prediction of burn injury in different layers of the skin depending on blood perfusion rate, thermal conductivity, power density, and exposure time. Even though the numerical simulations were focused on electromagnetic waves with frequencies of 10 GHz and below, this multi-layer model can be extended directly to MMW frequencies. These studies concluded that the rate of skin surface temperature increase in humans in response to brief, high-power MMW exposures can be mathematically modeled using the one-dimensional thermal conduction model.
In [9], a spherical heterogeneous model was developed to simulate the thermal effects (surface and subsurface) on a primate (monkey) head exposed to far-field radiation at 100 GHz. It was found that temperature (both surface and subsurface) increases as the energy level increases, whereas high-power millimeter-waves (HPMs) with power densities in the range of 1.0 to 3.0 W/cm 2 for 3 seconds have a higher peak temperature than low-power millimeter-waves (LPMs) in the range of 0.1 to 0.3 W/cm 2 for 30 seconds. The surface temperature increase is generally linear with applied energy density for HPMs except under combined conditions of high blood-flow rate and high-energy density. In contrast, with LPMs, the surface temperatures do not behave linearly with respect to incident energy. The simulations also showed that the subsurface (i.e., mid-scalp and mid-skull) temperature increases are substantially damped compared to the surface (i.e., scalp) temperature.
Far-field exposure of the human face to a linearly polarized plane wave at frequencies from 6 to 100 GHz and with exposure durations of 100 milliseconds to 10 seconds was modeled in [10]. The Maxwell equations were solved using a finite-difference time-domain solver. The Pennes' bioheat transfer equation extended with a term representing the electromagnetic power absorption was used to model body temperature and it was numerically solved. The electromagnetic and thermal simulations revealed highly non-uniform frequency-dependent patterns of temperature rise.
According to a review published in 2016 [11], there are very limited data in the literature related to skin temperature rises or thermal sensation thresholds in humans exposed to MMWs. Most data involves brief exposures to MMWs around 94 GHz [12]. For example, experiments with 13 male adults were conducted in San Antonio, TX to quantify the thermal and behavioral effects of exposure to 30 kW and 95 GHz MMW [13]. The experimental results indicated that the power density necessary to achieve pain intolearability in 90% of the human subject population is 3.03 W/cm 2 when the duration of the exposures was fixed at 3 seconds. The effects of variable spot size (beam diameter) on human exposure to 95 GHz MMW were investigated in [14] with a small number of volunteers. The experiments confirmed that repel times of stationary subjects decreased with increasing beam size, although the strength of this relationship varied with power density.
In this paper, we study the mathematical framework for human subjects' behavioral responses when exposed to millimeter wave radiation. The electromag-netic energy absorbed in the skin increases the skin temperature via dielectric heating. High temperature activates the heat-sensitive nociceptors which produce a stimulus that is sent to the Central Nervous System (CNS) [15]. When the stimulus is sufficiently strong, the withdrawal reflex is initiated and the subject moves away from the exposure. The withdrawal reflex is also called nociceptive flexion reflex and is a spinal reflex intended to protect the body from noxious stimulations [16]. In the ADT CHEETEH-E model [17] recently developed at the Institute for Defense Analyses, the process from the power emitted at the radiation antenna to the subject's behavioral response is divided into multiple sub-steps, and model components were proposed to describe individual sub-steps [17]. In this study, we follow these model components as general guidelines and use only the physical principles of these model components. We construct a concise dose-response relation for predicting the occurrence of withdrawal reflex based on the spatial temperature profile. The resulting concise model is independent of the specific function forms and parameters used in the sub-steps of ADT CHEETEH-E [17]. In our formulation, the dose quantity is defined as the volume of activated region where the temperature is above the activation temperature of nociceptors. The dose-response model for withdrawal reflex is specified by only two parameters: the activation temperature of nociceptors ( act T ) and the critical threshold on the dose quantity ( c z ). When the spatial temperature profile at reflex is measurable, the two model parameters act T and c z can be determined from test data. The inference formulation for calculating act T and c z is based solely on the measured temperature profiles. It does not depend on any symmetry of temperature profile or electromagnetic heating or heat conduction. It does not require any input parameter.
The dose-response model maps a given spatial temperature profile to the binary outcome of subject's reflex response. We connect it to a time evolution model governing the temperature increase for electromagnetic heating in the skin. The result is a composite model that takes as input the electromagnetic beam deposited on the skin, which is characterized by the beam spot area and the power density. For the temperature evolution, we first consider the simple case of uniform temperature over beam cross-section with no heat conduction (which we shall call idealized case A). The exact solution of temperature is proportional to the time and proportional to the applied power density. Its spatial shape is determined by that of the electromagnetic heating. In order to assess the validity of no-conduction approximation, we consider the case of uniform temperature over the beam cross-section with heat conduction in the depth direction (which we shall call idealized case B). To pinpoint the effects of individual parameters and to facilitate the exact solution, we carry out non-dimensionalization. The temperature solution of the non-dimensional system is the product of the applied power density and a standardized parameter-free function, which we shall call the normalized temperature. The normalized temperature as a function of non-dimensional time and depth has no dependence on the power density or any other parame-wave toward the skin surface of a subject, located at a certain distance away. Let N be the number of radiators in the antenna and i ℘ be the power output of the i-th radiator. The power output of the antenna is described by the vector . We review the formulation components used in [17] for modeling sub-steps in the process from the electromagnetic power output of the antenna eventually to the subject's behavioral response.
where: where y is the coordinate of depth from the skin surface and 1 µ is the characteristic depth that the millimeter wave penetrates into the skin. In the formulation here r is a vector in 3  , restricted to the 2-D skin surface, describing the 2-D coordinates on the skin surface. In a local 3-D coordinate system with the depth direction selected as an axis, r is effectively a vector in 2  , and mathematically ( ) , y r represents the 3-D coordinates in the skin. 5) Temperature as a function of spatial coordinates and time The temperature distribution is governed by the heat conduction along the depth direction with a source term of electromagnetic heating. where ρ is the mass density of the subject's skin; p C is the specific heat capacity of the skin and; -K is the thermal conductivity of the skin.
In the initial boundary value problem above, an insulated boundary condition is imposed at the skin surface ( 0 y = ). 6) Number of heat-sensitive nociceptors activated in a local voxel Each small voxel either has all its heat-sensitive nociceptors activated or has none of them activated depending on the average temperature.
where -voxel j is a volume element in the computational discretization;

A Concise Model for the Withdrawal Reflex
This section focuses on the occurrence of withdrawal reflex, observed as the subject moving out of beam in tests. At a given time t, we consider the event that the spatial temperature profile results in withdrawal reflex. In terms of the subject's behavioral response, this event is simply described by To formulate a concise model, we follow the model components proposed in [17] as guidelines. We adopt the same assumption that the perceived pain level at time t is solely determined by the number of activated nociceptors at time t, which in turn is solely determined by the spatial temperature profile at time t. In other words, given the temperature profile at time t, the past history of the temperature profile does not affect the perceived pain level at time t. Our strategy is to use only the principles of these model components, not the specific function forms. The goal is to construct a formulation i) that is broadly applicable without assuming a particular function form for any internal sub-step and ii) that is concise with only a few model parameters. Specifically, we require the model to have 3 design features: [leftmargin = 1.5 cm].
• it maps a given spatial temperature profile to the corresponding binary outcome with regard to the occurrence of withdrawal reflex; • it is independent of the function forms adopted in [17] and; • it has only two parameters. Building upon the model components of [17], we proceed with the analysis steps below.
• Based on the relation between ( ) Using this expression of ( ) x t , we write the event in terms of ( ) , , T y t r : Result (1) The deterministic dose-response model is

Inferring Model Parameters from Spatial Temperature Profiles
We study the methodology of determining the two parameters ( c z and act T ) in model (4) from test data. We consider the hypothetical situation where the tem- The dose-response model described in (4)  T and c z . In experimental setup, two quantities are tunable: the power density deposited on the skin surface dep P and the beam spot area A. When the applied power density is uniform inside a geometric area and zero outside, the beam spot is naturally defined as that area. When the applied power density is a function of the 2-D coordinates on the skin surface: ( ) dep dep P P = r , the scalar beam spot area A refers to a characteristic area of the beam cross-section. For a Gaussian beam with radius w, the beam spot may refer to the circle of radius w around the beam center, which is the region of ( ) ( ) or equivalently the region of ( ) ( ) . Alternatively, the beam spot area may refer to the peak effective area, which is defined as In test data, information on the unknown parameters act T and c z is contained in the reflex time ref t and in the spatial temperature profile at reflex

Theoretical Behaviors of the System in the Case of No Heat Conduction
In this section, we study the case of no heat conduction. We examine the behaviors of the temperature distribution, the time until reflex, the latency in withdrawal reflex, and the constraint function on model parameters based on the spatial temperature profile at reflex. To facilitate the analysis in a simple theoretical setting, we first introduce idealized cases.

Idealized Case U and Case A
Idealized case U is characterized by the two assumptions below: 1) At any given time, the temperature is uniform over the beam cross-section A (i.e., independent of r ) and is a decreasing function of depth y.
2) Outside the beam cross-section, the temperature is always below the nociceptor activation temperature (for example, at the normal body temperature).
One particular situation of case U is when the initial temperature of skin is a constant below the activation temperature and the applied power density is uniform over the beam cross-section A and zero outside. Mathematically, idealized case U is characterized by , , , , decreasing with , for Case U , , , for

T y t T y t y A T y t T A
In case U, it is mathematically more convenient to view the constraint func- Case U : , Equation (8) Case A satisfies the conditions of case U, and is a special situation of case U.

Selecting Test Conditions for
It follows from (8) Note that (12) is the constraint function for case A (the case of no heat conduction) and it is independent of the power density dep P , which specifies how fast the electromagnetic beam heats up the skin. When the effect of heat conduction along the depth direction is included, dep P will affect the steepness of spatial temperature profile at reflex. A smaller power density dep P heats up the skin slower and the heat conduction has more time before reflex to smooth out the temperature profile. Figure 1 compares constraint functions constructed from simulated test data In simulations, we used ( ) The choice of ( ) does not actually alter the essential shape of constraint functions. If we divide (12) by ( ) and use the scaled temperature as the dependent variable, then (12) is completely specified by A µ , independent of ( ) In summary, in case A, to determine model parameters act T and c z , we need to carry out tests with different beam spot areas: one with a moderate value of A µ to make the constraint function act T vs. c z slant, the other with a relatively large value of A µ to make act T vs. c z flat. Notice that the meanings of phrases "moderate value of A µ " and "relatively large value of A µ " are not precise since quantity A µ is not dimensionless. This issue will be addressed when we discuss non-dimensionalization.

Time until Withdrawal Reflex vs. Beam Spot Area
We study how the reflex time ref t varies with the beam spot area A while the power density dep P is maintained at a constant level. In case A (the case of no heat conduction), the analytical expression of ref t is given in (10). In expression (10), when the beam spot area A is fixed, the reflex time ref t is inversely proportional to the applied power density dep P . When dep P is fixed, as beam spot Eventually ref t reaches a constant level above zero for large A. This can be seen in the Taylor expansion of (10) as ( ) In all realistic situations, heat conduction is always present. Case A (the case of no heat conduction) corresponds to the situation where the effect of heat conduction is relatively small in comparison with others. In the non-dimensional analysis of the next section, we will see that the effect of heat conduction is negligible only in a region away from the skin surface and only over a short time.
Even during a short time, near skin surface, the heat conduction may still be significant or even dominant in the temperature evolution. Expression (10)  A with the effect of heat conduction will be carried out after the non-dimensionalization analysis. We will see that with the effect of heat conduction, the (1/A) term in expansion (13) is invalid.

Latency in Withdrawal Reflex
We hypothesize that when the number of nociceptors activated reaches a threshold, the stimulus sent to the Central Nervous System (CNS) is strong enough to initiate the withdrawal reflex. However, it takes time for the signal to travel to the CNS, for the CNS to send a signal to muscles, and for muscles to act upon the signal to carry out the withdrawal reflex before the actual withdrawal reflex action (i.e., the subject moving out of the beam) is observed in tests. To facilitate the discussion of biological latency in the observed withdrawal reflex, we first introduce proper mathematical notations for these time instances. Let the time until a stimulus strong enough for initiating withdrawal reflex is generated; we regard ref t as the true underlying reflex time in the sense that even if the beam power is turned off at that point, withdrawal reflex will still occur; We assume that the time delay del t is an intrinsic property of the test subject, not affected by parameters of the electromagnetic beam, such as power density ( dep P ) or beam spot area (A). In case A (the case of no heat conduction), the true underlying reflex time ref t is given by (10). The observed reflex time obs t has the expression: The plot of obs t vs

Non-Dimensionalization and Solution of Heat Conduction in the Depth Direction
We study the effect of heat conduction along the depth direction of the skin (the y-direction of the coordinate system). We first introduce case B which has all assumptions of case A except that it includes the effect of heat conduction in the depth direction with uniform thermal conductivity: uniform power density over beam spot : for and 0 for Case B heat conduction in with uniform conductivity uniform initial temperature : , 0 In case B, the temperature distribution is governed by , , 0, , , 0

Non-Dimensionalization
We introduce a temperature scale T ∆ , the characteristic difference between the initial temperature 0 T and the nociceptor activation temperature act T . Note that T ∆ serves only as the temperature scale for non-dimensionalization. T ∆ does not need to be the exact difference between 0 T and act T . We look at the physical dimensions of various quantities: We introduce length scale and time scale Here s y is the depth at which the beam decays by a factor of 1 e − , and s t is the time at which a δ function spreads to a Gaussian distribution of standard deviation s y . We use these characteristic scales to carry out non-dimensionalization. Non-dimensional depth and time: Non-dimensional temperature as a function of ( ) nd nd , y t : The non-dimensional temperature distribution offers two mathematical ad- Notice that the normalized Equation (20) is parameter-free.

Analytical Solution
To solve problem (20), we view the forcing term as the time integral of impulse forcing. Let ( ) , G y t denote the solution of a unit impulse forcing at 0 t = , which is governed by In Appendix B, we derive that ( ) , G y t has the expression ( ) The pre-scaling physical temperature distribution ( ) Equation (23)  which is parameter-free.

Effect of Heat Conduction
We investigate the consequence of neglecting heat conduction (i.e., completely turning off the heat conduction in Equation (20)). Specifically, we remove the which yields the no-conduction approximation of ( ) It is important to clarify the precise difference between Equations (20) and (24). They do not correspond to different regimes of heat conductivity. Rather, they are both obtained in non-dimensionalization using exactly the same parameters. The only difference is that at the end, in Equation (24) we discard the terms involving heat conduction.
To assess the validity of neglecting heat conduction, we compare Figure 2 at several time instances. Here y and t are the nondimensional depth and time. When t is small and y is away from 0, provides a good approximation to  For the physical temperature distribution, the no-conduction approximation is  Approximation (27) is valid when the non-dimensional time

Theoretical Behaviors of the System in the Case with Heat Conduction
In this section, we study the behaviors of case B (the case with heat conduction), based on its analytical solution. We examine the reflex time vs. beam spot area, the latency in withdrawal reflex, asymptotics of the normalized temperature for small t and for large t, the energy consumption at small and at large applied power density, the spatial temperature profile at reflex, and constraint functions constructed from spatial temperature profiles. In case U defined in (7), which includes case B, the governing equation for ref t is given in (8). Using the exact solution (23) for case B, we write the equation in terms of the normalized temperature distribution

Reflex Time vs. Beam Spot Area
Substituting expansion form (29) into Equation (28) and carrying out the Taylor expansion of ( ) Based on governing Equation (20) and expression (22) of ( ) , H y t , we can derive: Matching corresponding terms of where 0 t and A c are solved from (33) and are given by

Latency in Withdrawal Reflex
We study the behavior of biological latency (time delay) del t in observed reflex in case B (the case with heat conduction). We investigate the mathematical for- In the limit of c A z µ → ∞ , the observed reflex time obs t as a function of dep P has the form: ( ) ( ) ( ) In (38) and (39), the number of unknowns is reduced to two: 1 c and 2 c .
Next, we get rid of unknown 1 c . For that purpose, we consider quantity obs v defined as ( ) From (38) and (39), it follows that obs v is a function of 2 0 c p where function ( ) R u is defined in terms of ( ) (31) and is parameter-free. As a result, function ( ) R u introduced above is also parameter-free. Figure 4 demonstrates that ( ) R u is an increasing function. This observation indicates that the inverse function ( ) 1 R − ⋅ is well defined.
In Appendix C, we derive asymptotics of ( ) R u u u u The last line of (44) gives the predicted function of observed reflex time vs.
applied power density, based on the 3 data points. Algorithm (44) is parameter-free. It calculates the time delay del t using only the 3 data points. Even though algorithm (44) is derived based on the theoretical behavior of case B in the limit of large beam spot area, operationally it can be applied in any situation since it does not require any input parameter. In the next section, we will test it in the case of a Gaussian beam where the temperature is not uniform over the beam cross-section and the beam radius is finite.

Temperature at a Fixed Depth for Small Time and Large Time
In this subsection, we study the time evolution of the normalized non-dimensional temperature

( )
, H y t at a fixed y, respectively for small t and for large t.
Here the pair ( )   The linear portion of the temperature growth at depth y is attributed to the heating at y from the electromagnetic wave penetrating into the skin. The portion above the linear growth is caused by the positive net heat gain via conduction. For a small interval around y, the net heat gain via conduction is positive when the heat in-flow from the upstream is more than the heat out-flow to the downstream. The temperature growth caused by the electromagnetic heating is augmented by the conduction when the net heat gain via conduction is positive.
Mathematically, at depth y, the net heat gain via conduction is positive when ( ) 2 2 , 0 H y t y ∂ ∂ > , which is true for small t. We will show below that for large t, however, the net heat gain via conduction is negative. As a result, the temperature increase caused by electromagnetic heating is diminished by heat conduction, and the realized temperature increase for large t is much smaller than the linear growth.
We use the expansion of ( ) erfcx z given in (68) of Appendix C.
( ) Integrating the expansion of ( ) , G y t with respect to t, we obtain ( ) ( )    We look into the net heat gain/loss due to conduction, for small t and for large t. Differentiating (45) and (47) The results in (48) are valid only when (45) and (47) are valid, respectively, in the regimes of small t and large t. This is particularly evident in line 2 of (48), which suggests that at fixed depth y, the regime of large t has to satisfy 2 e y t > π.
Since there is no heat in-flow at depth 0 (skin surface), the conduction always leads to an overall heat loss for the interval [ ] 0, y . The local net heat loss is initially confined to near 0 and gradually propagates to the entire interval. For large y, it takes long time for the local net heat loss to reach depth y. It takes even longer time for the net heat loss to offset the previously accumulated heat gain at y and then to substantially neutralize the electromagnetic heating.

Reflex Times at Large and at Small Applied Power Densities
In case B, the temperature is uniform over the beam cross section and the activated region is a cylinder. Given the beam spot area A and the threshold c z for activated volume, withdrawal reflex occurs when the temperature at depth c z A reaches the activation temperature act T . In this subsection, we use the asymptotics of temperature for small time and for large time (obtained in subsection 7.3) to study the reflex time at large and at small applied power densities.
We also compare the results with the no-conduction approximation of the reflex time given in (10), which is inversely proportional to the applied power density. In case B, the reflex time ref , as described in (8). The non-dimensional version of the equation is given in (28). For our discussion here, we write the left side of (28) in terms of act,nd T and dep,nd P . The relative reduction in reflex time actually increases slightly when the applied power density is slightly decreased as long as it is still in the regime of large power density so that the small time approximation (45)   In Figure 6, we make an important observation that the energy consumption has a minimum with respect to applied power density. We now demonstrate the minimum analytically. Based on asymptotics (51)  large. In general, the range of moderately large dep,nd P is also a good compromise between inducing withdrawal reflex quickly and preventing burn injury of accidental over-heating. For large dep,nd P , burn injury may occur even when the power is turned off shortly after the internal initiation of withdrawal reflex.

Spatial Temperature Profile at Reflex
We examine the temperature profile at reflex and how it varies as the applied power density increases from small to large. We compare it with the no-conduction approximation, which is independent of the applied power density.
We follow the equation and notation for the reflex time used in subsection for several values of dep,nd P and its no-conduction approximation. As dep,nd P increases, the reflex time decreases to zero, and the heat conduction has little time to take its effect. Figure 7 demonstrates that at a fixed depth 0 y > , as dep,nd P increases the temperature at reflex converges to that in the case of no heat conduction. Therefore, at any fixed depth away from the skin surface, the effect of heat conduction is eventually negligible when the applied power density is large enough.

Selecting Test Conditions for Determining
The situation of inferring parameters ( ) act , c T z from temperature profiles at ref t is similar to that in the case of no heat conduction. Once the latency del t is determined, the true reflex time is calculated from the observed reflex time: τ τ (57) is a constraint function in the form of The solution τ is a function of ( )

The Case of a Gaussian Beam with Heat Conduction
We consider a Gaussian beam with power density where ( )  using (60). Following this procedure, we generate simulated data and then we use the data to test the theoretical predictions derived in the previous section for case B.
To facilitate the analysis, we introduce parameters 2 c and 1 c as in case B, and write temperature distribution (60) in a non-dimensional form, in terms of 2 c and 1 c .    Figure 10 as exact since it is numerically accurate to the computer precision.

Estimating the Latency
We study the methodology for determining the time delay in the observed with- Given 3 data points, we applied (37) to construct a constraint at each data point. Then we used the joint system of 3 constraints to derive algorithm (44) for determining del t . Algorithm (44) is based solely on the 3 data points. It does not require any input parameter. In the idealized situation above, del t solved from algorithm (44) is exact. In this subsection, we test the performance of algorithm (44) for inferring del t in the case of a Gaussian beam with a finite beam radius.
The purpose is to demonstrate the versatility and accuracy of algorithm (44) when the idealized conditions above are not met. Here we focus on the accuracy of estimating the time delay since del t is our objective while coefficients 2 c and 1 c are just auxiliary unknowns accompanying del t in algorithm (44).
The simulated curve of obs t vs ( ) 0 dep P for a Gausian beam is shown in Figure  10.  Figure 10. These 3 data points and only these 3 data points are then used in algorithm (44) to estimate del t . A Gaussian beam with finite beam radius does not satisfy the idealized conditions above. Nevertheless, we apply algorithm (44) to the 3 data points in Figure 10 to estimate ( ) . Operationally, algorithm (44) is universally applicable in any situation since it does not require any input parameter.
Once we obtain the estimated values

Scaling Properties of the Formulation for Estimating Latency
In addition to its high accuracy for estimating the time delay in situations where the idealized conditions are not met, formulation (44) and the predicted function based on it, have several important properties. We now discuss these properties. 3) The prediction is invariant with respect to a scaling of ( )  Therefore, for the purpose of determining del t , we only need to measure There is no need to measure ( )     In case B, we showed analytically (in subsection 7.4) that the energy consumption attains a minimum at an intermediate range of applied power density. The minimum energy consumption is attributed to the two opposite effects of heat conduction on temperature increase, respectively over short time and over long time. In case B, with a given beam spot area, withdrawal reflex occurs when the temperature at a fixed depth reaches the nociceptor activation temperature. Mathematically, we derived (in subsection 7.3) that at a fixed depth, over short time the net heat gain via conduction is positive. When the applied power is large, the reflex occurs in short time. Over that short time, the heat conduction augments the temperature increase of electromagnetic heating, speeds up the process of reaching the activation temperature, and reduces the energy consumption. In the regime of large applied power, increasing the applied power further makes the reflex time very short and leaves very little time for the conduction to take its effect in reducing the energy consumption. As a result, in the regime of large applied power, the energy consumption increases with the applied power. At a given depth, the positive net heat gain via conduction depends on the heat in-flow from upstream. At the skin surface, however, there is no heat in-flow. As time goes on, the effect of insulated boundary propagates in the depth direction in the form of attenuating the heat flow. Mathematically, at a fixed depth, eventually the net heat gain via conduction turns negative and the magnitude of net heat loss grows with time. When the applied power is small, it takes long time to induce the reflex. Over that long time, the heat conduction yields net heat loss, diminishes the temperature increase of electromagnetic heating, slows down the process of reaching the activation temperature, and drives up the energy consumption. In the regime of small applied power, reducing the applied power will give the conduction more time to take its effect in neutralizing the electromagnetic heating, and increase the energy consumption. Consequently, In the regime of small applied power, the energy consumption increases when the applied power is lowered. The transition between these two regimes produces a minimum for the energy consumption. Although the behaviors of these two regimes were derived for the idealized case B, in the case of a Gaussian beam, we observe the same behaviors in Figure 11: in the regime of small ( )   . For each data set, we apply formulation

Concluding Remarks
We studied theoretically the occurrence of heat-induced withdrawal reflex caused by exposure to an electromagnetic beam. We investigated several aspects of the problem, including i) non-dimensionalization to pinpoint the effects of parameters; ii) normalization of the non-dimensional temperature into a parameter-free function; iii) forward prediction of system behaviors given model American Journal of Operations Research Parameters; iv) asymptotic behaviors in certain regimes of parameters; and v) backward inference of model parameters based on measured data.
First, we reviewed the model components used in ADT CHEETEH-E [17], and we synthesized them to distill a concise dose-response model for predicting the occurrence of withdrawal reflex based on the given spatial temperature profile. In the dose-response relation, the dose quantity is defined as the volume of the activated region of heat-sensitive nociceptors. Under the assumption that the nociceptor density is uniform in the skin, the activated volume is proportional to the total number of nociceptors activated, and thus, increases monotonically with the perceived pain level, which controls the occurrence of withdrawal reflex. In a deterministic setting, withdrawal reflex occurs precisely when the activated volume reaches a critical threshold. The resulting dose-response relation described in (4) is specified by two parameters: the nociceptor activation temperature act T and the critical threshold c z on the activated volume for withdrawal reflex.
The concise model has the advantage that the two parameters can be determined from test data in the situation where the reflex time and the spatial temperature profile at reflex are measurable. Parameters act T and c z are constrained by function (5), which is constructed solely from the measured temperature profile at reflex. The inference method based on (5) is parameter-free, requiring no input parameter. The data set measured at each test condition provides only one constraint on ( ) act , c T z . To determine the two unknown parameters, several tests at various conditions are needed to produce distinct constraint functions. The intersection of distinct constraint curves yields the true value ( ) act , c T z .
The dose-response relation predicts the occurrence of withdrawal reflex from the given spatial temperature profile. We connected it with a time evolution model for the temperature increase of electromagnetic heating. The result is a composite model that takes as input the test condition described by two tunable quantities: the beam spot area and the applied power density. Other quantities, such as the initial temperature, the characteristic depth of the electromagnetic wave penetrating into the skin and the nociceptor activation temperature, are viewed as parameters. In the composite model, the time evolution model produces the temperature distribution from the test condition, and then the doseresponse relation uses the temperature distribution to determine the time of withdrawal reflex. We solved the composite model analytically in two idealized cases, first in the case of no heat conduction and uniform electromagnetic heating over beam cross-section (case A).
To examine the validity of the no-conduction assumption, we investigated the effect of heat conduction on temperature distribution. In the case of uniform electromagnetic heating over beam cross-section with heat conduction in the depth direction (case B), we carried out non-dimensionalization to facilitate the exact solution and to pinpoint the effects of model parameters. We solved the non-dimensional system analytically to obtain the exact solution given in (22).
As a function of non-dimensional time and depth, the non-dimensional temperature is the product of the non-dimensional power density and a parameter-free function, which is called the normalized temperature. The normalized temperature has no dependence on the power density or any other parameters. The effects of physical parameters are contained in the non-dimensional variables and quantities. Scaling temperature distribution into a parameter-free function significantly simplifies the structure of exact solution and reveals its dependence on parameters. We compared the normalized temperatures in the presence and in the absence of heat conduction. We found that the no-conduction approximation is valid only in the region away from the skin surface and only over a short time. Both conditions need to be satisfied in order to make the effect of heat conduction negligible. This result indicates that the no-conduction approximation is not an adequate tool for analyzing the temperature evolution near the skin surface.
Using the exact solution obtained in case B (the case with heat conduction), we studied theoretical behaviors of the system. Below are the findings.
• As the beam spot area A increases, the reflex time decreases rapidly to a positive value above zero. The convergence is described by 1/A 2 in (34).
• We define the true reflex time as the moment when the number of activated nociceptors is sufficiently large to produce a perceived pain level exceeding the tolerance and thus to initiate the withdrawal reflex. The observed reflex time is when the reflex response (i.e., the subject moving out of the beam) actually takes place. There is a latency (time delay) between the observed reflex time and the true reflex time. We assume that the latency is intrinsic to the subject being tested and is independent of the applied power density. We derived algorithm (44) for inferring the time delay from measured values of observed reflex times at a sequence of applied power density values. Once the time delay is obtained, the true reflex time is calculated from the observed reflex time.
• The true reflex time decreases when the applied power density is increased.
To examine the behaviors respectively at large and at small applied power density, we carried out asymptotic expansions of the normalized temperature respectively for small and large time. From the expansions of the normalized temperature, we derived asymptotic approximations of the reflex time, respectively, for large applied power density given in (51) and for small power density given in (52). At a fixed beam spot area, the product of beam power density and reflex time is proportional to the energy consumed for inducing withdrawal reflex. The asymptotic analysis indicates that when we increase the power density gradually in the small regime, initially the energy consumption decreases and then it attains a minimum; eventually in the large regime, further increase of power density leads to an increase in energy consumption. The two sides of the energy consumption curve surrounding the minimum correspond to the two opposite effects of heat conduction, respectively over short time and over long time. At a fixed depth, over short time, heat conduction augments the temperature increase of electromagnetic heating. It leads to a temperature growth faster than the no-conduction solution and reduces the energy consumption by shortening the reflex time. In contrast, over long time, heat conduction diminishes the temperature increase of electromagnetic heating. It results in a much slower temperature growth than the no-conduction solution and drives up the energy consumption needed for inducing withdrawal reflex.
• We examined the spatial temperature profile at reflex and how it varies with the applied power density. We found that at a fixed depth, for large power density, the temperature at reflex converges to the no-conduction solution.
This result indicates that it is appropriate to neglect the effect of heat conduction only when the location is away from the skin surface and when the (non-dimensional) power density is large.
• In the situation where the spatial temperature profile at reflex is measured in tests, we studied the methodology for determining ( ) act , c T z from test data.
( ) act , c T z is constrained by Equation (5), which is parameter-free and is constructed based solely on the measured spatial temperature profile. The intersection of two or more distinct constraint curves, calculated from test data at various test conditions, determines the true value ( ) act , c T z . We explored how to select test conditions to produce significantly distinct constraints on ( ) act , c T z . We found that the steepness of ( ) act c T z is negatively associated with the beam spot area. For the purpose of estimating ( ) act , c T z , we need to carry out tests with moderate beam spot area and with large beam spot area. This is to ensure that the constraint curves constructed from test data are substantially distinct and have a well-defined intersection.
We conducted a case study of a Gaussian beam with finite beam radius and with heat conduction in the depth direction. Since it has no closed-form solution, the case study was carried out using numerical simulations. The goal was to test the theoretical results predicted for the idealized case of uniform temperature over beam cross-section (case B). We first examined the performance of algorithm (44) for determining the unknown latency from measured values of observed reflex times at a sequence of applied power density values. Algorithm (44) was derived under the assumptions of case B and in the limit of beam spot area approaching infinity. (44) calculates the time delay solely from 3 data points of ( ) dep obs , P t where dep P is the applied power density and obs t is the corresponding observed reflex time. Algorithm (44) does not require any input parameter. Operationally, it is universally applicable in all situations where the observed reflex time is measured. We applied it to the case of a Gaussian beam with a finite beam radius. In Figure 10, the estimated time delay is very close to the true value, suggesting the wide applicability of algorithm (44) beyond the idealized case B.
In addition to its high accuracy, formulation (44) has an important scaling property: the prediction by (44)  In the case of a Gaussian beam, we examined numerically the energy consumption for inducing withdrawal reflex as a function of applied power density.
We first considered the situation where the power is turned off at the true reflex time when the activated volume reaches the threshold to start the internal initiation of withdrawal reflex. The observed reflex action (the subject moving out of beam) occurs with time delay del t after the internal initiation. We observed that a minimum energy consumption is attained at a moderately large beam center power density, similar to what we theoretically predicted in the idealized case B.
This observation suggests that the existence of minimum energy consumption may be a general phenomenon in all situations. In the situation where we keep the power on for a short period beyond the true reflex time, the energy consumption has a more pronounced minimum, which is attained at a more moderate value of the applied power density.
In the case of a Gaussian beam, we tested the strategy of selecting test conditions for determining parameters act decreases when the beam size is increased. To obtain substantially distinct constraint curves for ( ) act , c T z , we need to carry out tests with a moderate beam radius and with a large beam radius. In summary, in a deterministic setting, we studied mathematically the system from the electromagnetic power deposited on the subject's skin to the subject's behavioral response. Our theoretical findings provide insight into how the system behaves in various regimes of model parameters and how the system behavior varies in response to changes in parameters. Furthermore, through theoretical formulation and analysis in idealized cases, we established a mathematical framework for designing tunable parameters in tests to fully sample the effects of hidden parameters in test data. Properly sampling the effects of model parameters in test data is a vital step toward reliably determining the values of these parameters.
As t goes to infinity, ( ) h t increases unbounded. Thus, the range of ( ) h t is [ ) 0, +∞ .
Step 2: Using the results of ( ) h t obtained in Step 1, we see that