Hearing loss impacts about 360 million people worldwide (over 5 percent of the global population) . One cause of hearing loss is exposure to loud sounds . For example, impulse noise from explosive blasts and weapon firings can cause traumatic injuries to the auditory system   .
In  hearing loss injury experiments were conducted at the University of California, San Diego, using chinchilla as the animal subjects. Each chinchilla was individually exposed to impulse noises generated by a small shock tube placed at various distances to the animal to cover a range of intensity levels; after the exposure, the injury status of the animal was measured and recorded; results from multiple individual animals were assembled into data sets to mimic the collective response of a crowd in responding to impulse noises. The main reason to choose chinchilla as a surrogate animal model for human exposures is that hearing capabilities are very similar between chinchilla and humans. However, the threshold for susceptibility of injury for chinchilla is found to be lower than humans. Thus, the chinchilla data were then shifted by 28 dBA, representing scaling from chinchilla to humans. Based on the scaled data, an empirical injury model was developed for human exposed to multiple sound impulses of equal intensity.
In our recent work  we use the framework of immunity to interpret the empirical dose-response relation from  for exposure to multiple sound impulses. More specifically, from the observed dose-response relation, we derived the (negative) synergistic effect of a sequence of impulses in causing injury. We revealed that the phenomenological effect of a preceding sound exposure event on the subsequent sound exposure event is always immunity. The study of  was based on the assumption that all subjects in the crowd are statistically indistinguishable and thereby it excluded the effect of biological variability. In this paper we focus on the effect of biological variability instead. This approach will provide an alternative view of the empirical dose-response relation from a different theoretical angle.
Biological variability (or “biovariability” in short) is defined as “the natural variability in a lab parameter due to physiologic differences among subjects and within the same subject over time” . There are basically two types of biological variability, inter-individual and intra-individual. Inter-individual biovariability refers to “differences between subjects due to differences in diet, genetics or immune status” whereas intra-individual biovariability applies to “differences in the same subject over time due to diurnal cycles and other rhythms, biological repair mechanisms, dietary variables, aging, etc.” . In this paper, we study the inter-individual biovariability. Prominent examples of inter-individual biological variability include diversities in biologically inherited traits such as height, weight and skin color. Due to biological variability there is always a range of responses displayed in measurements collected from a crowd of animal or human subjects. Therefore, it is reasonable to expect that some people are biologically less or more vulnerable to hearing loss injury than others. For example, one of the physiologic causes for variability is that body temperature is known to affect sensitivity to noise .
The presence of biovariability prompts us to wonder: what kind of role does the biological variability play when a crowd of subjects is theoretically exposed to multiple acoustic impulses? To address this question, we investigate the feasibility of explaining the observed logistic dose-response relation in the framework of biovariability. We mathematically derive the distribution density of individual injury probability in a crowd that will reproduce the observed dose-response relation. For several parameter values relevant for applications, we show that the derived distribution is a proper density function. Further, we study the asymptotic approximations of the density function and their usage in numerical computation of density function in finite precision arithmetic. Our mathematical results indicate that it is theoretically possible to interpret the observed logistic dose-response relation in terms of biovariability. This study on biovariability, together with the previous study on immunity effects  , offers two complementary views of the observed dose-response relation. These two very different views represent the two extreme theoretical boundaries for understanding the observed dose-response relation.
2. Mathematical Formulation
In experiments  , the injury risk of a crowd caused by a sound exposure was measured by conducting separate tests on animal subjects individually and assembling individual results to mimic the collective response of a crowd. A sound exposure event is a complicated process, described by the acoustic waveform, intensity, and duration. In  , several candidates were compared as an effective single metric for characterizing the overall injury causing effect of a sound exposure event. It was found that 8-hour Equivalent A-weighted Level ( ) is the best predictor of injury  . In the injury model developed in  , SELA (A-weighted Sound Exposure Level) was selected as the dose, which is mathematically equivalent to , up to an additive constant. The observed injury risk follows the logistic dose-response relation:
For a crowd, the injury risk p is the average fraction of injured. In the logistic dose response relation (1), α is the shape parameter controlling the steepness of function, and is the median injury dose. For a crowd, the median injury dose is the dose level at which 50% of the crowd is expected to be injured. We consider the hearing loss injury risk for a crowd of subjects with a wide spectrum of individual injury probabilities due to biovariability. That is, some subjects are biologically less or more susceptible to hearing loss injury than others. At the apparent median injury dose for such a crowd, a particular subject's individual injury probability may be below or above 50% due to biovariability. For injury of permanent threshold shift (PTS) > 25 dB, the measured values of parameters α and are respectively, and dB . It was observed that parameter remains the same for PTS injuries of all cut-off levels while the median injury dose increases with the PTS cut-off level .
For a sequence of N identical impulses with SELA value = S, the effective combined dose is given by the dose combination rule  :
The injury risk of the crowd caused by the N impulses is
where parameters a and η are related to other parameters as
We express parameter a in terms of probability by solving ,
This expression states that parameter a is the odds of injury of a hypothetical subject in the crowd with the average injury probability, which is the ratio of average injury fraction to average non-injury fraction of the crowd. In the subsequent discussion, we shall refer to parameter a simply as the odds of injury of an average subject. We shall avoid the term “average odds” since it is not the average of odds. Rather it is the odds for a subject with the average injury probability.
In this study, we assume that N impulses act independently from each other in causing injury, regardless of whether one is preceded by another (i.e., no immunity effect). Under the assumption of independence, we explore the possibility of interpreting the observed logistic dose-response relation for a crowd in the framework of biovariability.
For mathematical convenience, we consider non-injury probability. When all events are independent of each other, the overall non-injury probability is simply the product of non-injury probabilities of individual events. Let denote the intrinsic non-injury probability of a random subject in the crowd, in responding to one impulse. Here is a random variable and we include symbol ω in notation explicitly to indicate the presence of biovariability. Let be the probability density of random variable . For one impulse, the average non-injury fraction of the crowd has the expression
For N impulses, the non-injury probability of an individual subject is and the average non-injury fraction of the crowd is the average of :
On the other hand, the injury risk follows the observed logistic dose-response relation (3), which relates to the average non-injury fraction as
Combining (6) and (7), we arrive at an equation for , the distribution density of non-injury probability of a random subject in the crowd.
In the next section, we derive an analytical expression for the distribution density based on Equation (8).
3. Analytical Results
We solve for density in Equation (8). To proceed, we make a change of variables . Since the non-injury probability is constrained to interval , variable has the domain . In terms of the new variable s, Equation (8) becomes
3.1. Power Series Solution
We expand the right-hand side (RHS) of Equation (9) as a power series of .
Recall the mathematical identity
where the gamma function is defined as
We use identity (11) to map each power of on the RHS of (10) back to a power of s.
Using this mapping, we express on the left-hand side (LHS) of (10) as a power series.
Function is the probability density of random variable . The derivation above tells us that the theoretical prediction matches the observed logistic dose-response relation if and only if the density of in the theoretical model has the power series in (13). It is not obvious at all whether or not function defined by power series (13) is a proper density function. Below we first show that function is well defined in its domain ( ). That is, power series (13) converges for all values of .
For large values of index k, the gamma function is well approximated by the Stirling formula   :
With the help of the Stirling formula, we calculate the radius of convergence using the ratio test:
It follows that function is well defined by power series (13) for all values of s. To be a proper density function, however, must satisfy two more properties:
We will rigorously demonstrate properties (16) and (17) for the cases of
and . In addition, we will derive asymptotic approximations of
for large s, which are essential for practical computation of in finite precision arithmetic. As we will see in the analysis, the proof approach differs quite significantly between these two cases. Thus, it is unlikely that the proof approach of either case can be directly extended to other values of η. Nevertheless, we conjecture that properties (16) and (17) are valid for all values of η.
To facilitate the analysis, we consider the cumulative distribution function (CDF) of random variable : . Integrating both sides of (13), we express as
We use scaling to get rid of parameter a. For simplicity and conciseness of presentation, we still denote the scaled variable by s and set . After scaling, the cumulative distribution function and the density function have the expressions
Our general strategy is to derive a differential equation for based on the power series around . From the differential equation, we solve for and , and then derive their asymptotic approximations at
large s. Building on the theoretical insight gained in the cases of and , we then extend the asymptotic analysis to the case of to derive an
essential numerical formula for practical computation of in finite precision arithmetic.
3.2. The Case of
Setting and differentiating (19) with respect to s, we have
which leads to a differential equation for :
Using the integrating factor method and applying the condition , we write the solution as
It is straightforward to verify that
• increases monotonically with s because
These two results lead directly to and .
Remark: It is fairly unusual for a solution of differential Equation (22) to converge to a constant level as s increases to . In the homogeneous version of differential Equation (22), solution exponentially diverges away from as s increases. It is the unique combination of the condition at and the RHS term in differential Equation (22) that makes . If we change either or both of these two, will diverge to infinity.
To derive the asymptotic approximations of for large s, we change the integration variable to z defined by . We write as
For the integral in (25), at large s, the domain of dominant contribution is a small neighborhood near . We expand around :
Each term in the power series of gives us
Putting these results together, we write out the asymptotics of for large s:
Differentiating with respect to s yields the asymptotics of for large s:
Interestingly, after we concluded , the asymptotics of for large s can be derived in a simple way, directly from differential Equation (22), by assuming that has an expansion of the form
Substituting expansion (28) iteratively into differential Equation (22) gives us
This set of equations yields the solution
The result of this simple method agrees with (26), which is rigorously derived without any additional assumption.
Remark: This simple method is quite powerful. By assuming the expansion form, we are able to derive the asymptotics of at large s, based upon the power series around .
Next, we study the case of . We will see that the differential equation for in the case of behaves very differently from that in the case of although the main conclusions regarding and remain the same.
3.3. The Case of
At , differentiating (19) with respect to s leads to
which gives us a differential equation for .
Using the integrating factor method and applying the condition , we write the solution as
In Appendix A, we prove that given in (31) satisfies
1) , and
2) increases monotonically.
These two results imply that is a mathematically proper density function.
To derive the asymptotics of for large s, we use a change of variable to rewrite solution (31) as
We expand and around :
For large s, each term in the power series of gives us
Substituting the power series into (32) and neglecting terms, we write out the asymptotics of for large s:
Differentiating with respect to s yields the asymptotics of for large s:
Again, the asymptotics of for large s can be derived directly from differential Equation (30) by assuming the expansion form (28). Below we use this approach to derive the asymptotics of for large s in the case of
. We select the parameter value to approximate the observed value
of in experiments.
3.4. Asymptotics at Large s in the Case of
At , differentiating (19) with respect to s leads to
which gives us a differential equation for :
We derive the asymptotics of for large s by assuming
• satisfies , and
• has the expansion form (28) for large s.
Substituting expansion (28) iteratively into differential Equation (36), we obtain
Differentiating with respect to s yields the asymptotic of for large s
3.5. Significance of Asymptotics at Large s for Computing
In this subsection, we discuss the necessity of using the asymptotic approximations in calculating density for large s. The power series of around , given in (20), has the nice property that it converges analytically for any s. In practical computation with finite precision arithmetic, however, the numerical convergence of (20) is limited by the finite precision and the limitation is fatal at large s. To see this limitation, we examine the net sum of the power series and compare the net sum with individual terms in the power series. When the largest term is 1016 times the net sum or more, implementing power series (20) with double precision arithmetic will lose all numerical accuracy in computing the net sum.
To simplify the discussion, we write power series (20) in a concise and abstract form
We study the case of as an example. The power series converges for
any s and the exact net sum is , which has the asymptotic approximation (38) at large s. We use the leading term in (38) to estimate the net sum for large s:
This gives us the general magnitude and trend of . In particular, decreases to zero as s increases. In power series (39), we write the absolute value of the k-th term as
The largest term in power series (39) is found by maximizing function with respect to index k. For mathematical convenience, we approximate this discrete maximization as a continuous maximization process. We maximize with respect to z as a continuous variable. Using the Stirling approximation (14), we write as
We maximize with respect to z as a continuous variable. The derivative of is
The derivative indicates that for large s, the maximum of is attained
approximately at . Consequently, is a fairly tight lower
bound on . That is, the true maximum may be a bit larger but not
too much larger than :
The largest term (in absolute value) in power series (39) is approximately
In summary, for large s, we estimated the net sum and the largest term in power series (39):
• Net sum: ,
• Largest term: .
The ratio of largest to net sum is
At , for example, the largest term in (39) is approximately ; the net sum of (39) is ; and the ratio of largest to net sum is . Recall that the machine precision of the IEEE double precision is about . Clearly, at , computing function by numerically summing terms in power series (39) will not yield any accuracy when implemented in double precision arithmetic. In addition, as s increases, the ratio of largest to net sum, given by Equation (46), grows exponentially with s. As a result, adopting a higher precision arithmetic will only enable us to accommodate a slightly larger range of s. For example, at , the ratio of largest to net sum is which will completely wipe out the numerical accuracy of quadruple precision arithmetic. In conclusion, although power series (39) converges theoretically for any s, it is not a practical numerical method for computing function at large s. For large values of s, the asymptotic approximation (38) is absolutely essential in providing a viable way for computing function .
Let us calculate , the index value at which the largest term occurs in power series (39). For large s, index satisfies approximately
Solving for from this equation, we obtain
For and , we have .
Figure 1 plots the k-th term of power series (39) vs index k for and . The largest term occurs at indices k = 282 and k = 283, which correspond very well to the predicted value of .
Figure 1. Plot of vs k for and . The inset shows details of the plot near the location of maximum.
It is worthwhile to find , the cut-off index beyond which all terms in the power series are smaller than a threshold value specified by a relative error tolerance (ε). Specifically, we find index such that
where ε is the relative error tolerance. With the expression of given in (41), we cast it approximately into a continuous problem
Applying the Stirling approximation (42) on the LHS and applying the leading term asymptotics of for large s (40) on the RHS, the constraint on becomes
For large s, the largest term of power series (39) occurs at . We expect the cut-off threshold to be a small multiple of , which suggests us to re-write the inequality for approximately in terms of :
Function is well approximated by its quadratic expansion.
Using this quadratic approximation, we solve for from inequality (50).
The corresponding index is
For , , and , we have .
Note that summing the first terms in the power series (39) is only a necessary condition for achieving a relative accuracy of in calculating . It is not a sufficient condition. In the numerical computation, the result is polluted by the accumulation of round-off errors caused by finite precision arithmetic. Thus, the numerical accuracy is limited by the machine precision of the number representation system multiplied by the largest term in the power series. As we showed above, with double precision arithmetic, even at a moderate , the round-off errors completely wipe out the numerical accuracy.
3.6. Density of Individual Injury Probability in the Crowd
In this subsection, we write out , the distribution density of individual injury probability, which shows the biological variability of the crowd in its susceptibility to hearing loss injury.
In all analysis above, function is the density of after scaling . Density of before the scaling is given by
Going back to non-injury probability , we write out the density of individual non-injury probability in terms of .
The density of individual injury probability is
where is the density of the post-scaling random variable . Function is calculated as follows.
• For , is numerically computed using the partial sum of first terms in the power series around given in (20).
Here is the number of terms to include in summation to make the truncation error smaller than 10−16. is estimated using (51) and (52). Since the estimate of is not very reliable for small s, we use at least 500 terms in summation.
• For , is evaluated using the first groups of terms in asymptotic approximation (38) for large s. In (38), only the first 3 groups are explicitly shown. Now let us write out a general formula for the partial sum of the first groups of terms in the asymptotic approximation. The general formula is in the form of double summation and is suitable for numerical implementation.
The threshold ( ) for switching from power series around to asymptotic approximation at large s is discussed in next section.
We need to point out that in asymptotic approximations (38) and (57), at a fixed value of s, if we include more and more groups of terms in the approximation, eventually the sum diverges. This is caused by the divergence of infinite series
In this infinite series, the gamma function coefficient increases much faster than exponential growth. At any fixed value of s, as index n increases, the n-th
term behaves like , which diverges to infinity. At a fixed and moderate
value of , however, asymptotic approximation (57) provides a vital numerical tool for computing function at large s when the power series around loses numerical accuracy due to accumulation of round-off errors.
4. Numerical Results and Discussion
In this section, we study the numerical behaviors of the power series around and the asymptotic approximations at large s. Based on the numerical results, we develop a numerical procedure for computing that is accurate for all values of s, small or large. The procedure is already mentioned in the previous section and is based on switching between power series (56) and asymptotics (57). Using this numerical procedure, we examine the distribution density of the individual injury probability in the crowd. We also plot the cumulative distribution function (CDF) of the injury probability because the CDF illustrates the effect of parameter a on the distribution of individual injury probability. Here parameter a is the observed odds of injury of an average subject in the crowd.
All simulations in this section are for the case of , which is meant to
approximate the observed value of in experiments. We first study the errors in numerical implementation of power series (56) in IEEE double precision. When sufficiently large number of terms are included in the power series summation to keep the theoretical truncation error well below the round-off errors, the numerical errors are mainly caused by accumulation of round-off errors. To calculate the numerical errors, we need a very accurate solution. For that purpose, we implemented all operations and functions (including , and ) in a very high numerical precision by representing each number using an array of 200 integers to achieve 1600 significant decimal digits. Numerical error is calculated as the difference between the numerical solution in double precision and the high precision solution. Figure 2 plots the absolute error and relative error vs s when the power series around is implemented in IEEE double precision. In the previous section, we derived that the largest term in power series (56) is
. Based on that, we expect the absolute numerical error to grow
approximately proportional to . In Figure 2, the gray dashed line is a fitting of the form to the numerical absolute errors, confirming the theoretical prediction.
The exponential growth of the accumulated round-off errors with s implies that, for large s, the numerical summation of power series (56) is not a good tool for computing . For large s, we turn our attention to asymptotics (57) of . The numerical errors in using asymptotics (57) to calculate are mainly due to the approximation errors at finite s. In the previous section, we showed that at any fixed value of s, if we include more and more terms in the asymptotics, it eventually diverges. At a finite s, the approximation error of asymptotics cannot be made arbitrarily small by including more terms even if we work in exact arithmetic (i.e., no round-off error). The approximation errors can only be reduced by increasing s. Thus, asymptotics (57) provides a good approximation of only when we use a fixed number of terms and when s is at least moderately large. Figure 3 displays the approximation errors in asymptotics (57) vs s when respectively 6 groups of terms and 12 groups of terms are used in approximation. For , the 12-group approximation is more accurate while for , the 6-group approximation is more accurate. With a fixed number of groups of terms, the approximation is more accurate as s is increased.
Figure 2. Numerical errors vs s when power series (56) of is implemented in IEEE double precision. Here the numerical errors are caused by the accumulation of round-off errors. The errors grow exponentially with s because the largest term in summation grows exponentially. The gray dashed line is a fitting of the form .
These numerical findings agree with the theoretical predictions in the previous section.
From the results of Figure 2 and Figure 3, we see that a good numerical strategy for calculating is to use the power series around when s is not too large and switch to the asymptotics when s is large enough. Figure 4 compares the errors of these two formulas: 1) the numerical error when implementing power series summation (56) in IEEE double precision, and 2) the approximation error in asymptotics (57). The accumulation of round-off errors in power series grows with s while the approximation error in asymptotics decays as s is increased. Figure 4 suggests as an optimal threshold for switching from power series (56) to asymptotics (57). When this switching strategy is implemented in IEEE double precision, it provides a numerical procedure for calculating that is accurate to 10−8 for all values of s in .
With the switching mechanism developed above for calculating , we study the distribution density and cumulative distribution function of the individual injury probability. In Figure 5, we plot density of the individual injury probability p in the crowd. Density is given in equation (55), expressed in terms of function . Figure 5 indicates that density diverges to infinity at and . Approximately, distribution can be viewed as the mix of 3 component distributions:
1) a δ-function distribution at ,
2) a nearly uniform distribution in , and
3) a δ-function distribution at .
Because of the divergence of density at and at , the weights of the two δ-function components cannot be readily read out from the plot of density
Figure 3. Approximation errors in asymptotics (57) vs s when respectively 6 groups of terms and 12 groups of terms are used in approximation. With each fixed number of groups of terms, the asymptotic approximation becomes more accurate at larger s.
Figure 4. Comparison of two errors: 1) the error in numerical summation of power series around , and 2) the approximation error in asymptotics for large s. The intersection of the two suggests a threshold for switching from one to the other.
. To see the weights, we examine , the cumulative distribution function (CDF) of p, which clearly illustrates the weight of a δ-function component as a jump in function value of the same magnitude.
The CDF, , has the expression
where is CDF of the post-scaling random variable , expressed as a power series around in (19). While theoretically, power series (19) converges for any s, the practical numerical convergence is limited
Figure 5. Density of the individual injury probability p for where parameter a is the observed odds of injury of an average subject in the crowd. The logarithmic scale is used in the vertical direction because diverges to infinity at and .
by the precision of the number representation system. In IEEE double precision, we adopt a switching numerical procedure for calculating , similar to that for calculating .
• For , is numerically computed using the partial sum of first terms in the power series around given in (19).
• For , is evaluated using the first groups of terms in the asymptotics for large s given in (37). We rewrite it as a general formula below.
Figure 6 compares plots of respectively for , , and . Parameter a is the observed odds of injury of an average subject in the crowd. As parameter a increases from small to one to large, the δ-function component at decreases and the δ-function component at increases by the same amount while the distribution in between is approximately unchanged.
In this study, we examined the observed logistic dose-response relation for hearing loss injury caused by exposure to multiple acoustic impulses, in the mathematical framework of biovariability. Previously, we interpreted the
Figure 6. Plots of , CDF of the individual injury probability, respectively for , , and . Parameter a is the observed odds of injury of an average subject in the crowd. As parameter a increases, its most prominent effect is to shift F(p) down nearly uniformly while keeping and .
observed dose-response relation based solely on the immunity effect under the assumption that all individual subjects are statistically identical (i.e., no biovariability). In the current study, we view the problem from a completely different angle; we explored the possibility of understanding the observed dose-response relation based solely on the biovariability of the crowd under the assumption that all sound exposure events act independently from each other in causing injury regardless of whether one event is preceded by another. We found that theoretically the biovariability alone is indeed capable of explaining the observed dose-response relation. We derived an analytical expression for the distribution density of individual injury probability in the crowd that produces the observed results. We also derived asymptotic approximations to the distribution density and the cumulative distribution function. The asymptotic approximations make it possible to compute the distribution density in finite precision arithmetic, in parameter regions where the analytical expression suffers catastrophically from loss of accuracy due to the accumulation of round-off errors. A robust numerical procedure for calculating the distribution density was developed based on switching between the analytical expression and the asymptotics. The resulting numerical procedure is accurate in all parameter regions, providing a practical tool for computing the distribution density. The mathematical framework constructed in the current study, along with the theoretical results and the numerical tools obtained, paves a pathway for further investigating the presence and effects of biovariability.
The authors would like to thank the Joint Non-Lethal Weapons Directorate of US Department of Defense 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 US Government.
We prove that given in (31) satisfies
1) , and
2) increases monotonically.
Item 1) is obtained directly by applying L’Hospital’s rule.
To prove 2), we introduce function
We write in terms of and define function :
We only need to show . The derivative of has the expression
From the definition, and satisfy
The sign of tells us that
In other words, attains its minimum at . Thus, to show , we only need to show . We use integration by parts to rewrite as
Using the sign of , we have
Combining these results, we conclude that and thus, function increases monotonically.