Effect of Depth-Dependent Nociceptor Density on the Heat-Induced Withdrawal Reflex

Hongyun Wang¹, Wesley A. Burgei², Hong Zhou^3*

¹Department of Applied Mathematics, University of California, Santa Cruz, CA, USA.
²U.S. Department of Defense, Joint Intermediate Force Capabilities Office, Quantico, VA, USA.
³Department of Applied Mathematics, Naval Postgraduate School, Monterey, CA, USA.

DOI: 10.4236/am.2020.118053                

Abstract

Previously we introduced a concise dose-response model for the heat-induced withdrawal reflex caused by millimeter wave radiation. The model predicts the occurrence of withdrawal reflex from the given spatial temperature profile. It was formulated on the assumption that the density of nociceptors in skin is uniform, independent of the depth. The model has only two parameters: the activation temperature of heat-sensitive nociceptors and the critical threshold on the activated volume for triggering withdrawal reflex. In this study, we consider the case of depth-dependent nociceptor density in skin. We use a general parametric form with a scaling parameter in the depth direction to represent the nociceptor density. We analyze system behaviors for four density types of this form. Based on the theoretical results, we develop a methodology for 1) identifying from test data the density form of nociceptors distribution, 2) finding from test data the scaling parameter in the density form, and 3) determining from test data the activation temperature of nociceptors.

Keywords

High-Energy Millimeter Wave Radiation, Heat-Induced Pain, Depth-Dependent Nociceptor Density In Skin

1. Introduction

The rapid development of applications such as wireless communications, security scanning, tissue diagnosis, and non-lethal weapons for crowd control or perimeter security has considerably increased human exposure to high-frequency millimeter waves (MMW) ranging from 30 to 300 gigahertz (GHz). For the purpose of biological risk assessment, it is vital to understand the effects of this irradiation on humans.

Many experiments have shown that exposure to MMW at sufficiently high intensities primarily produces fast heating near the surface of the skin [1] - [7]. The transmitted MMW power is absorbed in the skin to a depth of less than 0.5 mm at 100 GHz [8] and is attenuated exponentially as a function of skin depth. The skin generally consists of three different layers, namely epidermis, dermis, and hypodermis [9]. These layers have different thickness depending on the location of the skin. In particular, the epidermis is the outermost layer of skin containing both living and dead cells with thickness 0.075 - 0.15 mm. The dermis lies beneath the epidermis and is much thicker (1 - 4 mm). There are blood vessels and nerves in the dermis. The third layer is the hypodermis, which is composed of mainly subcutaneous fat. The hypodermis is about 1.1 - 5.6 mm in thickness. Studies performed at 60 GHz demonstrated that while the maximum value of the power density and specific absorption rate occurs at the epidermis, up to 60% of the incident power reaches the dermis, and only 10% gets to the hypodermis [10] [11]. Absorption of the MMW energy causes the local temperature of the skin to rise and can activate nociceptors [12] and consequently lead to a sensation of pain [13] [14].

Nociceptors are sensory nerve cells that respond to painful stimuli by sending out signals to the spinal cord via a chain of nerve fibers. When the collective signal becomes strong enough, the withdrawal reflex is triggered and the subject moves away from the exposure [15] [16]. Previously, we formulated a dose-response model for the heat-induced withdrawal reflex from MMV radiation [17]. The concise model predicts the occurrence of withdrawal reflex from a given spatial temperature profile of the skin. A prominent feature of the model is that it contains only two parameters. One key assumption in the concise model is that the nociceptor density in skin does not vary with the depth. In this paper we extend our earlier study by relaxing this assumption. Our goal is to determine the effect of depth-dependent nociceptor density on the heat-induced withdrawal reflex.

2. Mathematical Formulation in the Case of Depth Dependent Nociceptor Density

In this section, we study the mathematical formulation when the nociceptor density is a function of depth in the skin. We start by introducing proper mathematical notations:

· y: the depth coordinate (the skin surface is $y=0$ ).

· $r$ : 2-D coordinates on the skin surface; $\left(r\mathrm{,}y\right)$ is the 3-D coordinates.

· $T\left(r\mathrm{,}y\right)$ : 3-D spatial temperature profile of the skin.

· $T_{\text{act}}$ : activation temperature of nociceptors; given $T\left(r\mathrm{,}y\right)$, the activation status of a nociceptor at $\left(r\mathrm{,}y\right)$ is governed by the indicator function

$I_{\left(T\ge T_{\text{act}}\right)}\left(r\mathrm{,}y\right)\equiv \{\begin{array}{ll}\mathrm{1,}\hfill & \text{if}T\left(r\mathrm{,}y\right)\ge T_{\text{act}}\hfill \\ \mathrm{0,}\hfill & \text{if}T\left(r\mathrm{,}y\right)<T_{\text{act}}\hfill \end{array}$

· $\rho \left(y\right)$ : nociceptor density at depth y (number per volume).

· $X_{\text{act}}$ : total number of activated nociceptors in the skin,

$X_{\text{act}}={\displaystyle \int \rho \left(y\right)I_{\left(T\ge T_{\text{act}}\right)}\left(r\mathrm{,}y\right)\text{d}r\text{d}y}$

· ${\rho }_0$ : characteristic reference nociceptor density.

When the nociceptor density is uniform, $\rho \left(y\right)={\rho }_0$, the activated volume is proportional to the number of activated nociceptors. In this case, we adopted the activated volume as the single metric predictor variable (the input dose) for predicting withdrawal reflex [17].

Input dose in the case of uniform density:

$z\equiv {\displaystyle \int I_{\left(T\ge T_{\text{act}}\right)}\left(r\mathrm{,}y\right)\text{d}r\text{d}y}$ (1)

The advantages of this input dose are 1) it makes the dose quantity independent of the nociceptor density ${\rho }_0$, and 2) it shifts the effect of ${\rho }_0$ into the dose threshold $z_c$. In the dose response model, withdrawal reflex occurs when $z\ge z_c$ where the effects of ${\rho }_0$ and all other factors are reflected the single metric quantity $z_c$.

In the case of non-uniform nociceptor density, the activated volume is no longer proportional to the number of activated nociceptors, and thus is no longer a valid candidate for the input dose. We like to define the dose such that it contains (1) as a special case. For that purpose, we define the equivalent activated volume $z_{\text{equ}}$ based on reference density ${\rho }_0$, and use $z_{\text{equ}}$ as the input dose.

Input dose in the case of non-uniform density:

$z_{\text{equ}}\equiv \frac{1}{{\rho }_0}{X}_{\text{act}}={\displaystyle \int }\frac{\rho \left(y\right)}{{\rho }_0}{I}_{\left(T\ge T_{\text{act}}\right)}\left(r\mathrm{,}y\right)\text{d}r\text{d}y$ (2)

The dose response relation has the same form as in the case of uniform density.

$\text{Outcome}\left(z_{\text{equ}}\right)=\{\begin{array}{ll}1\left(\text{withdrawal}\text{reflex}\right)\mathrm{,}\hfill & \text{if}{z}_{\text{equ}}\ge z_c\hfill \\ 0\left(\text{no}\text{withdrawal}\text{reflex}\right)\mathrm{,}\hfill & \text{if}{z}_{\text{equ}}<z_c\hfill \end{array}$

In this study, we consider the hypothetical situation where the time of withdrawal reflex and the spatial temperature profile at reflex are measurable in experiments. With these two measurable entities, we explore the behaviors of several parameterized nociceptor density types. The objectives of the study are 1) to distinguish these candidate density types from each other based on the measurable entities, and 2) to infer the parameter values.

The calculation of $z_{\text{equ}}$ defined in (2) requires only the relative density $\rho \left(y\right)/{\rho }_0$. The effect of ${\rho }_0$ is contained in the dose threshold $z_c$. When $\rho \left(y\right)/{\rho }_0$ is given, the dose response model has only two unknown parameters: $T_{\text{act}}$ and $z_c$. In a test, the measured temperature profile $T\left(r\mathrm{,}y\right)$ at reflex provides a constraint on $\left(T_{\text{act}}\mathrm{,}{z}_c\right)$. Mathematically, we construct constraint function $z_c\left(T_{\text{act}}\right)$ as follows. For any value of $T_{\text{act}}$, by definition, the corresponding value of $z_c$ is the value of $z_{\text{equ}}$ calculated based on $T_{\text{act}}$ and $T\left(r\mathrm{,}y\right)$.

$z_c\left(T_{\text{act}}\right)={\displaystyle \int }\frac{\rho \left(y\right)}{{\rho }_0}{I}_{\left(T\ge T_{\text{act}}\right)}\left(r\mathrm{,}y\right)\text{d}r\text{d}y$ (3)

when $\rho \left(y\right)/{\rho }_0$ is given, function $z_c\left(T_{\text{act}}\right)$ is completely determined by the measured $T\left(r\mathrm{,}y\right)$. Constraint functions $z_c\left(T_{\text{act}}\right)$ based on $T\left(r\mathrm{,}y\right)$ measured at different test conditions are potentially distinct from each other. All these constraint functions have one common intersection, which gives the true values of $\left(T_{\text{act}}\mathrm{,}{z}_c\right)$.

When the true relative density $\rho \left(y\right)/{\rho }_0$ is unknown, we work with a trial relative density $r^{\left(\text{try}\right)}\left(y\right)$. We use $r^{\left(\text{try}\right)}\left(y\right)$ to replace $\rho \left(y\right)/{\rho }_0$ in (3) and construct trial constraint function $z_c\left(T_{\text{act}}\right)$ from the test data. Note that the true values of $\left(T_{\text{act}}\mathrm{,}{z}_c\right)$ satisfy only the true constraint function calculated using the true $\rho \left(y\right)/{\rho }_0$. When $r^{\left(\text{try}\right)}\left(y\right)$ deviates from the true $\rho \left(y\right)/{\rho }_0$, the true values of $\left(T_{\text{act}}\mathrm{,}{z}_c\right)$ are not on the trial constraint curve calculated using $r^{\left(\text{try}\right)}\left(y\right)$. Consequently, for a pair of trial constraint functions calculated using $r^{\left(\text{try}\right)}\left(y\right)$ (based on measured $T\left(r\mathrm{,}y\right)$ of two distinct test conditions), their intersection is not at the true values of $\left(T_{\text{act}}\mathrm{,}{z}_c\right)$, and the intersection varies with the test conditions of the pair. The test-condition-dependence of the intersection serves as an indication that the trial density $r^{\left(\text{try}\right)}\left(y\right)$ is incorrect. The specific behavior of test-condition-dependence of the intersection provides a venue for us to tune $r^{\left(\text{try}\right)}\left(y\right)$ toward the true $\rho \left(y\right)/{\rho }_0$.

We examine several types of parameterized density. We study the test-condition dependence of 1) the reflex time and 2) the intersection of a pair of trial constraint functions. The goal is to identify system behaviors that a) help us distinguish these density types from each other and b) guide us to tune the trial parameter toward its true value.

We carry out the analysis in the idealized situation where

· the electromagnetic heating is uniform over the beam cross-section (with area A) and it decays exponentially with depth y;

· the initial temperature is uniform everywhere;

· the heat conduction is included only in the depth direction.

This is the same as case B in our previous study [17]. At any given time, the temperature inside the beam cross-section is a function of depth y only and it decreases with y. The time evolution of temperature distribution $T\left(y,t\right)$ is governed by

