^{1}

^{2}

^{3}

^{*}

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.

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 [

Considerable work has been done in MMW interactions with the human body [^{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 [^{2} at 94 GHz.

Heating of tissues by electromagnetic waves has also been studied in [

In [

In [^{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 [

According to a review published in 2016 [^{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 [

In this paper, we study the mathematical framework for human subjects’ behavioral responses when exposed to millimeter wave radiation. The electromagnetic 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) [

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 parameters. The effects of physical parameters are contained in the non-dimensional variables and in the non-dimensional power density (the multiplier coefficient). Comparing the solutions, respectively, in the presence and in the absence of heat conduction, we find that the no-conduction approximation is appropriate only in the region away from the skin surface and only over a short time. With the analytical solution, we explore the theoretical behaviors of the system and discuss the issues listed below.

· We derive the asymptotic expansion of reflex time at large beam spot area, and compare the result with that in the case of no heat conduction.

· We investigate the biological latency (time delay) in the observed withdrawal reflex. We formulate an algorithm for determining the time delay in the observed withdrawal reflex. The algorithm is based solely on three data points of observed reflex time vs. applied power density. The algorithm is parameter-free, and thus, is operationally applicable in all situations.

· We carry out asymptotic analysis on the normalized temperature respectively for small time and for large time. Building on the asymptotic behaviors of normalized temperature, we construct asymptotic approximations of the reflex time respectively for large and for small applied power density.

· Using the asymptotic results obtained, we examine the energy consumption by the time of reflex, as a function of applied power density. We find that the energy consumption attains a minimum at a moderately large value of applied power density.

· We examine the spatial temperature profile at reflex and demonstrate that it converges to the no-conduction solution at large power density.

· We study how to select test conditions for determining model parameters ( T act , z c ) from the measured temperature profiles at reflex. Each data set yields only one constraint on ( T act , z c ) . To determine both T act and z c simultaneously, we need to obtain distinct constraints on ( T act , z c ) , which is achieved by carrying out tests both with large beam spot area and with moderate beam spot area.

Finally, we do a case study of Gaussian beams using numerical simulations. We revisit the theoretically predicted behaviors of the system derived in the idealized case B (uniform temperature over beam cross-section). For Gaussian beams (which do not satisfy the assumptions of case B), 1) we evaluate the performance of the algorithm for estimating the biological latency from a sequence of observed reflex times; 2) we examine the existence of energy consumption minimum with respect to the applied power density; and 3) we test the strategy of selecting optimal test conditions for determining model parameters T act and z c .

The ADT CHEETEH-E model was proposed in [

1) Electric field near the skin of the subject

E ( r ) ︸ EElectric field on skin at r = ∑ i = 0 N − 1 β i 2 η 0 ℘ i 1 4 π R i e ( − α 2 R i − j k R i ) f ( R ^ i )

where:

- r is the position vector of a point on the skin surface of the subject;

- r i is the position vector of the i-th radiator of the antenna;

- R i = | r − r i | and R ^ i = r − r i | r − r i | are respectively the magnitude and the unit direction of vector ( r − r i ) ;

- α is the atmospheric attenuation coefficient;

- k = 2 π / λ is the wave number of the electromagnetic wave; λ is the wave length; and j = − 1 .

2) Power incident per area on the skin

P incident ( r ) ︸ Power incident per area on skin at r = | E ( r ) | 2 2 η 0

3) Power deposited per area on the skin

P dep ( r ) ︸ Power deposited per area on skin at r = ( 1 − γ ) P incident (r)

4) Power absorbed per volume into the skin

q ( r , y ) ︸ Power absorbed per volume into skin at depth y = P dep ( r ) μ e − μ y

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 ( r , y ) 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.

{ ρ C p ∂ T ( r , y , t ) ∂ t = ∂ ∂ y ( K ∂ T ( r , y , t ) ∂ y ) + q ( r , y ) ∂ T ( r , y , t ) ∂ y | y = 0 = 0, T ( r , y ,0 ) = T 0 (y)

where

- ρ is the mass density of the subject’s skin;

- C p 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 ( y = 0 ).

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.

x j ( t ) ︸ NNumber of nociceptors activated in voxel j = { V j ρ noc , if T j ¯ ( t ) ≥ T act 0, if T j ¯ ( t ) < T act

where

- voxel j is a volume element in the computational discretization;

- V j is the volume of voxel j;

- ρ noc is the density of heat-sensitive nociceptors in the subject’s skin;

- T j ¯ ( t ) is the average temperature of voxel j at time t and;

- T act is the activation temperature of heat-sensitive nociceptors.

7) Total number of heat-sensitive nociceptors activated at time t

x ( t ) ︸ Total number of nociceptors activated = ∑ j = 0 M − 1 x j (t)

where M is the number of voxels in the computational discretization.

8) Perceived pain level

h Dol ( t ) ︸ Perceived pain level in Dol scale at time t = a 1 + e x p ( − ( l n x ( t ) − b ) c )

9) Motivation-modulated perceived pain level

h ^ Dol ( t ) ︸ Motivation modulated perceived pain level in Dol scale at time t = h Dol ( t ) − m m 0

10) Subject’s behavioral response

g ( t ) ︸ Subjec ′ sbehavioral response at time t = { 0 ( no movement ) , if h ^ Dol ( t ) < Y Low 1 ( flinch but remain in beam ) , if Y Low ≤ h ^ Dol ( t ) < Y High 2 ( move out of beam ) , if h ^ Dol ( t ) ≥ Y High

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 Event { g ( t ) = 2 } . To formulate a concise model, we follow the model components proposed in [

· 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 [

· it has only two parameters.

Building upon the model components of [

· Based on the relation between g ( t ) and h ^ Dol ( t ) (the motivation-modulated perceived pain level) described in model component 10 above, we write.

Event { g ( t ) = 2 } = Event { h ^ Dol ( t ) ≥ Y High }

· h ^ Dol ( t ) is a monotonically increasing function of h Dol ( t ) (the perceived pain level). In model component 9 above, this function is set to: h ^ Dol ( t ) = h Dol ( t ) − m / m 0 . Here we write generally h ^ Dol ( t ) = F 1 ( h Dol ( t ) ) and require only that F 1 ( ⋅ ) be an increasing function. Applying F 1 − 1 ( ⋅ ) to h ^ Dol ( t ) ≥ Y High , we write:

Event { g ( t ) = 2 } = Event { h Dol ( t ) ≥ F 1 − 1 ( Y High ) }

· h Dol ( t ) increases monotonically with x ( t ) (the total number of nociceptors activated). In model component 8 above, h Dol ( t ) is set to:

h Dol ( t ) = a 1 + e x p ( − ( l n ( x ( t ) ) − b ) c ) . Again, we write generally h Dol ( t ) = F 2 ( x ( t ) ) and require only that F 2 ( ⋅ ) be an increasing function. Applying F 2 − 1 ( ⋅ ) to h Dol ( t ) ≥ F 1 − 1 ( Y High ) , we write:

Event { g ( t ) = 2 } = Event { x ( t ) ≥ F 2 − 1 ( F 1 − 1 ( Y High ) ) }

· At a given time t, x ( t ) is determined from the spatial temperature profile T ( r , y , t ) :

x ( t ) = ρ noc ∫ I ( T ≥ T act ) ( r , y , t ) d r d y ︸ Activatedvolume

where the indicator function I ( T ≥ T act ) is defined as:

I ( T ≥ T act ) ( r , y , t ) ≡ { 1, if T ( r , y , t ) ≥ T act 0, if T ( r , y , t ) < T act

Using this expression of x ( t ) , we write the event in terms of T ( r , y , t ) :

Event { g ( t ) = 2 } = Event { ∫ I ( T ≥ T act ) d r d y ≥ 1 ρ noc F 2 − 1 ( F 1 − 1 ( Y High ) ) } (1)

Result (1) leads to a deterministic dose-response relation for withdrawal reflex. Given spatial temperature profile T ( r , y ) , we select the activated volume as the single metric predictor variable (the dose quantity) for predicting withdrawal reflex.

z ≡ z ( { T } , T act ) ︸ Dose ≡ ∫ I ( T ≥ T act ) d r d y ︸ Activatedvolume (2)

Here we use { T } as a concise notation for spatial temperature profile T ( r , y ) . Equation (1) tells us that the critical threshold on dose z for withdrawal reflex is

z c ≡ 1 ρ noc F 2 − 1 ( F 1 − 1 ( Y High ) ) (3)

The deterministic dose-response model is

Outcome ( z ) = { 1 ( withdrawal reflex ) , if z ≥ z c 0 ( no withdrawal reflex ) , if z < z c (4)

Model (4) is completely specified by 2 parameters: T act and z c .

· Activation temperature T act is used in calculating dose z in (2).

· Threshold z c is used in determining the binary outcome corresponding to dose z.

Note that threshold z c varies with many internal parameters of ADT CHEETEH-E [

· Y High in the subject behavioral response function g ( t ) ;

· m and m 0 in the motivation modulated perceived pain level h ^ Dol ( t ) ;

· a, b and c in the perceived pain level h Dol ( t ) and;

· ρ noc in the total number of nociceptors activated x ( t ) .

In addition, threshold z c depends on function forms of F 1 ( ⋅ ) and F 2 ( ⋅ ) . The advantage of model (4) is that the effect of all these model parameters and function forms is captured in a single parameter z c . Once the values of z c and T act are known, dose-response model (4) is completely specified.

We study the methodology of determining the two parameters ( z c and T act ) in model (4) from test data. We consider the hypothetical situation where the temperature of the skin as a function of the 3-D coordinates and the time is measurable in experiments.

The dose-response model described in (4) is based on the assumption that the occurrence of withdrawal reflex at time t ref is solely attributed to the spatial temperature profile at t ref . As a result, in the inference method, only T ( r , y , t ref ) is relevant for determining T act and z c . In experimental setup, two quantities are tunable: the power density deposited on the skin surface P dep 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: P dep = P dep ( 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

P dep ( r ) ≥ 1 e 2 P dep ( 0 ) or equivalently the region of

In test data, information on the unknown parameters

Function

Let us summarize and clarify what information regarding model parameters

· At a given test condition

· As described in (5), constraint function

· In the subsequent sections, we will analyze the behaviors of

· If test data from different test conditions provide two distinct constraint equations, then parameters

In the next section, we will explore theoretically what test conditions are likely to provide substantially distinct constraint equations for

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 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

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

In case U, it is mathematically more convenient to view the constraint function in the form of

Equation (8) serves various purposes depending on which are known/unknown. On one hand, given the reflex time

It should be mentioned that case U does not exclude heat conduction in the depth direction. We now add the assumption of no-heat-conduction and introduce case A.

Case A satisfies the conditions of case U, and is a special situation of case U.

Case A is solved analytically in Appendix A. The temperature distribution and the reflex time are given by (see Equations (65) and (66) of Appendix A).

It follows from (8) that the constraint function

Let

We like to get rid of

we view (11) as a constraint function on

We cast the constraint function into the form of

Constraint function (12) is specified by two parameters:

in case A, with two different beam spot areas:

In summary, in case A, to determine model parameters

We study how the reflex time

Case A: no heat conduction

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) for

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

·

·

·

We assume that the time delay

The plot of

Case A: no heat conduction

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:

In case B, the temperature distribution is governed by

We introduce a temperature scale

We introduce length scale and time scale

Here

Non-dimensional depth and time:

Non-dimensional temperature as a function of

Non-dimensional reflex time:

Non-dimensional power density deposited on skin surface:

Non-dimensional beam spot area:

Non-dimensional activation temperature:

The governing equation for

The non-dimensional temperature distribution offers two mathematical advantages.

· For

·

We define the normalized non-dimensional temperature

Here we have denoted

Notice that the normalized Equation (20) is parameter-free.

To solve problem (20), we view the forcing term as the time integral of impulse forcing. Let

In Appendix B, we derive that

where

We integrate

Function

(i) At each fixed y,

(ii) At each fixed t,

(iii) At a fixed y, for small t,

The third property of

In summary, in case B (the case with heat conduction),

The pre-scaling physical temperature distribution

Equation (23) expresses the physical temperature distribution

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

We examine the relative error in approximating

For the physical temperature distribution, the no-conduction approximation is

Approximation (27) is valid when the non-dimensional time

small and the non-dimensional depth

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 subsection 5.3, we studied the asymptotic behavior of reflex time for large beam spot in case A (the case of no heat conduction). The asymptotic result of

In case U defined in (7), which includes case B, the governing equation for

which leads to a non-dimensional equation for

where

the non-dimensional reflex time as defined in subsection 6.1. As the 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

Using these results in expansion (30), we obtain

Matching corresponding terms of

Using expansion (29) and

Case B: with heat conduction

where

Here ^{2} in (34) (the case with heat conduction) is faster than the convergence of 1/A in (13) (the case of no heat conduction).

We study the behavior of biological latency (time delay)

The observed reflect time is the sum of the true reflex time and the unknown time delay:

time

In the limit of

In (37), besides the unknown time delay

· measure three data points of

· then solve for

Specifically, we carry out tests at 3 values of

We examine the difference in the observed reflex time. The difference is independent of

In (38) and (39), the number of unknowns is reduced to two:

From (38) and (39), it follows that

where function

Function

In Appendix C, we derive asymptotics of

These asymptotic results establish analytically the invertibility of function

We now describe the method of determining

Case B: with heat conduction

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

In this subsection, we study the time evolution of the normalized non-dimensional temperature

In solution (22),

The abbreviation

In the no-conduction approximation given in (25),

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

To expand

We use the expansion of

Integrating the expansion of

The constant term

and (47) are based on

are invalid. In the regime of fixed t and large y, the net heat gain at y via conduction is positive. Expansions (46) and (47), however, tell us the behavior in a different regime: at a fixed

We compare the normalized temperature

We look into the net heat gain/loss due to conduction, for small t and for large t. Differentiating (45) and (47) twice with respect to y, we have

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

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

In case B, the reflex time

where

The no-conduction version of Equation (49) is obtained by replacing

Case A: no heat conduction

For large

which yields

Result (51) indicates that heat conduction reduces the reflex time when the applied power density is large. This corresponds to the asymptotic result in subsection 7.3 that at a fixed depth and over short time, heat conduction enhances the temperature increase of electromagnetic heating. We measure the relative reduction in reflex time by

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) is valid.

For small

which becomes a quadratic equation for

The second solution of the quadratic equation satisfies

Result (53) indicates that for small applied power density

We compare reflex times

In

For large

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 7.4. At reflex, the non-dimensional temperature has the expression

The no-conduction approximation of

which is independent of the applied power density

The situation of inferring parameters

As we discussed in section 4, constraint function

directly from the measured temperature profile at reflex

The non-dimensional version of (8) is Equation (28). Recall that (28) was derived as the governing equation for the reflex time

In particular, the non-dimensional beam spot area

independent variable

For simplicity and clarity, we use

With the proper notations clarified above, we write constraint (28) as

The true value

(57) is a constraint function in the form of

The solution

We examine constraint (57) in the form of

the parameter space of

inferring

operationally realistic since

density

In summary, to reliably determine model parameters

We consider a Gaussian beam with power density

where

where

To facilitate the analysis, we introduce parameters

First, we discuss how to generate simulated data of reflex time in the case of a Gaussian beam with finite radius. At time t, the activated domain is described by

Withdrawal reflex occurs when the volume of the activated domain reaches the critical threshold

We non-dimensionalize spatial coordinates

Note that

Equation (64) for

With these parameters, we solve for the true reflex time

regard the simulated

We study the methodology for determining the time delay in the observed withdrawal reflex:

In subsection 7.2, a joint constraint was constructed for 3 unknowns:

· the applied power density is uniform over the beam cross-section and is zero outside; heat conduction is included in the depth direction (case B);

· the beam spot area is large and approaching infinity (

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

The simulated curve of

Once we obtain the estimated values

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.

1) The prediction is based solely on the 3 data points. Formulation (44) and the predicted function

2) The prediction is invariant with respect to a shift in

In the equation above, all quantities on the right hand side are independent of

3) The prediction is invariant with respect to a scaling of

The predicted observed reflex time as a function of Q satisfies

Therefore, for the purpose of determining

It is worthwhile to point out that all of the properties above are attributed to the formulation form of (44). They are independent of the data on which algorithm (44) is applied. In particular, they are not affected by the true model governing the data generation in simulations or in experiments. Property 3 above is very powerful in applications. To estimate the time delay

Let

The proportionality constant

Let

We treat

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

In subsection 7.6, we discussed selecting test conditions for producing distinct constraint functions on

large

We use the temperature distribution given in (61). To mimic the situation of real applications, we work with

Here, for clarity, we reserve

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

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 [

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

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 dose-response 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

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

In addition to its high accuracy, formulation (44) has an important scaling property: the prediction by (44) is invariant with respect to a scaling in the independent variable

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

In the case of a Gaussian beam, we tested the strategy of selecting test conditions for determining parameters

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.

The authors thank Dr. John Biddle and Dr. Stacy Teng of the Institute for Defense Analyses for introducing the ADT CHEETEH-E model and for many helpful discussions. The authors acknowledge the Joint Non-Lethal Weapons Directorate of U.S. Department of Defense and the Naval Postgraduate School for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government.

The authors declare no conflicts of interest regarding the publication of this paper.

Wang, H.Y., Burgei, W.A. and Zhou, H. (2020) A Concise Model and Analysis for Heat-Induced Withdrawal Reflex Caused by Millimeter Wave Radiation. American Journal of Operations Research, 10, 31-81. https://doi.org/10.4236/ajor.2020.102004

In case A defined in (9), the temperature distribution

Integrating with respect to t yields the temperature distribution

Since case A is a special situation of case U, the reflex time

Solving for

Recall that

To solve for

Using the fundamental solution of the heat equation, we write

Completing the square in the exponent of term

Similarly, we can derive

Substituting

Recall that function

We carry out analysis in steps below to derive the asymptotics of

Step 1: We first show that

Noticing that

where

Using (68), we write out the expansion of

As t goes to infinity,

Step 2: Using the results of

We notice that u as a function of

Based on this, we write out an iterative formula for expanding

Starting the iteration with

Squaring both sides, we obtain the expansion of

Step 3: In this step, we expand

Step 4: In this step, we expand

Based on this, we write out an iterative formula for expanding

Starting the iteration with

Thus, near

Step 5: In this step, we expand

The expansion of

The insulated boundary condition