Applied Mathematics
Vol.11 No.08(2020), Article ID:102614,37 pages
10.4236/am.2020.118053
Effect of Depth-Dependent Nociceptor Density on the Heat-Induced Withdrawal Reflex
Hongyun Wang1, Wesley A. Burgei2, Hong Zhou3*
1Department of Applied Mathematics, University of California, Santa Cruz, CA, USA
2U.S. Department of Defense, Joint Intermediate Force Capabilities Office, Quantico, VA, USA
3Department of Applied Mathematics, Naval Postgraduate School, Monterey, CA, USA
Copyright © 2020 by author(s) and Scientific Research Publishing Inc.
This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).
http://creativecommons.org/licenses/by/4.0/
Received: June 17, 2020; Accepted: August 28, 2020; Published: August 31, 2020
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 ).
· : 2-D coordinates on the skin surface; is the 3-D coordinates.
· : 3-D spatial temperature profile of the skin.
· : activation temperature of nociceptors; given , the activation status of a nociceptor at is governed by the indicator function
· : nociceptor density at depth y (number per volume).
· : total number of activated nociceptors in the skin,
· : characteristic reference nociceptor density.
When the nociceptor density is uniform, , 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:
(1)
The advantages of this input dose are 1) it makes the dose quantity independent of the nociceptor density , and 2) it shifts the effect of into the dose threshold . In the dose response model, withdrawal reflex occurs when where the effects of and all other factors are reflected the single metric quantity .
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 based on reference density , and use as the input dose.
Input dose in the case of non-uniform density:
(2)
The dose response relation has the same form as in the case of uniform density.
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 defined in (2) requires only the relative density . The effect of is contained in the dose threshold . When is given, the dose response model has only two unknown parameters: and . In a test, the measured temperature profile at reflex provides a constraint on . Mathematically, we construct constraint function as follows. For any value of , by definition, the corresponding value of is the value of calculated based on and .
(3)
when is given, function is completely determined by the measured . Constraint functions based on 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 .
When the true relative density is unknown, we work with a trial relative density . We use to replace in (3) and construct trial constraint function from the test data. Note that the true values of satisfy only the true constraint function calculated using the true . When deviates from the true , the true values of are not on the trial constraint curve calculated using . Consequently, for a pair of trial constraint functions calculated using (based on measured of two distinct test conditions), their intersection is not at the true values of , 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 is incorrect. The specific behavior of test-condition-dependence of the intersection provides a venue for us to tune toward the true .
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 is governed by
where
· is the mass density of skin,
· is the specific heat capacity of the skin,
· K is the thermal conductivity of the skin,
· is the absorption coefficient of the skin,
· is the beam power density deposited on (absorbed into) the skin, and
· is the initial temperature of the skin.
In this idealized situation, the region of activated nociceptors ( ) is a cylinder with depth governed by . By definition, , and the reflex time are constrained by the temperature distribution, which we assume is measurable.
(4)
We consider a general 1-parameter form for the relative density of nociceptors
where can be any positive function. The activated depth and the dose both vary with the activation temperature , and are related by Equation (2) as
(5)
(6)
where A is the beam spot area, , and is the inverse function of . Since is positive, function is monotonically increasing and the inverse function is well-defined over the range of . By definition, functions and satisfy and . Recall that the dose threshold is defined as . Equation (6) gives
(7)
Equation (4) in combination with (7) provides a constraint on ,, and , which can be used for different purposes, depending on which parameters are known. When and are measured, (4) gives us a constraint on , and . On the other hand, when , and are given, (4) can be viewed as a governing equation for . This is useful, for example, for examining the behavior of vs. A.
In the analysis of subsequent sections, we need the expansions of and , and their derivatives. We now derive these expansions. We first write out the Taylor expansion of around .
The expansion of is derived from that of using an iterative method. It depends on which term in the expansion of is the leading non-zero term. Let . The inverse function maps u back to s. We discuss two cases.
· Case 1: . Function has the expansion
Based on that, we built an iteration formula:
The iteration gives us expansions of and
(8)
(9)
· Case 2: but . Function has the expansion
Based on that, we construct an iteration formula for :
which yields the expansion of in terms of u
Taking square roots of both sides, we obtain
(10)
(11)
with the mathematical results above, we study the behavior of vs A.
3. Analysis of Reflex Time vs. Beam Spot Area
When , and are given, (4) with (7) governs vs A. In our previous study [17], we scaled and shifted the physical temperature distribution to the normalized non-dimensional temperature .
(12)
where the non-dimensional depth and the non-dimensional time are defined as
The normalized temperature has the expression
(13)
It is important to notice that the normalized temperature is parameter-free. The non-dimensional version of (4) with (7) has the form
(14)
In the above, is the non-dimensional beam spot area. In the limit of , we have and the equation for becomes
where . Let . We have . We examine the asymptotic behavior of this convergence. We seek an expansion of the form
(15)
Expanding around and substituting (15), we get
(16)
Exponent and coefficient are determined using the leading term expansion of . The result depends on whether the nociceptor density at the skin surface is zero.
· Case 1: .
The expansion of given in (8). Substituting it into (16), we have
(17)
Here we have used and derived in our previous study.
· Case 2: but
The expansion of given in (10). Substituting it into (16), we arrive at
(18)
Returning to the physical quantities before the non-dimensionalization, the reflex time vs beam spot area is given by
(19)
4. Analysis of Constraint Functions on
When and are measured, (4) with (7) serves as a constraint on , and with beam spot area A as a parameter describing the test condition. We denote the constraint function as and use (4) to write it as
(20)
where is the temperature profile at reflex with beam spot area . We represent the effect of A via variable since is a smooth function of as . 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.
· : true value of coefficient .
· : a trial value of coefficient .
· : true values of model parameters and .
· : constraint function (20) calculated using trial value ; in , denotes the independent variable, and the dependent variable.
Note that the data is generated with the true value . In the calculation constraint function (20), we use the measured data and a trial value since is unknown. When , the constraint function shares one common intersection for all values of A:
(21)
when , in general is not on constraint curve , 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 around . For conciseness, we introduce and write (20) as:
We apply the chain rule to calculate partial derivatives.
(22)
(23)
Using these derivatives, we write out the expansion of .
(24)
The slope of vs. is , which is given in (22). In Appendix A, we show that the slope converges to zero as . It follows that
(25)
(25) shows that in the limit of , the constraint curve is a horizontal line at , independent of and . For finite A, and in (24) varies with and . We consider the intersection of and . The -coordinate of the intersection is . Let denote the -coordinate of the intersection. Solving for from (24) and (25) leads to
(26)
The dependence of on A is contained in given in (23). Using the expansion of derived in (73) in Appendix A, we investigate the behavior of at large A.
· Case 1: .
In case 1, as , the intersection given by (26) converges to the true value regardless of trial value .
The residual in convergence is proportional to and has the same sign as . More specifically, we have
(27)
· Case 2: but
(28)
In case 2, as , the intersection given by (26) converges to
which is proportional to trial value . The residual in convergence is proportional to and has the same sign as . Asymptotically, is
(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:
For type 1 nociceptor density, the relative density takes the parametric form
(30)
The graph of type 1 is plotted in Figure 17. It has the properties
5.1. Reflex Time
The reflex time vs beam spot area is given by (19) and (17), namely,
(31)
As A increases, converges to its limit with the residual proportional to 1/A2:
Figure 1 plots the relation between and A in two ways. Left panel: vs A. Right panel: vs. 1/A. In particular, the right panel confirms that the residual in the convergence of decays faster than 1/A for large A, as predicted in the analysis above.
Figure 1. The relation between reflex time ( ) and beam spot area (A) for type 1 nociceptor density (30) with . Left panel: vs A. Right panel: vs 1/A.
It is interesting to notice that expansion (31) is independent of . Consequently, the measurements of vs. A do not contain any information for estimating .
5.2. Constraint Function
We consider the intersection of the pair and . The -coordinate of the intersection, , is generally described by (26). For type 1 density, we have , and , and the specific expression of is given by (27).
(32)
At , the intersection is independent of A. When , the trend of vs. A tells us whether or .
· For , ascends toward from below as A increases.
· For , descends toward from above as A increases.
Figure 2 displays simulated for several values of A, respectively for and for . Here constraint function is based on test data (which is generated with true value ) and is calculated using formulation (20) with trial value . The trend of 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 down or up toward the true value depending on whether the calculated 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
Figure 2. Constraint curves for type 1 density (30). is based on test data generated with true value , and calculated using (20) with trial value . Left panel: . Right panel: .
Figure 3. vs. A, respectively for and for . Here is the intersection of and 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 denote the constraint function based on the test data (which is generated using type 1 density (30) with true value ) and calculated using framework (20) with the uniform trial density . Notice that the uniform density is a member of type 1 family (30) with . Thus, the two constraint functions and are related by . Let denote the -coordinate of the intersection of the pair and . Setting in (32), we obtain
(33)
Since the uniform density is parameter-free, the calculation of constraint function and intersection 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, and can be calculated. When the true underlying density is type 1 given in (30), result (33) predicts that converges to as with the difference proportional to 1/A. Figure 4 compares simulated and . The simulation results validate the theoretical prediction. In particular, varies linearly with 1/A. We fit function to data of vs 1/A. The fitting coefficients give us
Before we end this section, we clarify that result (33) predicts the behavior of constraint function 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 and is to identity the type of nociceptor density from the observed behavior of 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:
For type 2 nociceptor density, the relative density has the parametric form
(34)
Figure 4. Simulated results of for type 1 density given in (30). Here is based on the data from the true density (with ) but is calculated using the uniform trial density. Left panel: vs A. Right panel: vs. 1/A.
The graph of type 2 is shown in Figure 17. It is straightforward to see that
6.1. Reflex Time
The reflex time vs beam spot area is described by (19) and (18).
(35)
As A increases, converges to its limit with the difference proportional to 1/A.
Figure 5 plots the relation between and A in two ways. Left panel: vs A. Right panel: vs 1/A. In particular, the right panel confirms that 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/A2 for type 1 nociceptor density (30). To distinguish between these two density forms, we introduce an auxiliary quantity Q
(36)
The theoretical prediction above tells us
(37)
In (35), coefficient does contain . However, in (35) is tangled with
Figure 5. The relation between reflex time ( ) and beam spot area (A) for type 2 nociceptor density (34) with . Left panel: vs A. Right panel: vs 1/A.
other parameters. It is not possible to extract the value of unless all other parameters are known.
6.2. Constraint Function
We consider the intersection of the pair and . The -coordinate of the intersection, , is generally described by (26). For type 2 density, we have ,, and , and the specific expression of is given by (29).
At , the intersection is independent of beam spot area A. When , the trend of vs A tells us whether or .
· For , ascends toward as A increases.
· For , descends toward as A increases.
The increase/decrease trend of 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) , which varies with trial value , and 2) the residual in convergence is proportional to , instead of 1/A.
Figure 6 shows simulated for several values of A, respectively for and for . Here constraint function is based on test data (which is generated with true value ) and is calculated using formulation (20) with trial value . The trend of 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). is based on test data generated with true value , and calculated using (20) with trial value . Left panel: . Right panel: .
Figure 7. vs. A, respectively for and for . Here is the intersection of and 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 down or up toward the true value depending on whether the calculated increases or decreases with A.
6.3. Constraint Function Calculated Using the Uniform Density
When the type of parametric form of 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 denote the constraint function based on the test data (which is generated using type 2 density (34) with true value ) and calculated using framework (20) with the uniform trial density . 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 is not a special case of . To study the behavior of , we try to connect it to , the true constraint function. Both and 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:
· is the constraint function calculated using the correct density type (34) with the true parameter value . Operationally, the calculation of is not realistic unless we know . The theoretical advantage of is that it shares a common intersection for all values of A.
· is the constraint function calculated using the uniform trial density. Operationally, we can always calculate from test data. However, the point in general is not on the curve .
The two constraint functions above are unified in via , as described in (4).
(38)
The mapping between and , however, depends on the trial density. As a result, variable is related to differently in the two constraint functions. For clarity of the discussion, we use to denote the variable in , and use for the variable in . The difference lies in the trial density used in the formulation. For the true density, is related to in (5). For the uniform density, . Expressing in terms of or , we obtain
(39)
Combining (39) and (38), we write out and
(40)
Since is calculated using the true density, it satisfies
Taking the limit as yields
In (40), letting and using , we obtain
Let denote the intersection of the pair and . Both and are calculated by first mapping to and then mapping to . Given the measured temperature profile at reflex, Equation (38) maps to a unique , independent of the trial density. It follows that in (39), line 1 with and line 2 with are both equal to .
Therefore, and are related by
(41)
For type 2 density (34), the expansion of is given in (10) with and . Substituting the expansion of into (41) yields
(42)
Figure 8 compares simulated and . The simulation results confirm the theoretical prediction. In particular, varies linearly with . We fit function to the data of vs . The fitting coefficients give us
is the intersection of two constraint functions calculated from test data using the uniform trial density. Given the test data, the process of calculating is parameter-free. However, 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 increases linearly with , unbounded as . Result (42) for type 2 density is in sharp contrast
Figure 8. Simulated results of for type 2 density given in (34). Here is based on the data from the true density (with ) but is calculated using the uniform trial density. Left panel: vs A. Right panel: vs. .
with the situation for type 1 density, described in (33) of Section 5, where converges to the true value as . This difference in the behavior of 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:
We write type 3 density in the same parametric form as that of types 1 and 2.
(43)
where is the Heaviside step function.
(44)
The graph of type 3 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 at depth . Because of the discontinuity, the Taylor expansions of and derived in Section 3 are no longer valid. We write out and directly
The mapping between and is described by (7) and has the specific expression
(45)
7.1. Reflex Time
With the expression of in (45), Equation (14) for becomes
In the limit of , we have and the equation for converges to
Consider as a function of t with q as a parameter. Let be the inverse of and . As , we have
By definition, satisfies . We expand function around and write the equation for as
(46)
Substituting into (46) to calculate and , and then mapping back to the physical quantities, we obtain
As A increases, converges to its limit with the residual proportional to 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 and A in two ways. Left panel: vs A. Right panel: vs 1/A. In particular, the right panel confirms that is linear with respect to 1/A for large A, as predicted in the analysis above.
Figure 9. The reflex time ( ) as a function of beam spot area (A) for type 3 nociceptor density (43) with . Left panel: vs A. Right panel: vs. 1/A.
7.2. Constraint Function
For type 3 density, it is more sensible to use as the parameter since it has the clear physical meaning of the depth at which the nociceptor density jumps from 0 to . In the unified parametric form in (43), the generic parameter is . Mathematically, working with 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 and .
We adopt the general convention of using to denote the true value and to represent the variable. The general form of constraint function with trial value is given in (20). Using and , we write it as
(47)
In the limit of , converges to a horizontal line
At , the constraint function shares the common intersection for all A.
In particular, we have . We expand around and expand around to write them respectively as
(48)
(49)
Let denote the -coordinate of the intersection of and . is governed by equating the RHSs of (48) and (49).
Solving for yields
(50)
We expand with respect to into the form
(51)
Substituting the expansion into (50), we obtain
(52)
In Appendix B, we show that the coefficients satisfy
(53)
At , the intersection is independent of beam spot area A. When , the trend of vs A implies whether or .
· For , ascends toward as A increases.
· For , descends toward as A increases.
In terms of , the increase/decrease trend of 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: .
Figure 10 depicts simulated for several values of A, respectively for and for . Here constraint function is based on test data (which is generated with true value ) and is calculated using formulation (47) with trial value . The trend of 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 up or down toward the true value
Figure 10. Constraint curves for type 3 density (43). is based on test data generated with true value , and calculated using (47) with trial value . Left panel: . Right panel: . Notice that for type 3 density, is a horizontal line, whose height varies with .
Figure 11. vs A, respectively for and for . Here is the intersection of and from Figure 10.
depending on whether the calculated increases or decreases with A.
7.3. Constraint Function Calculated Using the Uniform Density
When the type of parametric form of 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 denote the constraint function based on the test data (which is generated using type 3 density (43) with true value ) and calculated using framework (47) with the uniform trial density . For type 3 parametric family (43), the uniform density is a special case at . However, expansions (48) and (49) are only for the case of near , away from the skin surface. At , the expansions will be different because of the insulated boundary condition . At , we have
We expand around . It follows that
Let denote the intersection of and . It satisfies
All derivatives are evaluated at . Solving for , we obtain
(54)
Result (54) predicts that when the true underlying density is type 3 given in (43), increases linearly with , unbounded as . This behavior is similar to that for type 2 density. Both type 3 and type 2 share the common feature of . Figure 12 compares simulated and . The simulation results confirm the theoretical prediction.
8. Type 4 Nociceptor Density:
We represent type 4 density in the same parametric form as that of types 1 - 3.
(55)
where is the Heaviside step function defined in (44). The graph of type 4 is shown in Figure 17. Type 4 relative density (55) has value 1 in , and drops down to value 0.5 for . 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 , 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 ( ).
Figure 12. Simulated results of for type 3 density given in (43). Here is based on the data from the true density (with ) but is calculated using the uniform trial density. Left panel: vs. A. Right panel: vs .
8.1. Reflex Time
Because of the discontinuity in density profile (8), the Taylor expansions of and derived in Section 3 are no longer valid. We write out and directly.
(56)
The activated depth at reflex given in (7) has the expression
(57)
The non-dimensional reflex time is related to A via in Equation (14).
(58)
Conventionally, we view 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 . We study function for type 4 density and connect it to its counterpart for the uniform density. Let
· : function A vs. for type 4 density
· : function A vs. for the uniform density
· : function vs. 1/A for type 4 density
· : function vs. 1/A for the uniform density.
For any given value of , the corresponding is completely determined by Equation (58), independent of the nociceptor density. We use to connect and . For type 4 density, is given in (57). Its inverse function is
(59)
For the uniform density, ,, and is given by
(60)
Combining (59) and (60), we can express in terms of .
(61)
(61) described the relation between and . Now we treat as the dependent variable and view it as a function of 1/A. (61) leads to
(62)
Equation (62) reveals that is obtained from by a piecewise linear scaling on the independent variable 1/A. Function is smooth. After the piecewise linear scaling, has a discontinuity in derivative. Figure 13 plots the relation between and A in two ways. Left panel: vs. 1/A. Right panel: derivative of vs. 1/A. In particular, the right panel verifies that vs. 1/A has a discontinuity in derivative, as predicted in the analysis above.
8.2. Constraint Function
For type 4 density, it is more sensible to choose as the parameter since it has the clear physical meaning of the depth at which the nociceptor density drops sharply from to . In the unified parametric form in (55), the generic parameter is . We adopt the general convention of using to denote the true value and to represent the variable. The general form of constraint function with trial value is given in (20). Using and given in (56), we write it as
Figure 13. The relation between reflex time ( ) and beam spot area (A) for type 4 density (55) with . Left panel: vs 1/A. Right panel: derivative of vs. 1/A.
(63)
when beam spot area A is sufficiently large to make , (57) gives 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 and is calculated with trial value . For any positive trial value and a finite interval of near in consideration, when A is sufficiently large, we have and
Here is the constraint function based on the data for the uniform density and is calculated using the uniform trial density. is independent of trial value , and passes through for all values of A.
In the limit of , we have . Let denote the -coordinate of the intersection of and . Our analysis above shows that
· For any , when A is sufficiently large.
· For any ,.
· For any , when A is sufficiently large.
More specifically, to have , we only need for near . The condition becomes . We conclude that
(64)
Thus, for large A, is independent of . To probe the position of trial value relative to true value , we need to examine the behavior of for . In this range of A, constraint function (63) takes the form
At true value , shares the common intersection for all A.
We expand around .
Setting the LHS to and solving for , we obtain
(65)
For and , the activated depth at reflex is . Substituting the expression of into (57) and solving for , we have
(66)
For and , the activated depth at reflex is . This may or may not be beyond depth . The case of corresponds to . Substituting the expression of into (57) and solving for yields
(67)
Similarly, the case of corresponds to . We get
(68)
Summarizing results for various cases, we write out a complete description of .
The case of ,
The case of ,
(69)
At , we have for all A. When , has non-trivial dependence on A. The trend of vs A tells us whether or .
· When , from small A to large A, starts below , decreases further below , and then reverts rapidly back to and stays there.
· When , from small A to large A, starts above , increases further above , and then reverts rapidly back to and remains there.
Figure 14 shows simulated for several values of A, respectively for and for . Here constraint function is based on test data (which is generated with true value ) and is calculated using formulation (63) with trial value . The trend of 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 up or down toward the true value depending on the behavior of the calculated vs. A.
8.3. Constraint Function Calculated Using the Uniform Density
When the type of parametric form of is unknown, we use the uniform
Figure 14. Constraint curves for type 4 density (55). is based on test data generated with true value , and calculated using (63) with trial value . Left panel: . Right panel: .
Figure 15. vs A, respectively for and for . Here is the intersection of and 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 denote the constraint function based on the test data (which is generated using type 4 density (55) with true value ) and calculated using framework (63) with the uniform trial density . For type 4 parametric family (55), the uniform density is a special case with large .
It follows that the behavior of is the same as that of for large , which we analyzed in the previous subsection.
Let denote the intersection of the pair and . Intersection has the same behavior as for large . Based on result (69) for in the previous subsection, we conclude for that
(70)
Result (70) predicts that when the true underlying density is type 4 given in (55), starts at slightly below for small A; decreases linearly toward as A increases; arrives at at and stays there for . The qualitative behavior of converging to for large A is similar for type 1 and type 4 densities. Figure 16 compares simulated and . The simulation results confirm the theoretical prediction.
Figure 16. Simulated results of for type 4 density given in (55). Here is based on the data from the true density (with ) but is calculated using the uniform trial density.
9. Summary
In this study, we investigate the nociceptor density of the form: , and its effect on heat-induced withdrawal reflex. We examine 4 types of illustrated in Figure 17. Each shown has a characteristic length 1. When scaled in the depth direction by parameter , each 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 for heat-sensitive nociceptors,
· the critical threshold on the equivalent activated volume, and
· the parameter 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 and the true scaling parameter are known for the nociceptor density, the test data allow us to construct true constraint functions on . Multiple constraint functions, obtained from tests at different values of beam spot area A, share a common point at the true value of . Parameters and are determined by finding the intersection of these distinct constraint functions.
Problem 2:
When the density type is given but the true scaling parameter is unknown, we construct trial constraint functions using trial values of . Let denote the trial constraint function that is based on the measured data (which is generated with true value ) and that is calculated using the given parametric form with trial value . When the trial value is different from the true value , in general, the true value 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 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: and . Let denote the intersection of the pair, which is affected by the trial value used in calculating the trial constraint functions and by the beam spot area A used in tests. Given a trial value , we examine the trend of vs A. For each of the 4 density types, as the beam spot area A increases, vs A demonstrates distinct trend of increasing or decreasing, respectively when and when . The trend of 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 vs. A, we tune the trial value accordingly to approach the true value . At the true value , is independent of A. Once we arrive at the true value , the subsequent procedure of finding parameters and is the same as described in Problem 1 above: the intersection of and gives us the true value of .
Problem 3:
When both the density type and the true scaling parameter are unknown, we construct trial constraint functions using the uniform trial density , which is parameter-free. The trial constraint function using the uniform density is denoted by . When the true density is not the uniform density, the true value is not on trial constraint function , and for different values of A trial constraint functions do not share a common point. Similar to what we did on trial constraint functions in Problem 2 (where the density type is given), here we look at the intersection of a pair of trial constraint functions: and . Let 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 vs. A. The trend behavior of vs. A is displayed for the 4 density types, respectively in Figure 4, Figure 8, Figure 12, and Figure 16. The behaviors of 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 vs. A, we also check the convergence pattern of reflex time 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 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 vs. A and 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 vs. A or the trend of 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 and 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.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
Cite this paper
Wang, H.Y., Burgei, W.A. and Zhou, H. (2020) Effect of Depth-Dependent Nociceptor Density on the Heat-Induced Withdrawal Reflex. Applied Mathematics, 11, 788-824. https://doi.org/10.4236/am.2020.118053
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
Appendix A
Derivation of in (25)
We show the convergence in three steps.
Step 1: We first look at the factor in the expression of in (22). As we have . We expand as . In both cases 1 and 2, we have . It follows that
(71)
Step 2: Next we look at the factor in the expression of in (22). With and the insulated boundary condition , we have
(72)
Taking the limit as in (22), and using results (71) and (72), we conclude
Step 3: We show that the term in (24) is bounded as . In the expression of in (23), we look at , the component that varies with . Using expansions (8) and (9) in case 1 or expansions (10) and (11) in case 2, we get
(73)
In both cases, we have .
Combining the results of three steps above, we conclude that for any trial value , as , the constraint function converges to , as described in (25) in the main text.
Appendix B
Derivation of property (53)
Property (53) describes coefficients in the expansion of with respect to . Recall that where . Function as given in (12) and (13), has the properties below
1) is always negative.
2) At any given depth y, the absolute value of increases with time t.
3) At any given time t, the absolute value of start with 0 value at , increases with y until the inflection point at , and then decreases to zeros beyond the inflection point. The depth of inflection point increases with t.
We consider the situation of . That is, the nociceptor density jump occurs before the inflection point of temperature profile at reflex. Combining properties 1 and 3 of with , we have
(74)
As beam spot area A increases decreases and is bounded from below by . In turn, is bounded by , reflecting that the presence of a skin layer with no nociceptor delays the occurrence of withdrawal reflex. Thus, we have
Applying property 2 of to (74) with , we obtain
(75)
The first part of (53), , follows from inequality (75) and definition of in (51). The second part of (53), , is verified in simulations.