$\{\begin{array}{l}{\rho }_{\text{m}}{C}_p\frac{\partial T\left(y\mathrm{,}t\right)}{\partial t}=K\frac{{\partial }^{2}T\left(y\mathrm{,}t\right)}{\partial y^2}+P_{\text{dep}}\mu {\text{e}}^{-\mu y}\hfill \\ {\frac{\partial T\left(y\mathrm{,}t\right)}{\partial y}|}_{y=0}=0\mathrm{,}T\left(y\mathrm{,0}\right)=T_0\hfill \end{array}$

where

· ${\rho }_{\text{m}}$ is the mass density of skin,

· $C_p$ is the specific heat capacity of the skin,

· K is the thermal conductivity of the skin,

· $\mu$ is the absorption coefficient of the skin,

· $P_{\text{dep}}$ is the beam power density deposited on (absorbed into) the skin, and

· $T_0$ is the initial temperature of the skin.

In this idealized situation, the region of activated nociceptors ( $T\left(y\right)\ge T_{\text{act}}$ ) is a cylinder with depth $y_{\text{act}}$ governed by $T\left(y_{\text{act}}\right)=T_{\text{act}}$. By definition, $T_{\text{act}}$, ${y_{\text{act}}|}_{\left(\text{reflex}\right)}$ and the reflex time $t_{\text{ref}}$ are constrained by the temperature distribution, which we assume is measurable.

$T_{\text{act}}=T\left({y_{\text{act}}|}_{\left(\text{reflex}\right)}\mathrm{,}{t}_{\text{ref}}\right)$ (4)

We consider a general 1-parameter form for the relative density of nociceptors

$\frac{\rho \left(y\right)}{{\rho }_0}=f(\; \beta \; y\; )$

where $f\left(s\right)$ can be any positive function. The activated depth $y_{\text{act}}$ and the dose $z_{\text{equ}}$ both vary with the activation temperature $T_{\text{act}}$, and are related by Equation (2) as

$z_{\text{equ}}=A{\displaystyle {\int }_0^{y_{\text{act}}}}f\left(\beta y\right)\text{d}y=\frac{A}{\beta }F\left(\beta y_{\text{act}}\right)$ (5)

$y_{\text{act}}=\frac{1}{\beta }{F}^{-1}\left(\frac{\beta z_{\text{equ}}}{A}\right)$ (6)

where A is the beam spot area, $F\left(s\right)\equiv {\displaystyle {\int }_0^s}f\left(z\right)\text{d}z$, and $F^{-1}\left(u\right)$ is the inverse function of $F\left(s\right)$. Since $f\left(s\right)$ is positive, function $F\left(s\right)$ is monotonically increasing and the inverse function $F^{-1}\left(u\right)$ is well-defined over the range of $F\left(s\right)$. By definition, functions $F\left(s\right)$ and $F^{-1}\left(u\right)$ satisfy $F\left(0\right)=0$ and $F^{-1}\left(0\right)=0$. Recall that the dose threshold is defined as $z_c\equiv {z_{\text{equ}}|}_{\left(\text{reflex}\right)}$. Equation (6) gives

${y_{\text{act}}|}_{\left(\text{reflex}\right)}=\frac{1}{\beta }{F}^{-1}\left(\frac{\beta z_c}{A}\right)$ (7)

Equation (4) in combination with (7) provides a constraint on $T_{\text{act}}$, $z_c$, $\beta$ and $t_{\text{ref}}$, which can be used for different purposes, depending on which parameters are known. When $t_{\text{ref}}$ and $T\left(y\mathrm{,}t\right)$ are measured, (4) gives us a constraint on $z_c$, $T_{\text{act}}$ and $\beta$. On the other hand, when $z_c$, $T_{\text{act}}$ and $\beta$ are given, (4) can be viewed as a governing equation for $t_{\text{ref}}$. This is useful, for example, for examining the behavior of $t_{\text{ref}}$ vs. A.

In the analysis of subsequent sections, we need the expansions of $f\left(s\right)$ and $F\left(s\right)$, $F^{-1}\left(u\right)$ and their derivatives. We now derive these expansions. We first write out the Taylor expansion of $f\left(s\right)$ around $s=0$.

$f\left(s\right)=a_0+a_{1}s+a_2\frac{s^2}{2!}+\cdots \text{where}{a}_k=f^{\left(k\right)}(\; 0\; )$

$F\left(s\right)\equiv {\displaystyle {\int }_0^s}f\left(z\right)\text{d}z=a_{0}s+a_1\frac{s^2}{2!}+a_2\frac{s^3}{3!}+\cdots$

The expansion of $F^{-1}\left(u\right)$ is derived from that of $F\left(s\right)$ using an iterative method. It depends on which term in the expansion of $F\left(s\right)$ is the leading non-zero term. Let $u\equiv F\left(s\right)$. The inverse function $F^{-1}\left(u\right)$ maps u back to s. We discuss two cases.

· Case 1: $f\left(0\right)\ne 0$. Function $u=F\left(s\right)$ has the expansion

$u=a_{0}s+a_1\frac{s^2}{2!}+\cdots \text{where}{a}_k=f^{\left(k\right)}(\; 0\; )$

Based on that, we built an iteration formula:

$s_{j+1}=\frac{1}{a_0}\left(u-a_1\frac{1}{2!}{s}_j^2\right),s_0=0$

The iteration gives us expansions of $F^{-1}\left(u\right)$ and ${\left(F^{-1}\right)}^{\prime }(\; u\; )$

$F^{-1}\left(u\right)=\frac{1}{a_0}u-\frac{a_1}{2a_0^3}{u}^2+\cdots \text{for}f\left(0\right)\ne 0$ (8)

${\left(F^{-1}\right)}^{\prime }\left(u\right)=\frac{1}{a_0}-\frac{a_1}{a_0^3}u+\cdots \text{for}f\left(0\right)\ne 0$ (9)

· Case 2: $f\left(0\right)=0$ but $f^{\prime }\left(0\right)\ne 0$. Function $u=F\left(s\right)$ has the expansion

$u=a_1\frac{s^2}{2!}+a_2\frac{s^3}{3!}+\cdots \text{where}{a}_k=f^{\left(k\right)}(\; 0\; )$

Based on that, we construct an iteration formula for $s^2$ :

$s_{j+1}^2=\frac{2}{a_1}\left(u-a_2\frac{1}{3!}{\left(s_j^2\right)}^{3/2}\right),s_0^2=0$

which yields the expansion of $s^2$ in terms of u

$s^2=\frac{2}{a_1}u\left(1-\frac{a_2}{6}\cdot \frac{2}{a_1}{\left(\frac{2}{a_1}u\right)}^{1/2}+\cdots \right)$

Taking square roots of both sides, we obtain

$F^{-1}\left(u\right)=\frac{\sqrt{2}}{\sqrt{a_1}}{u}^{1/2}-\frac{a_2}{3a_1^2}u+\cdots \text{for}f\left(0\right)=0\text{but}{f}^{\prime }\left(0\right)\ne 0$ (10)

${\left(F^{-1}\right)}^{\prime }\left(u\right)=\frac{1}{\sqrt{2a_1}}{u}^{-1/2}-\frac{a_2}{3a_1^2}+\cdots \text{for}f\left(0\right)=0\text{but}{f}^{\prime }\left(0\right)\ne 0$ (11)

with the mathematical results above, we study the behavior of $t_{\text{ref}}$ vs A.

3. Analysis of Reflex Time vs. Beam Spot Area

When $T_{\text{act}}$, $z_c$ and $\beta$ are given, (4) with (7) governs $t_{\text{ref}}$ vs A. In our previous study [17], we scaled and shifted the physical temperature distribution to the normalized non-dimensional temperature $H\left(y_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$.

$\frac{\left(T\left(y\mathrm{,}t\right)-T_0\right)K\mu }{P_{\text{dep}}}=H\left(y_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$ (12)

where the non-dimensional depth $y_{\text{nd}}$ and the non-dimensional time $t_{\text{nd}}$ are defined as

$y_{\text{nd}}\equiv \mu y\mathrm{,}{t}_{\text{nd}}\equiv t\frac{K{\mu }^2}{{\rho }_{\text{m}}{C}_p}$

The normalized temperature $H\left(y_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$ has the expression

$H\left(y\mathrm{,}t\right)\equiv {\displaystyle {\int }_0^t}G\left(y\mathrm{,}s\right)\text{d}s$

$G\left(y\mathrm{,}t\right)\equiv \frac{1}{2}\text{erfc}\left(\frac{2t-y}{\sqrt{4t}}\right){\text{e}}^{t-y}+\frac{1}{2}\text{erfc}\left(\frac{2t+y}{\sqrt{4t}}\right){\text{e}}^{t+y}$ (13)

It is important to notice that the normalized temperature $H\left(y\mathrm{,}t\right)$ is parameter-free. The non-dimensional version of (4) with (7) has the form

$\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}=H\left(y_{\text{act},\text{nd}}\mathrm{,}{t}_{\text{ref},\text{nd}}\right)$ (14)

$y_{\text{act},\text{nd}}\equiv \frac{\mu }{\beta }{F}^{-1}\left(\frac{\beta }{\mu A_{\text{nd}}}\right)\mathrm{,}{A}_{\text{nd}}\equiv \frac{A}{\mu z_c}\mathrm{,}{t}_{\text{ref},\text{nd}}\equiv t_{\text{ref}}\frac{K{\mu }^2}{{\rho }_{\text{m}}{C}_p}$

In the above, $A_{\text{nd}}\equiv \frac{A}{\mu z_c}$ is the non-dimensional beam spot area. In the limit of $A\to \infty$, we have $y_{\text{act},\text{nd}}\to 0$ and the equation for $t_{\text{ref},\text{nd}}$ becomes

$\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}\mathrm{<mo>\; =\; </mo>}h\left({t_{\text{ref},\text{nd}}|}_{A\to \infty }\right)$

where $h\left(t_{\text{nd}}\right)\equiv H\left(\mathrm{0,}{t}_{\text{nd}}\right)$. Let $t_0\equiv h^{-1}\left(\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}\right)$. We have ${t_{\text{ref},\text{nd}}|}_{A\to \infty }=t_0$. We examine the asymptotic behavior of this convergence. We seek an expansion of the form

$t_{\text{ref},\text{nd}}\left(A_{\text{nd}}\right)=t_0\left[1+c_A{\left(\frac{1}{A_{\text{nd}}}\right)}^{\alpha }+\cdots \right]$ (15)

Expanding $H\left(y_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$ around $\left(\mathrm{0,}{t}_0\right)$ and substituting (15), we get

$\mathrm{0\; <mo>\; =\; </mo>}\frac{\partial H\left(\mathrm{0,}{t}_0\right)}{\partial t}{t}_0{c}_A{\left(\frac{1}{A_{\text{nd}}}\right)}^{\alpha }+\frac{1}{2}\frac{{\partial }^{2}H\left(\mathrm{0,}{t}_0\right)}{\partial y^2}{\left(\frac{\mu }{\beta }{F}^{-1}\left(\frac{\beta }{\mu A_{\text{nd}}}\right)\right)}^2$ (16)

Exponent $\alpha$ and coefficient $c_A$ are determined using the leading term expansion of $F^{-1}\left(u\right)$. The result depends on whether the nociceptor density at the skin surface is zero.

· Case 1: $f\left(0\right)\ne 0$.

The expansion of $F^{-1}\left(u\right)$ given in (8). Substituting it into (16), we have

$\alpha =\mathrm{2,}{c}_A=\frac{1}{2t_0{\left(f\left(0\right)\right)}^2}\cdot \frac{1-\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}{\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}$ (17)

Here we have used $\frac{\partial H\left(\mathrm{0,}{t}_0\right)}{\partial t}=\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}$ and $\frac{{\partial }^{2}H\left(\mathrm{0,}{t}_0\right)}{\partial y^2}=\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}-1$ derived in our previous study.

· Case 2: $f\left(0\right)=0$ but $f^{\prime }\left(0\right)\ne 0$

The expansion of $F^{-1}\left(u\right)$ given in (10). Substituting it into (16), we arrive at

$\alpha =\mathrm{1,}{c}_A=\frac{\mu }{t_0{f}^{\prime }\left(0\right)\beta }\cdot \frac{1-\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}{\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}$ (18)

Returning to the physical quantities before the non-dimensionalization, the reflex time vs beam spot area is given by

$t_{\text{ref}}\left(A\right)=\frac{t_0{\rho }_{\text{m}}{C}_p}{K{\mu }^2}\left[1+c_A{\left(\frac{\mu z_c}{A}\right)}^{\alpha }+\cdots \right]\mathrm{,}{t}_0=h^{-1}\left(\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}\right)$ (19)

4. Analysis of Constraint Functions on $\left(T_{act},z_c\right)$

When $t_{\text{ref}}$ and $T\left(y\mathrm{,}{t}_{\text{ref}}\right)$ are measured, (4) with (7) serves as a constraint on $T_{\text{act}}$, $z_c$ and $\beta$ with beam spot area A as a parameter describing the test condition. We denote the constraint function as $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and use (4) to write it as

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)=\Phi \left(\frac{1}{\beta }{F}^{-1}\left(\frac{\beta z_c}{A}\right)\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ (20)

where $\Phi \left(y\mathrm{<mo>\; ;\; </mo>}v\right)T\left(y\mathrm{,}{t}_{\text{ref}}\left(1/v\right)\right)$ is the temperature profile at reflex with beam spot area $A=1/v$. We represent the effect of A via variable $v=1/A$ since $t_{\text{ref}}$ is a smooth function of $v=1/A$ as $A\to \infty$. To facilitate the discussion, we introduce two sets of mathematical notations. These two sets of notations are used to distinguish a quantity’s true value from its role as a variable in a function.

· ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ : true value of coefficient $\beta$.

· $\beta$ : a trial value of coefficient $\beta$.

· $\left(T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{z}_c^{\mathrm{<mo>\; *\; </mo>}}\right)$ : true values of model parameters $T_{\text{act}}$ and $z_c$.

· $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ : constraint function (20) calculated using trial value $\beta$ ; in $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$, $z_c$ denotes the independent variable, and $T_{\text{act}}$ the dependent variable.

Note that the data is generated with the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$. In the calculation constraint function (20), we use the measured data and a trial value $\beta$ since ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ is unknown. When $\beta ={\beta }^{\mathrm{<mo>\; *\; </mo>}}$, the constraint function shares one common intersection for all values of A:

$T_{\text{act}}\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\text{for}\text{all}A$ (21)

when $\beta \ne {\beta }^{\mathrm{<mo>\; *\; </mo>}}$, in general $\left(T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{z}_c^{\mathrm{<mo>\; *\; </mo>}}\right)$ is not on constraint curve $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$, and the intersection of a pair of constraint functions varies with the test conditions (values of A). We study the behavior of the intersection vs A. We make use of (21) and expand $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ around $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)$. For conciseness, we introduce $\epsilon \equiv \frac{\beta z_c}{A}$ and write (20) as:

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)=\Phi \left(\xi \mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\mathrm{,}\xi \left(\beta \mathrm{,}\epsilon \right)\equiv \frac{1}{\beta }{F}^{-1}\left(\epsilon \right)\mathrm{,}\epsilon \left(\beta \mathrm{,}{z}_c\right)\equiv \frac{\beta z_c}{A}$

We apply the chain rule to calculate partial derivatives.

$\frac{\partial T_{\text{act}}}{\partial z_c}=\frac{\partial }{\partial y}\Phi \left(\xi \mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right){\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\frac{1}{A}$ (22)

$\begin{array}{c}\frac{\partial T_{\text{act}}}{\partial \beta }=\frac{\partial }{\partial y}\Phi \left(\xi \mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right){\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\frac{1}{A}\cdot \frac{z_c}{\beta }\left[-\frac{F^{-1}\left(\epsilon \right)}{\epsilon {\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)}+1\right]\\ \equiv \frac{\partial T_{\text{act}}}{\partial z_c}\cdot \eta \text{where}\eta \equiv \frac{z_c}{\beta }\left[-\frac{F^{-1}\left(\epsilon \right)}{\epsilon {\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)}+1\right]\end{array}$ (23)

Using these derivatives, we write out the expansion of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$.

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}+{\frac{\partial T_{\text{act}}}{\partial z_c}|}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}\cdot \left[\left(z_c-z_c^{\mathrm{<mo>\; *\; </mo>}}\right)+{\eta |}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)\right]$ (24)

The slope of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ vs. $z_c$ is ${\frac{\partial T_{\text{act}}}{\partial z_c}|}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}$, which is given in (22). In Appendix A, we show that the slope converges to zero as $A\to \infty$. It follows that

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A=\infty \right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\text{for}\text{all}\beta$ (25)

(25) shows that in the limit of $A\to \infty$, the constraint curve is a horizontal line at $T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$, independent of $z_c$ and $\beta$. For finite A, $\frac{\partial T_{\text{act}}}{\partial z_c}\ne 0$ and $T_{\text{act}}$ in (24) varies with $z_c$ and $\beta$. We consider the intersection of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. The $T_{\text{act}}$ -coordinate of the intersection is $T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. Let $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ denote the $z_c$ -coordinate of the intersection. Solving for $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ from (24) and (25) leads to

$z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}-{\eta |}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)$ (26)

The dependence of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ on A is contained in ${\eta |}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}$ given in (23). Using the expansion of $\frac{F^{-1}\left(\epsilon \right)}{\epsilon {\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)}$ derived in (73) in Appendix A, we investigate the behavior of $\eta$ at large A.

· Case 1: $f\left(0\right)\ne 0$.

${\eta |}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}={\frac{z_c}{\beta }\left(-\frac{a_1}{2a_0^2}\epsilon \right)|}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}=\frac{-{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\right)}^2{f}^{\prime }\left(0\right)}{2{\left(f\left(0\right)\right)}^2}\cdot \frac{1}{A}\text{for}\text{large}A$

In case 1, as $A\to \infty$, the intersection $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ given by (26) converges to the true value $z_c^{\mathrm{<mo>\; *\; </mo>}}$ regardless of trial value $\beta$.

$\underset{A\to \infty }{lim}{z}_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}\text{for}\text{all}\beta$

The residual in convergence is proportional to $\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)/A$ and has the same sign as $f^{\prime }\left(0\right)\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)$. More specifically, we have

$z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}+\frac{{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\right)}^2{f}^{\prime }\left(0\right)}{2{\left(f\left(0\right)\right)}^2}\cdot \frac{\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{A}\text{for}\text{large}A$ (27)

· Case 2: $f\left(0\right)=0$ but $f^{\prime }\left(0\right)\ne 0$

$\begin{array}{c}{\eta |}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}={\frac{z_c}{\beta }\left(-1-\frac{\sqrt{2}{a}_2}{3a_1^{3/2}}{\epsilon }^{1/2}\right)|}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}\\ =\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}\left(-1-\frac{\sqrt{2{\beta }^{\mathrm{<mo>\; *\; </mo>}}{z}_c^{\mathrm{<mo>\; *\; </mo>}}}{f}^{″}\left(0\right)}{3{\left(f^{\prime }\left(0\right)\right)}^{3/2}}\cdot \frac{1}{\sqrt{A}}\right)\to \frac{-z_c^{\mathrm{<mo>\; *\; </mo>}}}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}\text{for}\text{large}A\end{array}$ (28)

In case 2, as $A\to \infty$, the intersection $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ given by (26) converges to

$\underset{A\to \infty }{lim}{z}_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}-\left(\frac{-z_c^{\mathrm{<mo>\; *\; </mo>}}}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}\right)\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}\frac{\beta }{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}$

which is proportional to trial value $\beta$. The residual in convergence is proportional to $\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)/\sqrt{A}$ and has the same sign as $f^{″}\left(0\right)\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)$. Asymptotically, $z_c^{\left(I\right)}$ is

$z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}\frac{\beta }{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}+\frac{{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\right)}^{3/2}}{{\left({\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}^{1/2}}\cdot \frac{\sqrt{2}{f}^{″}\left(0\right)}{3{\left(f^{\prime }\left(0\right)\right)}^{3/2}}\cdot \frac{\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{\sqrt{A}}\text{for}\text{large}A$ (29)

with the analytical preparations above, we examine four types of parametric form for nociceptor density vs depth, depicted in Figure 17.

5. Type 1 Nociceptor Density: $\rho \left(y\right)={\rho }_0{e}^{-\beta y}$

For type 1 nociceptor density, the relative density takes the parametric form

$\frac{\rho \left(y\right)}{{\rho }_0}=f\left(\beta y\right)\mathrm{,}\text{where}f\left(s\right)={\text{e}}^{-s}$ (30)

The graph of type 1 $f\left(s\right)$ is plotted in Figure 17. It has the properties

$f\left(0\right)=1\ne 0,f^{\prime }\left(0\right)=-1,\underset{0\le s<\infty }{\max }f\left(s\right)=1$

5.1. Reflex Time

The reflex time vs beam spot area is given by (19) and (17), namely,

$t_{\text{ref}}\left(A\right)=\frac{t_0{\rho }_{\text{m}}{C}_p}{K{\mu }^2}\left[1+c_A{\left(\frac{\mu z_c}{A}\right)}^2+\cdots \right]\mathrm{,}{c}_A=\frac{1-\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}{2t_0\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}$ (31)

As A increases, $t_{\text{ref}}$ converges to its limit with the residual proportional to 1/A²:

$t_{\text{ref}}\left(A\right)-t_{\text{ref}}\left(\infty \right)~\frac{1}{A^2}$

Figure 1 plots the relation between $t_{\text{ref}}$ and A in two ways. Left panel: $t_{\text{ref}}$ vs A. Right panel: $t_{\text{ref}}$ vs. 1/A. In particular, the right panel confirms that the residual in the convergence of $t_{\text{ref}}$ decays faster than 1/A for large A, as predicted in the analysis above.

Figure 1. The relation between reflex time ( $t_{\text{ref}}$ ) and beam spot area (A) for type 1 nociceptor density (30) with ${\beta }^{\mathrm{<mo>\; *\; </mo>}}=1$. Left panel: $t_{\text{ref}}$ vs A. Right panel: $t_{\text{ref}}$ vs 1/A.

It is interesting to notice that expansion (31) is independent of ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$. Consequently, the measurements of $t_{\text{ref}}\left(A\right)$ vs. A do not contain any information for estimating ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$.

5.2. Constraint Function $T_{act}\left(z_c;\beta ,A\right)$

We consider the intersection of the pair $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. The $z_c$ -coordinate of the intersection, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$, is generally described by (26). For type 1 density, we have $f\left(s\right)={\text{e}}^{-s}$, $f\left(0\right)=1\ne 0$ and $f^{\prime }\left(0\right)=-1$, and the specific expression of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ is given by (27).

$z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}-\frac{{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\right)}^2}{2}\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)\cdot \frac{1}{A}$ (32)

At $\beta ={\beta }^{*}$, the intersection $z_c^{\left(I\right)}\left({\beta }^{*},A\right)=z_c^{*}$ is independent of A. When $\beta \ne {\beta }^{\mathrm{<mo>\; *\; </mo>}}$, the trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs. A tells us whether $\beta >{\beta }^{*}$ or $\beta <{\beta }^{*}$.

· For $\beta >{\beta }^{*}$, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ ascends toward $z_c^{\mathrm{<mo>\; *\; </mo>}}$ from below as A increases.

· For $\beta <{\beta }^{*}$, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ descends toward $z_c^{\mathrm{<mo>\; *\; </mo>}}$ from above as A increases.

Figure 2 displays simulated $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ for several values of A, respectively for $\beta >{\beta }^{*}$ and for $\beta <{\beta }^{*}$. Here constraint function $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ is based on test data (which is generated with true value ${\beta }^{*}=1$ ) and is calculated using formulation (20) with trial value $\beta$. The trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs A is illustrated in Figure 3. The simulation results in Figure 2 and Figure 3 confirm the theoretically predicted trend above.

When it is known that the nociceptor density has the parametric form of type 1 given in (30), we can tune the trial value $\beta$ down or up toward the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ depending on whether the calculated $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ increases or decreases with A.

5.3. Constraint Function Calculated Using the Uniform Density

We now consider the situation where the type of parametric form of $\rho \left(y\right)/{\rho }_0$

Figure 2. Constraint curves for type 1 density (30). $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ is based on test data generated with true value ${\beta }^{*}=1$, and calculated using (20) with trial value $\beta$. Left panel: $\beta >{\beta }^{*}$. Right panel: $\beta <{\beta }^{*}$.

Figure 3. $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs. A, respectively for $\beta >{\beta }^{*}$ and for $\beta <{\beta }^{*}$. Here $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ is the intersection of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ from Figure 2.

is unknown. With no information on the type of density form, we use the uniform density as the trial density in calculating the constraint function. Let $T_{\text{act},\text{uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ denote the constraint function based on the test data (which is generated using type 1 density (30) with true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ ) and calculated using framework (20) with the uniform trial density $\rho \left(y\right)/{\rho }_0\equiv 1$. Notice that the uniform density is a member of type 1 family (30) with $\beta =0$. Thus, the two constraint functions $T_{\text{act},\text{uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ are related by $T_{\text{act},\text{uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)=T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta =\mathrm{0,}A\right)$. Let $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ denote the $z_c$ -coordinate of the intersection of the pair $T_{\text{act},\text{uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A=\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. Setting $\beta =0$ in (32), we obtain

$z_{c,\text{uni}}^{\left(I\right)}\left(A\right)=z_c^{*}+\frac{{\left(z_c^{*}\right)}^2}{2}{\beta }^{*}\cdot \frac{1}{A}$ (33)

Since the uniform density is parameter-free, the calculation of constraint function $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and intersection $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is based solely on the test data. It does not require any input parameter or knowledge of the function form of the true density. Once the test data is available, $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ can be calculated. When the true underlying density is type 1 given in (30), result (33) predicts that $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ converges to $z_c^{\mathrm{<mo>\; *\; </mo>}}$ as $A\to \infty$ with the difference proportional to 1/A. Figure 4 compares simulated $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ and $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$. The simulation results validate the theoretical prediction. In particular, $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ varies linearly with 1/A. We fit function $c_0+c_1/A$ to data of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs 1/A. The fitting coefficients give us

$z_c^{*}=c_0,{\beta }^{*}=\frac{2c_1}{c_0^2}$

Before we end this section, we clarify that result (33) predicts the behavior of constraint function $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ calculated using the uniform trial density when the true underlying density affecting the test data is type 1 given in (30). When the true underlying density is of a different type, the behavior will be different. One objective of examining $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is to identity the type of nociceptor density from the observed behavior of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs A, based on the theoretically predicted behaviors for a list of density types. In the subsequent sections, we will study more density types.

6. Type 2 Nociceptor Density: $\rho \left(y\right)={\rho }_0\cdot \left(\beta y\right)e^{1-\beta y}$

For type 2 nociceptor density, the relative density has the parametric form

$\frac{\rho \left(y\right)}{{\rho }_0}=f\left(\beta y\right)\mathrm{,}\text{where}f\left(s\right)=s{\text{e}}^{1-s}$ (34)

Figure 4. Simulated results of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ for type 1 density given in (30). Here $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is based on the data from the true density (with ${\beta }^{*}=1$ ) but is calculated using the uniform trial density. Left panel: $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs A. Right panel: $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. 1/A.

The graph of type 2 $f\left(s\right)$ is shown in Figure 17. It is straightforward to see that

$f\left(0\right)=0,f^{\prime }\left(0\right)=\text{e}\ne 0,f^{″}\left(0\right)=-2\text{e},\underset{0\le s<\infty }{\max }f\left(s\right)=1$

6.1. Reflex Time

The reflex time vs beam spot area is described by (19) and (18).

$t_{\text{ref}}\left(A\right)=\frac{t_0{\rho }_{\text{m}}{C}_p}{K{\mu }^2}\left[1+c_A\cdot \frac{\mu z_c}{A}+\cdots \right]\mathrm{,}{c}_A=\frac{\mu \left(1-\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}\right)}{t_0\text{e}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\text{erfc}\left(\sqrt{t_0}\right){\text{e}}^{t_0}}$ (35)

As A increases, $t_{\text{ref}}$ converges to its limit with the difference proportional to 1/A.

$t_{\text{ref}}\left(A\right)-t_{\text{ref}}\left(\infty \right)~\frac{1}{A}$

Figure 5 plots the relation between $t_{\text{ref}}$ and A in two ways. Left panel: $t_{\text{ref}}$ vs A. Right panel: $t_{\text{ref}}$ vs 1/A. In particular, the right panel confirms that $t_{\text{ref}}$ is linear with respect to 1/A for large A, as predicted in the analysis above. This is in contrast with the convergence of 1/A² for type 1 nociceptor density (30). To distinguish between these two density forms, we introduce an auxiliary quantity Q

$Q\equiv \frac{t_{\text{ref}}\left(A\right)-t_{\text{ref}}\left(2A\right)}{t_{\text{ref}}\left(2A\right)-t_{\text{ref}}\left(4A\right)}$ (36)

The theoretical prediction above tells us

$Q=\{\begin{array}{ll}4\hfill & \text{if}\rho \left(y\right)={\rho }_0{\text{e}}^{-\beta y}\left(\text{type}1\right)\hfill \\ 2\hfill & \text{if}\rho \left(y\right)={\rho }_0\cdot \left(\beta y\right){\text{e}}^{1-\beta y}\left(\text{type}2\right)\hfill \end{array}$ (37)

In (35), coefficient $c_A$ does contain ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$. However, in (35) ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ is tangled with

Figure 5. The relation between reflex time ( $t_{\text{ref}}$ ) and beam spot area (A) for type 2 nociceptor density (34) with ${\beta }^{\mathrm{<mo>\; *\; </mo>}}=1$. Left panel: $t_{\text{ref}}$ vs A. Right panel: $t_{\text{ref}}$ vs 1/A.

other parameters. It is not possible to extract the value of ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ unless all other parameters are known.

6.2. Constraint Function $T_{act}\left(z_c;\beta ,A\right)$

We consider the intersection of the pair $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. The $z_c$ -coordinate of the intersection, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$, is generally described by (26). For type 2 density, we have $f\left(s\right)=s{\text{e}}^{1-s}$, $f\left(0\right)=0$, $f^{\prime }\left(0\right)=\text{e}\ne 0$ and $f^{″}\left(0\right)=-2\text{e}$, and the specific expression of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ is given by (29).

$z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}\frac{\beta }{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}-\frac{{\left(2z_c^{\mathrm{<mo>\; *\; </mo>}}\right)}^{3/2}}{3\sqrt{\text{e}{\beta }^{\mathrm{<mo>\; *\; </mo>}}}}\left(\beta -{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)\cdot \frac{1}{\sqrt{A}}$

At $\beta ={\beta }^{\mathrm{<mo>\; *\; </mo>}}$, the intersection $z_c^{\left(I\right)}\left({\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}$ is independent of beam spot area A. When $\beta \ne {\beta }^{\mathrm{<mo>\; *\; </mo>}}$, the trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs A tells us whether $\beta >{\beta }^{*}$ or $\beta <{\beta }^{*}$.

· For $\beta >{\beta }^{*}$, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ ascends toward $z_c^{\mathrm{<mo>\; *\; </mo>}}\frac{\beta }{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}$ as A increases.

· For $\beta <{\beta }^{*}$, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ descends toward $z_c^{\mathrm{<mo>\; *\; </mo>}}\frac{\beta }{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}$ as A increases.

The increase/decrease trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs. A is qualitatively the same as that for type 1 density given in (30). There are two differences. For type 2 density (34), we have 1) $li{m}_{A\to \infty }{z}_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}\frac{\beta }{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}$, which varies with trial value $\beta$, and 2) the residual in convergence is proportional to $1/\sqrt{A}$, instead of 1/A.

Figure 6 shows simulated $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ for several values of A, respectively for $\beta >{\beta }^{*}$ and for $\beta <{\beta }^{*}$. Here constraint function $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ is based on test data (which is generated with true value ${\beta }^{*}=1$ ) and is calculated using formulation (20) with trial value $\beta$. The trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs A is shown in Figure 7. The simulation results in Figure 6 and Figure 7 confirm the theoretically predicted trend above.

Figure 6. Constraint curves for type 2 density (34). $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ is based on test data generated with true value ${\beta }^{*}=1$, and calculated using (20) with trial value $\beta$. Left panel: $\beta >{\beta }^{*}$. Right panel: $\beta <{\beta }^{*}$.

Figure 7. $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs. A, respectively for $\beta >{\beta }^{*}=1$ and for $\beta <{\beta }^{*}$. Here $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ is the intersection of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ from Figure 6.

When it is known that the nociceptor density has the parametric form of type 2 given in (34), we can tune the trial value $\beta$ down or up toward the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ depending on whether the calculated $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ increases or decreases with A.

6.3. Constraint Function Calculated Using the Uniform Density

When the type of parametric form of $\rho \left(y\right)/{\rho }_0$ is unknown, we like to design a method to identify the true underlying density type among a set of candidate density types based on test data measured in experiments. In the analysis of reflex time above, we used quantity Q defined in (36), with behaviors described in (37), to distinguish between type 1 and type 2 densities. Now we look into using constraint functions as a tool for that purpose. With no information on the density type, we use the uniform density as the trial density in calculating the constraint function.

Let $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ denote the constraint function based on the test data (which is generated using type 2 density (34) with true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ ) and calculated using framework (20) with the uniform trial density $\rho \left(y\right)/{\rho }_0\equiv 1$. We select the uniform density to make the calculation of constraint function parameter-free, based solely on the test data.

Note that the uniform density is not a member of type 2 family (34). Consequently, constraint function $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ is not a special case of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$. To study the behavior of $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$, we try to connect it to $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$, the true constraint function. Both $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ are constructed from the same test data generated with the true underlying density (34). Each of them is calculated in framework (20) using a different type of trial density, and has different features:

· $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$ is the constraint function calculated using the correct density type (34) with the true parameter value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$. Operationally, the calculation of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$ is not realistic unless we know ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$. The theoretical advantage of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$ is that it shares a common intersection $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ for all values of A.

· $T_{\text{act,uni}}\left(z_{c,\text{uni}}\mathrm{<mo>\; ;\; </mo>}A\right)$ is the constraint function calculated using the uniform trial density. Operationally, we can always calculate $T_{\text{act,uni}}\left(z_{c,\text{uni}}\mathrm{<mo>\; ;\; </mo>}A\right)$ from test data. However, the point $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ in general is not on the curve $T_{\text{act,uni}}\left(z_{c,\text{uni}}\mathrm{<mo>\; ;\; </mo>}A\right)$.

The two constraint functions above are unified in $\Phi \left(y\mathrm{<mo>\; ;\; </mo>}1/A\right)$ via $y_{\text{act}}$, as described in (4).

$T_{\text{act}}=\Phi \left(y_{\text{act}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ (38)

The mapping between $y_{\text{act}}$ and $z_c$, however, depends on the trial density. As a result, variable $z_c$ is related to $y_{\text{act}}$ differently in the two constraint functions. For clarity of the discussion, we use $z_{c\text{,uni}}$ to denote the variable in $T_{\text{act,uni}}\left(z_{c,\text{uni}}\mathrm{<mo>\; ;\; </mo>}A\right)$, and use $z_c$ for the variable in $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$. The difference lies in the trial density used in the formulation. For the true density, $z_c$ is related to $y_{\text{act}}$ in (5). For the uniform density, $z_{c\text{,uni}}=A\cdot y_{\text{act}}$. Expressing $y_{\text{act}}$ in terms of $z_c$ or $z_{c\text{,uni}}$, we obtain

$y_{\text{act}}=\{\begin{array}{ll}\frac{1}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{F}^{-1}\left(\frac{{\beta }^{\mathrm{<mo>\; *\; </mo>}}{z}_c}{A}\right)\hfill & \text{for}\rho \left(y\right)/{\rho }_0=\left({\beta }^{\mathrm{<mo>\; *\; </mo>}}y\right){\text{e}}^{1-{\beta }^{\mathrm{<mo>\; *\; </mo>}}y}\hfill \\ \frac{z_{c,\text{uni}}}{A}\hfill & \text{for}\rho \left(y\right)/{\rho }_0=1\hfill \end{array}$ (39)

Combining (39) and (38), we write out $T_{\text{act,uni}}\left(z_{c,\text{uni}}\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)=\Phi \left(\frac{1}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{F}^{-1}\left(\frac{{\beta }^{\mathrm{<mo>\; *\; </mo>}}{z}_c}{A}\right)\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$

$T_{\text{act,uni}}\left(z_{c,\text{uni}}\mathrm{<mo>\; ;\; </mo>}A\right)=\Phi \left(\frac{z_{c,\text{uni}}}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ (40)

Since $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$ is calculated using the true density, it satisfies

$T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}=\Phi \left(\frac{1}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{F}^{-1}\left(\frac{{\beta }^{\mathrm{<mo>\; *\; </mo>}}{z}_c^{\mathrm{<mo>\; *\; </mo>}}}{A}\right)\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\text{for}\text{all}A$

Taking the limit as $A\to \infty$ yields

$T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}=\underset{A\to \infty }{lim}\Phi \left(\frac{1}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{F}^{-1}\left(\frac{{\beta }^{\mathrm{<mo>\; *\; </mo>}}{z}_c^{\mathrm{<mo>\; *\; </mo>}}}{A}\right)\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=\Phi \left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)$

In (40), letting $A\to \infty$ and using $\Phi \left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$, we obtain

$T_{\text{act,uni}}\left(z_{c\text{,uni}}\mathrm{<mo>\; ;\; </mo>}\infty \right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$

Let $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ denote the intersection of the pair $T_{\text{act,uni}}\left(z_{c\text{,uni}}\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act,uni}}\left(z_{c\text{,uni}}\mathrm{<mo>\; ;\; </mo>}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. Both $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ and $z_c^{\left(I\right)}\left({\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)\equiv z_c^{\mathrm{<mo>\; *\; </mo>}}$ are calculated by first mapping $T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ to $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ and then mapping to $z_c^{\left(I\right)}$. Given the measured temperature profile at reflex, Equation (38) maps $T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ to a unique $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$, independent of the trial density. It follows that in (39), line 1 with $z_c^{\left(I\right)}\left({\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)\equiv z_c^{\mathrm{<mo>\; *\; </mo>}}$ and line 2 with $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ are both equal to $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$.

$\frac{1}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{F}^{-1}\left(\frac{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{A}{z}_c^{\mathrm{<mo>\; *\; </mo>}}\right)=y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}=\frac{z_{c\text{,uni}}^{\left(I\right)}\left(A\right)}{A}$

Therefore, $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ and $z_c^{\mathrm{<mo>\; *\; </mo>}}$ are related by

$z_{c,\text{uni}}^{\left(I\right)}\left(A\right)=\frac{A}{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{F}^{-1}\left(\frac{{\beta }^{\mathrm{<mo>\; *\; </mo>}}}{A}{z}_c^{\mathrm{<mo>\; *\; </mo>}}\right)$ (41)

For type 2 density (34), the expansion of $F^{-1}\left(u\right)$ is given in (10) with $a_1=f^{\prime }\left(0\right)=\text{e}$ and $a_2=f^{″}\left(0\right)=-2\text{e}$. Substituting the expansion of $F^{-1}\left(u\right)$ into (41) yields

$z_{c,\text{uni}}^{\left(I\right)}\left(A\right)=\frac{\sqrt{2A}}{\sqrt{\text{e}{\beta }^{\mathrm{<mo>\; *\; </mo>}}}}\sqrt{z_c^{\mathrm{<mo>\; *\; </mo>}}}+\frac{2}{3\text{e}}{z}_c^{\mathrm{<mo>\; *\; </mo>}}$ (42)

Figure 8 compares simulated $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ and $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$. The simulation results confirm the theoretical prediction. In particular, $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ varies linearly with $\sqrt{A}$. We fit function $c_0+c_1\sqrt{A}$ to the data of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs $\sqrt{A}$. The fitting coefficients give us

$z_c^{\mathrm{<mo>\; *\; </mo>}}=\frac{3\text{e}}{2}{c}_0\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}=\frac{3c_0}{c_1^2}$

$z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is the intersection of two constraint functions calculated from test data using the uniform trial density. Given the test data, the process of calculating $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is parameter-free. However, $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ depends on the test data, which is affected by the underlying true density. When the true density has the type 2 form given in (34), result (42) predicts that $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ increases linearly with $\sqrt{A}$, unbounded as $A\to \infty$. Result (42) for type 2 density is in sharp contrast

Figure 8. Simulated results of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ for type 2 density given in (34). Here $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is based on the data from the true density (with ${\beta }^{*}=1$ ) but is calculated using the uniform trial density. Left panel: $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs A. Right panel: $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. $\sqrt{A}$.

with the situation for type 1 density, described in (33) of Section 5, where $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ converges to the true value $z_c^{\mathrm{<mo>\; *\; </mo>}}$ as $A\to \infty$. This difference in the behavior of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A provides another mechanism of distinguishing between type 1 and type 2 densities, in addition to the method of examining quantity Q, described in (37).

7. Type 3 Nociceptor Density: $\rho \left(y\right)={\rho }_0\cdot \theta \left(y-y_0\right)$

We write type 3 density in the same parametric form as that of types 1 and 2.

$\frac{\rho \left(y\right)}{{\rho }_0}=f\left(\beta y\right),f\left(s\right)=\theta \left(s-1\right),\beta \equiv \frac{1}{y_0}$ (43)

where $\theta \left(s\right)$ is the Heaviside step function.

$\theta \left(s\right)=\{\begin{array}{ll}0,\hfill & s<0\hfill \\ 1,\hfill & s\ge 0\hfill \end{array}$ (44)

The graph of type 3 $f\left(s\right)$ is shown in Figure 17. Type 3 density (43) is different from types 1 and 2 in that the nociceptor density jumps from 0 to ${\rho }_0$ at depth $y_0$. Because of the discontinuity, the Taylor expansions of $F\left(s\right)$ and $F^{-1}\left(u\right)$ derived in Section 3 are no longer valid. We write out $F\left(s\right)$ and $F^{-1}\left(u\right)$ directly

$F\left(s\right)={\displaystyle {\int }_0^s}\theta \left(s-1\right)\text{d}s=\{\begin{array}{ll}0,\hfill & s<1\hfill \\ s-1,\hfill & s\ge 1\hfill \end{array}$

$F^{-1}\left(u\right)=u+1\text{for}u>0$

The mapping between $y_{\text{act}}$ and $z_c$ is described by (7) and has the specific expression

${y_{\text{act}}|}_{\left(\text{reflex}\right)}=\frac{1}{\beta }{F}^{-1}\left(\frac{\beta z_c}{A}\right)=\frac{z_c}{A}+y_0\mathrm{,}\beta \equiv \frac{1}{y_0}$ (45)

7.1. Reflex Time

With the expression of $y_{\text{act}}$ in (45), Equation (14) for $t_{\text{ref,nd}}$ becomes

$\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}=H\left(y_{\text{act,nd}}\mathrm{,}{t}_{\text{ref,nd}}\right)$

$y_{\text{act,nd}}\equiv \mu y_{\text{act}}=\frac{1}{A_{\text{nd}}}+\mu y_0\mathrm{,}{A}_{\text{nd}}\equiv \frac{\mu z_c}{A}$

In the limit of $A\to \infty$, we have $y_{\text{act,nd}}\to \mu y_0$ and the equation for $t_{\text{ref,nd}}$ converges to

$\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}=H\left(\mu y_0\mathrm{,}{t_{\text{ref,nd}}|}_{A\to \infty }\right)$

Consider $h\left(t\mathrm{<mo>\; |\; </mo>}q\right)\equiv H\left(q\mathrm{,}t\right)$ as a function of t with q as a parameter. Let $h^{-1}\left(u\mathrm{<mo>\; |\; </mo>}q\right)$ be the inverse of $h\left(t\mathrm{<mo>\; |\; </mo>}q\right)$ and ${\tau }_0\left(q\right)\equiv h^{-1}\left(\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}|q\right)$. As $A\to \infty$, we have

${t_{\text{ref,nd}}|}_{A\to \infty }={\tau }_0\left(\mu y_0\right)$

By definition, ${\tau }_0\left(q\right)$ satisfies $H\left(q\mathrm{,}{\tau }_0\left(q\right)\right)=\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}$. We expand function $H\left(y_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$ around $\left(\mu y_0\mathrm{,}{\tau }_0\left(\mu y_0\right)\right)$ and write the equation for $t_{\text{ref,nd}}$ as

$0=\frac{\partial H\left(\mu y_0\mathrm{,}{\tau }_0\left(\mu y_0\right)\right)}{\partial t}\left(t_{\text{ref,nd}}-{\tau }_0\left(\mu y_0\right)\right)+\frac{\partial H\left(\mu y_0\mathrm{,}{\tau }_0\left(\mu y_0\right)\right)}{\partial y}\left(\frac{1}{A_{\text{nd}}}\right)$ (46)

Substituting $t_{\text{ref,nd}}\left(A_{\text{nd}}\right)={\tau }_0\left(\mu y_0\right)\left[1+c_A{\left(\frac{1}{A_{\text{nd}}}\right)}^{\alpha }+\cdots \right]$ into (46) to calculate $\alpha$ and $c_A$, and then mapping back to the physical quantities, we obtain

$\begin{array}{l}{t}_{\text{ref}}\left(A\right)=\frac{{\tau }_0\left(\mu y_0\right){\rho }_{\text{m}}{C}_p}{K{\mu }^2}\left[1+c_A\left(\frac{\mu z_c}{A}\right)+\cdots \right]\mathrm{,}\\ c_A=\frac{\frac{-\partial H\left(\mu y_0\mathrm{,}{\tau }_0\left(\mu y_0\right)\right)}{\partial y}}{{\tau }_0\left(\mu y_0\right)\frac{\partial H\left(\mu y_0\mathrm{,}{\tau }_0\left(\mu y_0\right)\right)}{\partial t}}\end{array}$

As A increases, $t_{\text{ref}}\left(A\right)$ converges to its limit with the residual proportional to 1/A.

$t_{\text{ref}}\left(A\right)-t_{\text{ref}}\left(\infty \right)~\frac{1}{A}$

This behavior is similar to that of type 2 density. For both type 2 and type 3, the nociceptor density is zero at the skin surface. Figure 9 plots the relation between $t_{\text{ref}}$ and A in two ways. Left panel: $t_{\text{ref}}$ vs A. Right panel: $t_{\text{ref}}$ vs 1/A. In particular, the right panel confirms that $t_{\text{ref}}$ is linear with respect to 1/A for large A, as predicted in the analysis above.

Figure 9. The reflex time ( $t_{\text{ref}}$ ) as a function of beam spot area (A) for type 3 nociceptor density (43) with $y_0^{\mathrm{<mo>\; *\; </mo>}}=\frac{1}{2}$. Left panel: $t_{\text{ref}}$ vs A. Right panel: $t_{\text{ref}}$ vs. 1/A.

7.2. Constraint Function $T_{act}\left(z_c;y_0,A\right)$

For type 3 density, it is more sensible to use $y_0$ as the parameter since it has the clear physical meaning of the depth at which the nociceptor density jumps from 0 to ${\rho }_0$. In the unified parametric form in (43), the generic parameter is $\beta \equiv 1/y_0$. Mathematically, working with $\beta \equiv 1/y_0$ allows us to use general results of the unified parametric form obtained in previous sections. In the analysis below, we will go back and forth between $y_0$ and $\beta$.

We adopt the general convention of using $y_0^{\mathrm{<mo>\; *\; </mo>}}$ to denote the true value and $y_0$ to represent the variable. The general form of constraint function with trial value $\beta$ is given in (20). Using $\beta =1/y_0$ and $F^{-1}\left(u\right)=u+1$, we write it as

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)=\Phi \left(\frac{z_c}{A}+y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ (47)

In the limit of $A\to \infty$, $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ converges to a horizontal line

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A=\infty \right)=\Phi \left(y_0\mathrm{<mo>\; ;\; </mo>0}\right)$

At $y_0^{\mathrm{<mo>\; *\; </mo>}}$, the constraint function shares the common intersection $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ for all A.

$T_{\text{act}}\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}{y}_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)=\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\text{for}\text{all}A$

In particular, we have $\Phi \left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>0}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. We expand $\Phi \left(\frac{z_c}{A}+y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ around $\left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ and expand $\Phi \left(y_0\mathrm{<mo>\; ;\; </mo>0}\right)$ around $\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,0}\right)$ to write them respectively as

$\Phi \left(\frac{z_c}{A}+y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}+\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\left(\frac{z_c-z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$ (48)

$\Phi \left(y_0\mathrm{<mo>\; ;\; </mo>0}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}+\frac{\partial }{\partial y}\Phi \left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>0}\right)\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$ (49)

Let $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ denote the $z_c$ -coordinate of the intersection of $T_{\text{act}}\left(z_c;y_0,A\right)$ and $T_{\text{act}}\left(z_c;y_0,\infty \right)$. $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ is governed by equating the RHSs of (48) and (49).

$\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\left(\frac{z_c-z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)=\frac{\partial }{\partial y}\Phi \left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>0}\right)\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$

Solving for $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ yields

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}+A\left(\frac{\frac{\partial }{\partial y}\Phi \left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>0}\right)}{\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)}-1\right)\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$ (50)

We expand $\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ with respect to $\frac{1}{A}$ into the form

$\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=\frac{\partial }{\partial y}\Phi \left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>0}\right)\left(1+\frac{b_1}{A}+\frac{b_2}{A^2}+\cdots \right)$ (51)

Substituting the expansion into (50), we obtain

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}+\left(-b_1+\left(b_1^2-b_2\right)\frac{1}{A}\right)\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$ (52)

In Appendix B, we show that the coefficients satisfy

$b_1>0\text{and}\left(b_1^2-b_2\right)>0$ (53)

At $y_0=y_0^{*}$, the intersection $z_c^{\left(I\right)}\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}$ is independent of beam spot area A. When $y_0\ne y_0^{\mathrm{<mo>\; *\; </mo>}}$, the trend of $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs A implies whether $y_0<y_0^{*}$ or $y_0>y_0^{*}$.

· For $y_0<y_0^{*}$, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ ascends toward $z_c^{*}-b_1\left(y_0-y_0^{*}\right)>z_c^{*}$ as A increases.

· For $y_0>y_0^{*}$, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ descends toward $z_c^{*}-b_1\left(y_0-y_0^{*}\right)<z_c^{*}$ as A increases.

In terms of $\beta \equiv 1/y_0$, the increase/decrease trend of $z_c^{\left(I\right)}$ vs A for type 3 density (43) resembles that for type 2 density (34). Both types of densities share the common feature that the nociceptor density is zero at skin surface: $\rho \left(0\right)=0$.

Figure 10 depicts simulated $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ for several values of A, respectively for $y_0<y_0^{*}$ and for $y_0>y_0^{*}$. Here constraint function $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ is based on test data (which is generated with true value $y_0^{*}=0.5$ ) and is calculated using formulation (47) with trial value $y_0$. The trend of $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs A is shown in Figure 11. The simulation results in Figure 10 and Figure 11 confirm the theoretically predicted trend above.

When it is known that the nociceptor density has the parametric form of type 3 given in (43), we can tune the trial value $y_0$ up or down toward the true value

Figure 10. Constraint curves for type 3 density (43). $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ is based on test data generated with true value $y_0^{*}=0.5$, and calculated using (47) with trial value $y_0$. Left panel: $y_0<y_0^{*}$. Right panel: $y_0>y_0^{*}$. Notice that for type 3 density, $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}\infty \right)$ is a horizontal line, whose height varies with $y_0$.

Figure 11. $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs A, respectively for $y_0<y_0^{*}=0.5$ and for $y_0>y_0^{*}$. Here $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ is the intersection of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}\infty \right)$ from Figure 10.

$y_0^{\mathrm{<mo>\; *\; </mo>}}$ depending on whether the calculated $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ increases or decreases with A.

7.3. Constraint Function Calculated Using the Uniform Density

When the type of parametric form of $\rho \left(y\right)/{\rho }_0$ is unknown, we apply the uniform trial density in calculating the constraint function and use it as a tool for probing the density type. This approach has the advantage of being operationally practical. Once the test data is available, the calculation of constraint function does not require any input parameters.

Let $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ denote the constraint function based on the test data (which is generated using type 3 density (43) with true value $y_0^{\mathrm{<mo>\; *\; </mo>}}$ ) and calculated using framework (47) with the uniform trial density $\rho \left(y\right)/{\rho }_0\equiv 1$. For type 3 parametric family (43), the uniform density is a special case at $y_0=0$. However, expansions (48) and (49) are only for the case of $y_0$ near $y_0^{\mathrm{<mo>\; *\; </mo>}}$, away from the skin surface. At $y_0=0$, the expansions will be different because of the insulated boundary condition $\frac{\partial }{\partial y}\Phi \left(0,t\right)=0$. At $y_0=0$, we have

$T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)=\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$

$T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\infty \right)=\Phi \left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)$

We expand $\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ around $\left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)$. It follows that

$\begin{array}{c}\Phi \left(\frac{z_c}{A};\frac{1}{A}\right)=\Phi \left(0;0\right)+\frac{\partial }{\partial v}\Phi \left(0;0\right)\frac{1}{A}+\frac{1}{2}\frac{{\partial }^2}{\partial y^2}\Phi \left(0;0\right)\frac{z_c^2}{A^2}\\ +\frac{1}{2}\frac{{\partial }^2}{\partial v^2}\Phi \left(0;0\right)\frac{1}{A^2}+\frac{{\partial }^2}{\partial y\partial v}\Phi \left(0;0\right)\frac{z_c}{A^2}+\cdots \end{array}$

Let $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ denote the intersection of $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\infty \right)$. It satisfies

${\Phi }_{yy}\left(z_{c,\text{uni}}^{\left(I\right)}\right){\left(A\right)}^2+2{\Phi }_{yv}{z}_{c,\text{uni}}^{\left(I\right)}\left(A\right)+2{\Phi }_{v}A+{\Phi }_{vv}+\cdots =0$

All derivatives are evaluated at $\left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)$. Solving for $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$, we obtain

$z_{c\text{,uni}}^{\left(I\right)}\left(A\right)=\sqrt{\frac{{\Phi }_v}{-{\Phi }_{yy}}}\cdot \sqrt{A}+\frac{2{\Phi }_{yv}}{-{\Phi }_{yy}}+\cdots$ (54)

Result (54) predicts that when the true underlying density is type 3 given in (43), $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ increases linearly with $\sqrt{A}$, unbounded as $A\to \infty$. This behavior is similar to that for type 2 density. Both type 3 and type 2 share the common feature of $\rho \left(0\right)=0$. Figure 12 compares simulated $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ and $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$. The simulation results confirm the theoretical prediction.

8. Type 4 Nociceptor Density: $\rho \left(y\right)={\rho }_0\left(1-\theta \left(y-y_0\right)/2\right)$

We represent type 4 density in the same parametric form as that of types 1 - 3.

$\frac{\rho \left(y\right)}{{\rho }_0}=f\left(\beta y\right)\mathrm{,}f\left(s\right)\mathrm{<mo>\; =\; </mo>1}-\frac{1}{2}\theta \left(s-1\right)\mathrm{,}\beta \equiv \frac{1}{y_0}$ (55)

where $\theta \left(s\right)$ is the Heaviside step function defined in (44). The graph of type 4 $f\left(s\right)$ is shown in Figure 17. Type 4 relative density (55) has value 1 in $\left[\mathrm{0,}{y}_0\right]$, and drops down to value 0.5 for $y>y_0$. When beam spot area A is sufficiently large, only a small depth of the skin needs to reach the activation temperature in order to trigger the withdrawal reflex. Thus, for large A, the behaviors of type 4 density are the same as those of the uniform density. For the purpose of revealing the effect of density jump at $y=y_0$, we examine the reflex time and the temperature profile at reflex in an intermediate range of A corresponding to the situation where the activated depth is around the density jump ( $y=y_0$ ).

Figure 12. Simulated results of $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ for type 3 density given in (43). Here $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ is based on the data from the true density (with $y_0^{*}=0.5$ ) but is calculated using the uniform trial density. Left panel: $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ vs. A. Right panel: $z_{c\text{,uni}}^{\left(I\right)}\left(A\right)$ vs $\sqrt{A}$.

8.1. Reflex Time

Because of the discontinuity in density profile (8), the Taylor expansions of $F\left(s\right)$ and $F^{-1}\left(u\right)$ derived in Section 3 are no longer valid. We write out $F\left(s\right)$ and $F^{-1}\left(u\right)$ directly.

$F\left(s\right)={\displaystyle {\int }_0^s}\left(1-\frac{1}{2}\theta \left(s-1\right)\right)\text{d}s=\{\begin{array}{ll}s,\hfill & s<1\hfill \\ \frac{1}{2}\left(s+1\right),\hfill & s\ge 1\hfill \end{array}$

$F^{-1}\left(u\right)=\{\begin{array}{ll}u,\hfill & u<1\hfill \\ 2u-1,\hfill & u\ge 1\hfill \end{array}$ (56)

The activated depth at reflex given in (7) has the expression

${y_{\text{act}}|}_{\left(\text{reflex}\right)}\left(A\right)=y_0{F}^{-1}\left(\frac{z_c}{y_{0}A}\right)=\{\begin{array}{ll}\frac{z_c}{A},\hfill & \frac{z_c}{A}<y_0\hfill \\ 2\frac{z_c}{A}-y_0,\hfill & \frac{z_c}{A}\ge y_0\hfill \end{array}$ (57)

The non-dimensional reflex time $t_{\text{ref,nd}}$ is related to A via $y_{\text{act}}\left(A\right)$ in Equation (14).

$\frac{\left(T_{\text{act}}-T_0\right)K\mu }{P_{\text{dep}}}=H\left(\mu y_{\text{act}}\left(A\right)\mathrm{,}{t}_{\text{ref,nd}}\right)$ (58)

Conventionally, we view $t_{\text{ref,nd}}$ as a function of A since in experiments the beam spot area is prescribed in test design and the corresponding reflex time is observed. To facilitate the mathematically analysis, we view A as a function of $t_{\text{ref,nd}}$. We study function $A\left(t_{\text{ref,nd}}\right)$ for type 4 density and connect it to its counterpart for the uniform density. Let

· $A^{\left(\text{tp4}\right)}\left(t_{\text{ref,nd}}\right)$ : function A vs. $t_{\text{ref,nd}}$ for type 4 density

· $A^{\left(\text{uni}\right)}\left(t_{\text{ref,nd}}\right)$ : function A vs. $t_{\text{ref,nd}}$ for the uniform density

· $t_{\text{ref,nd}}^{\left(\text{tp4}\right)}\left(1/A\right)$ : function $t_{\text{ref,nd}}$ vs. 1/A for type 4 density

· $t_{\text{ref,nd}}^{\left(\text{uni}\right)}\left(1/A\right)$ : function $t_{\text{ref,nd}}$ vs. 1/A for the uniform density.

For any given value of $t_{\text{ref,nd}}$, the corresponding $y_{\text{act}}\left(t_{\text{ref,nd}}\right)$ is completely determined by Equation (58), independent of the nociceptor density. We use $y_{\text{act}}\left(t_{\text{ref,nd}}\right)$ to connect $A^{\left(\text{tp4}\right)}\left(t_{\text{ref,nd}}\right)$ and $A^{\left(\text{uni}\right)}\left(t_{\text{ref,nd}}\right)$. For type 4 density, $y_{\text{act}}\left(A^{\left(\text{tp4}\right)}\right)$ is given in (57). Its inverse function is

$A^{\left(\text{tp4}\right)}\left(y_{\text{act}}\right)=\frac{z_c}{y_0}\cdot \frac{1}{F\left(\frac{y_{\text{act}}}{y_0}\right)}=\{\begin{array}{ll}\frac{z_c}{y_{\text{act}}}\mathrm{,}\hfill & y_{\text{act}}<y_0\hfill \\ \frac{2z_c}{y_{\text{act}}+y_0}\mathrm{,}\hfill & y_{\text{act}}\ge y_0\hfill \end{array}$ (59)

For the uniform density, $f\left(s\right)=1$, $F\left(s\right)=s$, and $A^{\left(\text{uni}\right)}\left(y_{\text{act}}\right)$ is given by

$A^{\left(\text{uni}\right)}\left(y_{\text{act}}\right)=\frac{z_c}{y_0}\cdot \frac{1}{F\left(\frac{y_{\text{act}}}{y_0}\right)}=\frac{z_c}{y_{\text{act}}}$ (60)

Combining (59) and (60), we can express $1/A^{\left(\text{uni}\right)}$ in terms of $1/A^{\left(\text{tp4}\right)}$.

$\frac{1}{A^{\left(\text{uni}\right)}}=\{\begin{array}{ll}\frac{1}{A^{\left(\text{tp4}\right)}}\mathrm{,}\hfill & \frac{1}{A^{\left(\text{tp4}\right)}}<\frac{y_0}{z_c}\hfill \\ \frac{2}{A^{\left(\text{tp4}\right)}}-\frac{y_0}{z_c}\mathrm{,}\hfill & \frac{1}{A^{\left(\text{tp4}\right)}}\ge \frac{y_0}{z_c}\hfill \end{array}$ (61)

(61) described the relation between $1/A^{\left(\text{tp4}\right)}\left(t_{\text{ref,nd}}\right)$ and $1/A^{\left(\text{uni}\right)}\left(t_{\text{ref,nd}}\right)$. Now we treat $t_{\text{ref,nd}}$ as the dependent variable and view it as a function of 1/A. (61) leads to

$t_{\text{ref,nd}}^{\left(\text{tp4}\right)}\left(1/A\right)=\{\begin{array}{ll}{t}_{\text{ref,nd}}^{\left(\text{uni}\right)}\left(1/A\right)\mathrm{,}\hfill & \frac{1}{A}<\frac{y_0}{z_c}\hfill \\ t_{\text{ref,nd}}^{\left(\text{uni}\right)}\left(\frac{2}{A}-\frac{y_0}{z_c}\right)\mathrm{,}\hfill & \frac{1}{A}\ge \frac{y_0}{z_c}\hfill \end{array}$ (62)

Equation (62) reveals that $t_{\text{ref,nd}}^{\left(\text{tp4}\right)}\left(1/A\right)$ is obtained from $t_{\text{ref,nd}}^{\left(\text{uni}\right)}\left(1/A\right)$ by a piecewise linear scaling on the independent variable 1/A. Function $t_{\text{ref,nd}}^{\left(\text{uni}\right)}\left(1/A\right)$ is smooth. After the piecewise linear scaling, $t_{\text{ref,nd}}^{\left(\text{tp4}\right)}\left(1/A\right)$ has a discontinuity in derivative. Figure 13 plots the relation between $t_{\text{ref}}^{\left(\text{tp4}\right)}$ and A in two ways. Left panel: $t_{\text{ref}}$ vs. 1/A. Right panel: derivative of $t_{\text{ref}}$ vs. 1/A. In particular, the right panel verifies that $t_{\text{ref}}$ vs. 1/A has a discontinuity in derivative, as predicted in the analysis above.

8.2. Constraint Function $T_{act}\left(z_c;y_0,A\right)$

For type 4 density, it is more sensible to choose $y_0$ as the parameter since it has the clear physical meaning of the depth at which the nociceptor density drops sharply from ${\rho }_0$ to ${\rho }_0/2$. In the unified parametric form in (55), the generic parameter is $\beta \equiv 1/y_0$. We adopt the general convention of using $y_0^{\mathrm{<mo>\; *\; </mo>}}$ to denote the true value and $y_0$ to represent the variable. The general form of constraint function with trial value $\beta$ is given in (20). Using $\beta =1/y_0$ and $F^{-1}\left(u\right)$ given in (56), we write it as

Figure 13. The relation between reflex time ( $t_{\text{ref}}$ ) and beam spot area (A) for type 4 density (55) with $y_0^{*}=0.25$. Left panel: $t_{\text{ref}}$ vs 1/A. Right panel: derivative of $t_{\text{ref}}$ vs. 1/A.

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)=\{\begin{array}{ll}\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\mathrm{,}\hfill & \frac{z_c}{A}<y_0\hfill \\ \Phi \left(2\frac{z_c}{A}-y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\mathrm{,}\hfill & \frac{z_c}{A}\ge y_0\hfill \end{array}$ (63)

when beam spot area A is sufficiently large to make $\frac{z_c^{*}}{A}<y_0^{*}$, (57) gives $y_{\text{act}}<y_0^{*}$ and the measured data for type 4 density is the same as that for the uniform density. Constraint function (63) is based on the data for true value $y_0^{\mathrm{<mo>\; *\; </mo>}}$ and is calculated with trial value $y_0$. For any positive trial value $y_0>0$ and a finite interval of $z_c$ near $z_c^{\mathrm{<mo>\; *\; </mo>}}$ in consideration, when A is sufficiently large, we have $\frac{\text{Interval}\text{of}{z}_c}{A}<min\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_0\right)$ and

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)=\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\text{for}\text{large}A$

Here $\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ is the constraint function based on the data for the uniform density and is calculated using the uniform trial density. $\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ is independent of trial value $y_0$, and passes through $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ for all values of A.

$\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\text{for}\text{all}A$

In the limit of $1/A\to 0$, we have $\Phi \left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. Let $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ denote the $z_c$ -coordinate of the intersection of $T_{\text{act}}\left(z_c;y_0,A\right)$ and $T_{\text{act}}\left(z_c;y_0,\infty \right)$. Our analysis above shows that

· For any $y_0$, $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)=\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ when A is sufficiently large.

· For any $y_0$, $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A=\infty \right)=\Phi \left(\mathrm{0\; <mo>\; ;\; </mo>0}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$.

· For any $y_0$, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}$ when A is sufficiently large.

More specifically, to have $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}$, we only need $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)=\Phi \left(\frac{z_c}{A}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ for $z_c$ near $z_c^{\mathrm{<mo>\; *\; </mo>}}$. The condition $\frac{\text{Interval}\text{of}{z}_c}{A}<min\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_0\right)$ becomes $A>\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{min\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_0\right)}$. We conclude that

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}\text{for}A>\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{min\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_0\right)}$ (64)

Thus, for large A, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ is independent of $y_0$. To probe the position of trial value $y_0$ relative to true value $y_0^{\mathrm{<mo>\; *\; </mo>}}$, we need to examine the behavior of $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ for $A<\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{max\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_0\right)}$. In this range of A, constraint function (63) takes the form

$T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)=\Phi \left(2\frac{z_c}{A}-y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$

At true value $y_0^{\mathrm{<mo>\; *\; </mo>}}$, $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$ shares the common intersection $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ for all A.

$\Phi \left(2\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}-y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$

We expand $\Phi \left(2\frac{z_c}{A}-y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ around $\left(\frac{2z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}-y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\frac{1}{A}\right)$.

$\Phi \left(2\frac{z_c}{A}-y_0\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}+\frac{\partial }{\partial y}\Phi \left(2\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}-y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)\left(2\frac{z_c-z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}-\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)\right)$

Setting the LHS to $T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ and solving for $z_c$, we obtain $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}+\frac{A}{2}\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)\text{for}A<\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{max\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}_0\right)}$ (65)

For $y_0<y_0^{*}$ and $\frac{z_c^{*}}{y_0^{*}}<A<\frac{z_c^{*}}{y_0}$, the activated depth at reflex is $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}=\frac{z_c^{*}}{A}>y_0$. Substituting the expression of $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ into (57) and solving for $z_c$, we have

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}-\frac{1}{2}\left(z_c^{\mathrm{<mo>\; *\; </mo>}}-A{y}_0\right)\text{for}\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0^{\mathrm{<mo>\; *\; </mo>}}}<A<\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0}$ (66)

For $y_0>y_0^{*}$ and $\frac{z_c^{*}}{y_0}<A<\frac{z_c^{*}}{y_0^{*}}$, the activated depth at reflex is $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}=2\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}-y_0^{\mathrm{<mo>\; *\; </mo>}}$. This $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ may or may not be beyond depth $y_0$. The case of $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}>y_0$ corresponds to $A<\frac{2z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0+y_0^{\mathrm{<mo>\; *\; </mo>}}}$. Substituting the expression of $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ into (57) and solving for $z_c$ yields

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}+\frac{A}{2}\left(y_0-y_0^{\mathrm{<mo>\; *\; </mo>}}\right)\text{for}\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0}<A<\frac{2z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0+y_0^{\mathrm{<mo>\; *\; </mo>}}}$ (67)

Similarly, the case of $y_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}<y_0$ corresponds to $A>\frac{2z_c^{*}}{y_0+y_0^{*}}$. We get

$z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}+\left(z_c^{\mathrm{<mo>\; *\; </mo>}}-A{y}_0^{\mathrm{<mo>\; *\; </mo>}}\right)\text{for}\frac{2z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0+y_0^{\mathrm{<mo>\; *\; </mo>}}}<A<\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0^{\mathrm{<mo>\; *\; </mo>}}}$ (68)

Summarizing results for various cases, we write out a complete description of $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$.

The case of $y_0<y_0^{*}$,

$z_c^{\left(I\right)}\left(y_0,A\right)=\{\begin{array}{ll}{z}_c^{*}+\frac{A}{2}\left(y_0-y_0^{*}\right)\hfill & \text{for}A<\frac{z_c^{*}}{y_0^{*}}\hfill \\ z_c^{*}-\frac{1}{2}\left(z_c^{*}-A{y}_0\right)\hfill & \text{for}\frac{z_c^{*}}{y_0^{*}}<A<\frac{z_c^{*}}{y_0}\hfill \\ z_c^{*}\hfill & \text{for}A>\frac{z_c^{*}}{y_0}\hfill \end{array}$

The case of $y_0>y_0^{*}$,

$z_c^{\left(I\right)}\left(y_0,A\right)=\{\begin{array}{ll}{z}_c^{*}+\frac{A}{2}\left(y_0-y_0^{*}\right)\hfill & \text{for}A<\frac{2z_c^{*}}{y_0^{*}+y_0}\hfill \\ z_c^{*}+\left(z_c^{*}-A{y}_0^{*}\right)\hfill & \text{for}\frac{2z_c^{*}}{y_0^{*}+y_0}<A<\frac{z_c^{*}}{y_0^{*}}\hfill \\ z_c^{*}\hfill & \text{for}A>\frac{z_c^{*}}{y_0^{*}}\hfill \end{array}$ (69)

At $y_0=y_0^{\mathrm{<mo>\; *\; </mo>}}$, we have $z_c^{\left(I\right)}\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)=z_c^{\mathrm{<mo>\; *\; </mo>}}$ for all A. When $y_0\ne y_0^{\mathrm{<mo>\; *\; </mo>}}$, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ has non-trivial dependence on A. The trend of $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs A tells us whether $y_0<y_0^{*}$ or $y_0>y_0^{*}$.

· When $y_0<y_0^{*}$, from small A to large A, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ starts below $z_c^{\mathrm{<mo>\; *\; </mo>}}$, decreases further below $z_c^{\mathrm{<mo>\; *\; </mo>}}$, and then reverts rapidly back to $z_c^{\mathrm{<mo>\; *\; </mo>}}$ and stays there.

· When $y_0>y_0^{*}$, from small A to large A, $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ starts above $z_c^{\mathrm{<mo>\; *\; </mo>}}$, increases further above $z_c^{\mathrm{<mo>\; *\; </mo>}}$, and then reverts rapidly back to $z_c^{\mathrm{<mo>\; *\; </mo>}}$ and remains there.

Figure 14 shows simulated $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ for several values of A, respectively for $y_0<y_0^{*}$ and for $y_0>y_0^{*}$. Here constraint function $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ is based on test data (which is generated with true value $y_0^{*}=0.25$ ) and is calculated using formulation (63) with trial value $y_0$. The trend of $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs A is shown in Figure 15. The simulation results in Figure 14 and Figure 15 confirm the theoretically predicted trend above.

When it is known that the nociceptor density has the parametric form of type 4 given in (55), we can tune the trial value $y_0$ up or down toward the true value $y_0^{\mathrm{<mo>\; *\; </mo>}}$ depending on the behavior of the calculated $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs. A.

8.3. Constraint Function Calculated Using the Uniform Density

When the type of parametric form of $\rho \left(y\right)/{\rho }_0$ is unknown, we use the uniform

Figure 14. Constraint curves for type 4 density (55). $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ is based on test data generated with true value $y_0^{*}=0.25$, and calculated using (63) with trial value $y_0$. Left panel: $y_0<y_0^{*}$. Right panel: $y_0>y_0^{*}$.

Figure 15. $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ vs A, respectively for $y_0<y_0^{*}=0.25$ and for $y_0>y_0^{*}$. Here $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ is the intersection of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}\infty \right)$ from Figure 14.

trial density in calculating the constraint function. The goal is to use the resulting constraint function as a tool to probe the underlying unknown density type.

Let $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ denote the constraint function based on the test data (which is generated using type 4 density (55) with true value $y_0^{\mathrm{<mo>\; *\; </mo>}}$ ) and calculated using framework (63) with the uniform trial density $\rho \left(y\right)/{\rho }_0\equiv 1$. For type 4 parametric family (55), the uniform density is a special case with large $y_0$.

It follows that the behavior of $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ is the same as that of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{y}_0\mathrm{,}A\right)$ for large $y_0$, which we analyzed in the previous subsection.

Let $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ denote the intersection of the pair $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\infty \right)\equiv T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$. Intersection $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ has the same behavior as $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ for large $y_0$. Based on result (69) for $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$ in the previous subsection, we conclude for $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ that

$z_{c,\text{uni}}^{\left(I\right)}\left(A\right)=\{\begin{array}{ll}{z}_c^{\mathrm{<mo>\; *\; </mo>}}+\left(z_c^{\mathrm{<mo>\; *\; </mo>}}-A{y}_0^{\mathrm{<mo>\; *\; </mo>}}\right)\hfill & \text{for}A<\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0^{\mathrm{<mo>\; *\; </mo>}}}\hfill \\ z_c^{\mathrm{<mo>\; *\; </mo>}}\hfill & \text{for}A>\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0^{\mathrm{<mo>\; *\; </mo>}}}\hfill \end{array}$ (70)

Result (70) predicts that when the true underlying density is type 4 given in (55), $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ starts at slightly below $2z_c^{\mathrm{<mo>\; *\; </mo>}}$ for small A; decreases linearly toward $z_c^{\mathrm{<mo>\; *\; </mo>}}$ as A increases; arrives at $z_c^{\mathrm{<mo>\; *\; </mo>}}$ at $A=\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{y_0^{\mathrm{<mo>\; *\; </mo>}}}$ and stays there for $A>\frac{z_c^{*}}{y_0^{*}}$. The qualitative behavior of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ converging to $z_c^{\mathrm{<mo>\; *\; </mo>}}$ for large A is similar for type 1 and type 4 densities. Figure 16 compares simulated $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ and $z_c^{\left(I\right)}\left(y_0\mathrm{,}A\right)$. The simulation results confirm the theoretical prediction.

Figure 16. Simulated results of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ for type 4 density given in (55). Here $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ is based on the data from the true density (with $y_0^{*}=0.25$ ) but is calculated using the uniform trial density.

9. Summary

In this study, we investigate the nociceptor density of the form: $\rho \left(y\right)={\rho }_{0}f\left(\beta y\right)$, and its effect on heat-induced withdrawal reflex. We examine 4 types of $f\left(s\right)$ illustrated in Figure 17. Each $f\left(s\right)$ shown has a characteristic length 1. When scaled in the depth direction by parameter $\beta$, each $f\left(s\right)$ yields a family of density profiles.

We consider the situation where the reflex time and the temperature profile at the reflex are measurable in tests. We build the mathematical formulation for extracting 3 key parameters from test data:

· the activation temperature $T_{\text{act}}$ for heat-sensitive nociceptors,

· the critical threshold $z_c$ on the equivalent activated volume, and

· the parameter $\beta$ in the relative density.

Our general strategy is to identify distinct behaviors for different densitytypes and distinct behaviors for different regions of parameter values. We compare these theoretical patterns with the observed patterns from test data to pinpoint the underlying unknown density type. We inspect the behavior calculated using a trial parameter value in the parametric form to determine whether the trial value is below or above the true value. Then we use that information to tune the trial value up or down accordingly toward the true value. The process is repeated with the new trial value until convergence. To best illustrate each key module individually in its own setting, we divide the task of finding the density type and the parameter value into stages and consider several problems. Mathematically, we proceed from the simplest problem to the realistic one in which both the density type and the parameter value are unknown. The solution of a simpler problem provides the building blocks for solving a more complicated problem.

Figure 17. Four types of nociceptor density vs. depth. Each type has been transformed into the standard form of maximum density 1 and characteristic depth 1.

Problem 1:

When both the density type $f\left(s\right)$ and the true scaling parameter ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ are known for the nociceptor density, the test data allow us to construct true constraint functions on $\left(z_c\mathrm{,}{T}_{\text{act}}\right)$. Multiple constraint functions, obtained from tests at different values of beam spot area A, share a common point at the true value of $\left(z_c\mathrm{,}{T}_{\text{act}}\right)$. Parameters $z_c$ and $T_{\text{act}}$ are determined by finding the intersection of these distinct constraint functions.

Problem 2:

When the density type $f\left(s\right)$ is given but the true scaling parameter ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ is unknown, we construct trial constraint functions using trial values of $\beta$. Let $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ denote the trial constraint function that is based on the measured data (which is generated with true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ ) and that is calculated using the given parametric form with trial value $\beta$. When the trial value $\beta$ is different from the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$, in general, the true value $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is not on trial constraint functions, and trial constraint functions at different values of beam spot area A do not share a common intersection. The behaviors of trial constraint functions $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ are demonstrated for the 4 density types, respectively in Figure 2, Figure 6, Figure 10, and Figure 14. We look at the intersection of a pair of trial constraint functions: $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}\infty \right)$. Let $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ denote the intersection of the pair, which is affected by the trial value $\beta$ used in calculating the trial constraint functions and by the beam spot area A used in tests. Given a trial value $\beta$, we examine the trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs A. For each of the 4 density types, as the beam spot area A increases, $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs A demonstrates distinct trend of increasing or decreasing, respectively when $\beta >{\beta }^{*}$ and when $\beta <{\beta }^{*}$. The trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs. A is illustrated for the 4 density types, respectively in Figure 3, Figure 7, Figure 11, and Figure 15. Depending on the given density type and the observed trend of $z_c^{\left(I\right)}\left(\beta \mathrm{,}A\right)$ vs. A, we tune the trial value $\beta$ accordingly to approach the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$. At the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$, $z_c^{\left(I\right)}\left({\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)\equiv z_c^{\mathrm{<mo>\; *\; </mo>}}$ is independent of A. Once we arrive at the true value ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$, the subsequent procedure of finding parameters $z_c$ and $T_{\text{act}}$ is the same as described in Problem 1 above: the intersection of $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}A\right)$ and $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\infty \right)$ gives us the true value of $\left(z_c\mathrm{,}{T}_{\text{act}}\right)$.

Problem 3:

When both the density type $f\left(s\right)$ and the true scaling parameter ${\beta }^{\mathrm{<mo>\; *\; </mo>}}$ are unknown, we construct trial constraint functions using the uniform trial density $\rho \left(y\right)/{\rho }_0\equiv 1$, which is parameter-free. The trial constraint function using the uniform density is denoted by $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$. When the true density is not the uniform density, the true value $\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{T}_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is not on trial constraint function $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$, and for different values of A trial constraint functions $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ do not share a common point. Similar to what we did on trial constraint functions $T_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)$ in Problem 2 (where the density type is given), here we look at the intersection of a pair of trial constraint functions: $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}A\right)$ and $T_{\text{act,uni}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\infty \right)$. Let $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ denote the intersection of the pair, which is affected only by the beam spot area A used in tests. Given a collection of measured data sets at several values of A, we examine the trend of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A. The trend behavior of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A is displayed for the 4 density types, respectively in Figure 4, Figure 8, Figure 12, and Figure 16. The behaviors of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A for types 1 and 4 are distinct from each other and are distinct from those for types 2 and 3. In addition to examining the trend of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A, we also check the convergence pattern of reflex time $t_{\text{ref}}$ as beam spot area A increases, which is shown for the 4 density types, respectively in Figure 1, Figure 5, Figure 9, and Figure 13. Again, the patterns of $t_{\text{ref}}$ vs. A for types 1 and 4 are distinct from each other and are distinct from those for types 2 and 3. Density types 2 and 3 have similar behaviors in both $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A and $t_{\text{ref}}$ vs. A. This is not completely surprising since both density types share the common feature of having zero nociceptor density at the skin surface. When we are presented with the task of identifying the underlying nociceptor density from among the 4 density types, it is still a challenge to distinguish between type 2 and type 3 based on the trend of $z_{c,\text{uni}}^{\left(I\right)}\left(A\right)$ vs. A or the trend of $t_{\text{ref}}$ vs. A. We need to explore other formulations and analytical tools for differentiating types 2 and 3. One possibility is to use alternative standardized parameter-free forms other than the uniform density as the trial density in calculating the trial constraint function. Another possibility is to use the trend of system behaviors vs. varying applied beam power, in addition to varying beam spot area. Once the density type is selected, the remaining task of determining the true values of parameters $\beta$ and $\left(z_c\mathrm{,}{T}_{\text{act}}\right)$ is the same as described in Problem 2 above.

The goal of our study is to develop a methodology that utilizes test data to 1) identify the density type for the underlying nociceptor density, 2) find by trial and error the true scaling parameter in the parametric form, and 3) determine the activation temperature of nociceptors and the critical threshold on the equivalent activated volume. The results of 1) and 2) basically specify the nociceptor density profile vs depth. We assume the test data available include measurements of the reflex time and of the spatial temperature profile at reflex. Combining the procedures outlined in Problems 1, 2 and 3 above, we obtain such a methodology exactly for this purpose.

Disclaimer

The authors thank Dr. John Biddle and Dr. Stacy Teng of Institute for Defense Analyses for introducing the ADT CHEETEH-E model and for many helpful discussions. The authors acknowledge the Joint Intermediate Force Capabilities Office 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.

Appendix A

Derivation of $li{m}_{A\to \infty }{T}_{\text{act}}\left(z_c\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}A\right)=T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$ in (25)

We show the convergence in three steps.

Step 1: We first look at the factor ${\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\frac{1}{A}$ in the expression of ${\frac{\partial T_{\text{act}}}{\partial z_c}|}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}$ in (22). As $A\to \infty$ we have $\epsilon \to 0$. We expand ${\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)$ as $\epsilon \to 0$. In both cases 1 and 2, we have $li{m}_{\epsilon \to 0}{\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\cdot \epsilon =0$. It follows that

$\underset{A\to \infty }{lim}{\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\frac{1}{A}=\underset{\epsilon \to 0}{lim}{\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\cdot \epsilon \frac{1}{\beta z_c}=0$ (71)

Step 2: Next we look at the factor $\frac{\partial }{\partial y}\Phi \left(\xi \mathrm{<mo>\; ;\; </mo>}1/A\right)$ in the expression of ${\frac{\partial T_{\text{act}}}{\partial z_c}|}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}$ in (22). With $li{m}_{A\to \infty }\xi =0$ and the insulated boundary condition $\frac{\partial T\left(\mathrm{0,}t\right)}{\partial y}=0$, we have

$\underset{A\to \infty }{lim}\frac{\partial }{\partial y}\Phi \left(\xi \mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)=\underset{\epsilon \to 0}{lim}\frac{\partial }{\partial y}T\left(\xi \mathrm{,}{t}_{\text{ref}}\left(A\right)\right)=0$ (72)

Taking the limit as $A\to \infty$ in (22), and using results (71) and (72), we conclude

$\underset{A\to \infty }{lim}\frac{\partial T_{\text{act}}}{\partial z_c}=\underset{A\to \infty }{lim}\frac{\partial }{\partial y}\Phi \left(\xi \mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right){\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)\frac{1}{A}=0$

Step 3: We show that the term ${\eta |}_{\left(z_c^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\beta }^{\mathrm{<mo>\; *\; </mo>}}\right)}$ in (24) is bounded as $A\to \infty$. In the expression of $\eta$ in (23), we look at $\frac{F^{-1}\left(\epsilon \right)}{\epsilon {\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)}$, the component that varies with $\epsilon$. Using expansions (8) and (9) in case 1 or expansions (10) and (11) in case 2, we get

$\frac{F^{-1}\left(\epsilon \right)}{\epsilon {\left(F^{-1}\right)}^{\prime }\left(\epsilon \right)}=\{\begin{array}{ll}1+\frac{a_1}{2a_0^2}\epsilon +\cdots \hfill & \text{for}\text{case}1\hfill \\ 2+\frac{\sqrt{2}{a}_2}{3a_1^{3/2}}{\epsilon }^{1/2}+\cdots \hfill & \text{for}\text{case}2\hfill \end{array}$ (73)

In both cases, we have $li{m}_{\epsilon \to 0}\eta =\text{finite}$.

Combining the results of three steps above, we conclude that for any trial value $\beta$, as $A\to \infty$, the constraint function converges to $T_{\text{act}}^{\mathrm{<mo>\; *\; </mo>}}$, as described in (25) in the main text.

Appendix B

Derivation of property (53)

Property (53) describes coefficients in the expansion of $\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)$ with respect to $1/A$. Recall that $\Phi \left(y\mathrm{<mo>\; ;\; </mo>}v\right)\equiv T\left(y\mathrm{,}{t}_{\text{ref}}\left(1/v\right)\right)$ where $v\equiv 1/A$. Function $T\left(y\mathrm{,}t\right)$ as given in (12) and (13), has the properties below

1) $\frac{\partial }{\partial y}T\left(y\mathrm{,}t\right)$ is always negative.

2) At any given depth y, the absolute value of $\frac{\partial }{\partial y}T\left(y\mathrm{,}t\right)$ increases with time t.

3) At any given time t, the absolute value of $\frac{\partial }{\partial y}T\left(y\mathrm{,}t\right)$ start with 0 value at $y=0$, increases with y until the inflection point at $y=\kappa \left(t\right)$, and then decreases to zeros beyond the inflection point. The depth of inflection point $y=\kappa \left(t\right)$ increases with t.

We consider the situation of $y_0^{*}<\kappa \left(t_0\right)$. That is, the nociceptor density jump occurs before the inflection point of temperature profile at reflex. Combining properties 1 and 3 of $T\left(y\mathrm{,}t\right)$ with $y_0^{*}<\kappa \left(t_0\right)$, we have

$\underset{\text{positive}}{\underset{︸}{-\frac{\partial }{\partial y}T\left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{t}_{\text{ref}}\left(A\right)\right)}}>\underset{\text{positive}}{\underset{︸}{-\frac{\partial }{\partial y}T\left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{t}_{\text{ref}}\left(A\right)\right)}}\text{for}\text{all}A$ (74)

As beam spot area A increases $t_{\text{ref}}\left(A\right)$ decreases and is bounded from below by $t_{\text{ref}}\left(A\right)\ge t_{\text{ref}}\left(A=\infty \right)={\tau }_0\left(\mu y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$. In turn, ${\tau }_0\left(\mu y_0^{\mathrm{<mo>\; *\; </mo>}}\right)$ is bounded by ${\tau }_0\left(\mu y_0^{*}\right)\ge {\tau }_0\left(0\right)=t_0>0$, reflecting that the presence of a skin layer with no nociceptor delays the occurrence of withdrawal reflex. Thus, we have

$t_{\text{ref}}\left(A\right)\ge t_{\text{ref}}\left(A=\infty \right)\ge t_0>0\text{independent}\text{of}A\text{and}{y}_0^{\mathrm{<mo>\; *\; </mo>}}$

Applying property 2 of $T\left(y\mathrm{,}t\right)$ to (74) with $t_{\text{ref}}\left(A\right)\ge t_{\text{ref}}\left(A=\infty \right)$, we obtain

$-\frac{\partial }{\partial y}\Phi \left(\frac{z_c^{\mathrm{<mo>\; *\; </mo>}}}{A}+y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>}\frac{1}{A}\right)>-\frac{\partial }{\partial y}\Phi \left(y_0^{\mathrm{<mo>\; *\; </mo>}}\mathrm{<mo>\; ;\; </mo>0}\right)>0$ (75)

The first part of (53), $b_1>0$, follows from inequality (75) and definition of $b_1$ in (51). The second part of (53), $\left(b_1^2-b_2\right)>0$, is verified in simulations.

Conflicts of Interest

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

References

[1]	Gandhi, O.P. and Riazi, A. (1986) Absorption of Millimeter Waves by Human Beings and Its Biological Implications. IEEE Transactions on Microwave Theory and Techniques, 34, 228-235. https://doi.org/10.1109/TMTT.1986.1133316
[2]	Romanenko, S., Begley, R., Harvey, A.R., Hool, L. and Wallace, V.P. (2017) The Interaction between Electromagnetic Fields at Megahertz, Gigahertz and Terahertz Frequencies with Cells, Tissues and Organisms: Risks and Potential. Journal of The Royal Society Interface, 14, Article ID: 20170585. https://doi.org/10.1098/rsif.2017.0585
[3]	Zhadobov, M., Chahat, N., Sauleau, R., Le Quement, C. and Le Drean, Y. (2011) Millimeter-Wave Interactions with the Human Body: State of Knowledge and Recent Advances. International Journal of Microwave and Wireless Technologies, 3, 237-247. https://doi.org/10.1017/S1759078711000122
[4]	Foster, K.R., Ziskin, M.C. and Balzano, Q. (2016) Thermal Response of Human Skin to Microwave Energy: A Critical Review. Health Physics, 111, 528-541. https://doi.org/10.1097/HP.0000000000000571
[5]	Ozen, S., Helhel, S. and Bilgin, S. (2011) Temperature and Burn Injury Prediction of Human Skin Exposed to Microwaves: A Model Analysis. Radiation and Environmental Biophysics, 50, 483-489. https://doi.org/10.1007/s00411-011-0364-y
[6]	Nelson, D.A., Nelson, M.T., Walters, T.J. and Mason, P.A. (2000) Skin Heating Effects of Millimeter-Wave Irradiation—Thermal Modeling Results. IEEE Transactions on Microwave Theory and Techniques, 48, 2111-2120. https://doi.org/10.1109/22.884202
[7]	Laakso, I., Morimoto, R., Heinonen, J., Jokela, K. and Hirata, A. (2017) Human Exposure to Pulsed Fields in the Frequency Range from 6 to 100 GHz. Physics in Medicine & Biology, 62, 6980-6992. https://doi.org/10.1088/1361-6560/aa81fe
[8]	Durney, C.H., Massoudi, H. and Iskander, M.F. (1986) Radiofrequency Radiation Dosimetry Handbook. 4th Edition, Brooks AFB, TX: USAF School Aerospace Med.
[9]	Xu, F. and Lu, T. (2011) Introduction to Skin Biothermomechanics and Thermal Pain. Springer, New York. https://doi.org/10.1007/978-3-642-13202-5
[10]	Gabriel, S., Lau, R.W. and Gabriel, C. (1996) The Dielectric Properties of Biological Tissues: II Measurements in the Frequency Range 10 Hz to 20 GHz. Physics in Medicine and Biology, 41, 2251-2269. https://doi.org/10.1088/0031-9155/41/11/002
[11]	Hwang, H., Yim, J., Cho, J.-W., Cheon, C. and Kwon, Y. (2003) 110 GHz Broadband Measurement of Permittivity on Human Epidermis Using 1 mm Coaxial Probe. International Microwave Symposium Digest, 1, 399-402.
[12]	Romanenko, S., Harvey, A.R., Hool, L., Fan, S. and Wallace, V.P. (2019) Millimeter Wave Radiation Activates Leech Nociceptors via TRPV1-Like Receptor Sensitization. Biophysical Journal, 116, 2331-2345. https://doi.org/10.1016/j.bpj.2019.04.021
[13]	Walters, T.J., Blick, D.W., Johnson, L.R., Adair, E.R. and Foster, K.R. (2000) Heating and Pain Sensation Produced in Human Skin by Millimeter Waves: Comparison to a Simple Thermal Model. Health Physics, 78, 259-267. https://doi.org/10.1097/00004032-200003000-00003
[14]	Plaghki, L., Decruynaere, C., Van Dooren, P. and Le Bars, D. (2010) The Fine Tuning of Pain Thresholds: A Sophisticated Double Alarm System. PLoS ONE, 5, e10269. https://doi.org/10.1371/journal.pone.0010269
[15]	Parker, J.E., Nelson, E.J., Beason, C.W. and Cook, M.C. (2017) Thermal and Behavioral Effects of Exposure to 30Kw, 95-GHz Millimeter Wave Energy. Technical Report, AFRL-RH-FS-TR-2017-0016.
[16]	Parker, J.E., Nelson, E.J., Beason, C.W. and Cook, M.C. (2017) Effects of Variable Spot Size on Human Exposure to 95-GHz Millimeter Wave Energy. Technical Report, AFRL-RH-FS-TQ-2017-0017.
[17]	Wang, H., 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