Risk of Hearing Loss Caused by Multiple Acoustic Impulses in the Framework of Biovariability

We consider the hearing loss injury among subjects in a crowd with a wide spectrum of individual intrinsic injury probabilities due to biovariability. For multiple acoustic impulses, the observed injury risk of a crowd vs the effective combined dose follows the logistic dose-response relation. The injury risk of a crowd is the average fraction of injured. The injury risk was measured in experiments as follows: each subject is individually exposed to a sequence of acoustic impulses of a given intensity and the injury is recorded; results of multiple individual subjects were assembled into data sets to mimic the response of a crowd. The effective combined dose was adjusted by varying the number of impulses in the sequence. The most prominent feature observed in experiments is that the injury risk of the crowd caused by multiple impulses is significantly less than the value predicted based on assumption that all impulses act independently in causing injury and all subjects in the crowd are statistically identical. Previously, in the case where all subjects are statistically identical (i.e., no biovariability), we interpreted the observed injury risk caused by multiple impulses in terms of the immunity effects of preceding impulses on subsequent impulses. In this study, we focus on the case where all sound exposure events act independently in causing injury regardless of whether one is preceded by another (i.e., no immunity effect). Instead, we explore the possibility of interpreting the observed logistic dose-response relation in the framework of biovariability of the crowd. Here biovariability means that subjects in the crowd have their own individual injury probabilities. That is, some subjects are biologically less or more susceptible to hearing loss injury than others. We derive analytically the distribution of individual injury probability that produces the observed logistic dose-response relation. For several parameter values, we prove that the derived distribution is mathematically a proper density function. We further study the asymptotic approximations for the density function and discuss their significance in practical numerical computation with finite precision arithmetic. Our mathematical analysis implies that the observed logistic dose-response relation can be theoretically explained in the framework of biovariability in the absence of immunity effect.

crowd is the average fraction of injured. The injury risk was measured in experiments as follows: each subject is individually exposed to a sequence of acoustic impulses of a given intensity and the injury is recorded; results of multiple individual subjects were assembled into data sets to mimic the response of a crowd. The effective combined dose was adjusted by varying the number of impulses in the sequence. The most prominent feature observed in experiments is that the injury risk of the crowd caused by multiple impulses is significantly less than the value predicted based on assumption that all impulses act independently in causing injury and all subjects in the crowd are statistically identical. Previously, in the case where all subjects are statistically identical (i.e., no biovariability), we interpreted the observed injury risk caused by multiple impulses in terms of the immunity effects of preceding impulses on subsequent impulses. In this study, we focus on the case where all sound exposure events act independently in causing injury regardless of whether one is preceded by another (i.e., no immunity effect). Instead, we explore the possibility of interpreting the observed logistic dose-response relation in the framework of biovariability of the crowd. Here biovariability means that subjects in the crowd have their own individual injury probabilities. That is, some subjects are biologically less or more susceptible to hearing loss injury than others. We derive analytically the distribution of individual injury probability that produces the observed logistic dose-response relation.

Introduction
Hearing loss impacts about 360 million people worldwide (over 5 percent of the global population) [1]. One cause of hearing loss is exposure to loud sounds [2].
For example, impulse noise from explosive blasts and weapon firings can cause traumatic injuries to the auditory system [3] [4] [5].
In [4] 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 [6] we use the framework of immunity to interpret the empirical dose-response relation from [3] 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 [6] 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 dif-H. Wang  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" [7]. 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

Mathematical Formulation
In experiments [3], 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 [8], 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 ( Aeq8hr L ) is the best predictor of injury [8] [9]. In the injury model developed in [3], SELA 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 50 D 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 50 D are respectively, 0.1 α = and 50 161 D = dB [3]. It was observed that parameter 0.1 α = remains the same for PTS injuries of all cut-off levels while the median injury dose 50 D increases with the PTS cut-off level [3].
For a sequence of N identical impulses with SELA value = S, the effective combined dose is given by the dose combination rule [3]: The injury risk of the crowd caused by the N impulses is We express parameter a in terms of probability 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 probabili- 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 ( ) q ω denote the intrinsic non-injury probability of a random subject in the crowd, in responding to one impulse. Here ( ) q ω is a random variable and we include symbol ω in notation explicitly to indicate the presence of biovariability. Let On the other hand, the injury risk follows the observed logistic dose-response relation (3), which relates to the average non-injury fraction as ( ) In the next section, we derive an analytical expression for the distribution

Analytical Results
We solve for density ( )

Power Series Solution
We expand the right-hand side (RHS) of Equation (9) Recall the mathematical identity ( ) where the gamma function is defined as We use identity (11) to map each power of 1 N on the RHS of (10) back to a power of s.
Using this mapping, we express ( ) ( ) ( ) With the help of the Stirling formula, we calculate the radius of convergence using the ratio test: It follows that function ( ) g s is well defined by power series (13) for all values of s. To be a proper density function, however, ( ) g s must satisfy two more properties: We will rigorously demonstrate properties (16) and (17) for the cases of 1 2 η = and 1 3 η = . In addition, we will derive asymptotic approximations of ( ) g s for large s, which are essential for practical computation of ( ) g s 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 Integrating both sides of (13), we express ( ) We use scaling

The
which leads to a differential equation for ( ) Using the integrating factor method and applying the condition These two results lead directly to 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) To derive the asymptotic approximations of ( ) G s for large s, we change the integration variable to z defined by ( ) Differentiating with respect to s yields the asymptotics of ( ) g s for large s: for large s can be derived in a simple way, directly from differential Equation 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 ( ) G s at large s, based upon the power series around 0 s = .
Next, we study the case of 1 3 η = . We will see that the differential equation which gives us a differential equation for ( ) Using the integrating factor method and applying the condition we write the solution ( ) In Appendix A, we prove that ( ) These two results imply that ( ) ( ) is a mathematically proper density function.
To derive the asymptotics of ( ) G s for large s, we use a change of variable For large s, each term in the power series of ( ) Differentiating with respect to s yields the asymptotics of ( ) g s for large s: ( ) 4  5  7  8  10  11  3  3  3  3  3  3   2  4  10  28  80  2  1  2  1  2  1  3  3  9  9  27  27  3  3  3  3  3  3   s  s  s  s  s  s Again, the asymptotics of ( )

Asymptotics at Large s in the
which gives us a differential equation for ( ) We derive the asymptotics of ( ) Differentiating with respect to s yields the asymptotic of ( )

Significance of Asymptotics at Large s for Computing ( ) g s
In this subsection, we discuss the necessity of using the asymptotic approximations in calculating density ( ) g s for large s. The power series of ( ) g s around 0 s = , 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 10 16 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 This gives us the general magnitude and trend of ( ) g s . In particular, ( ) g s 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 We maximize ( ) f z with respect to z as a continuous variable. The deriva- The derivative indicates that for large s, the maximum of ( ) In summary, for large s, we estimated the net sum and the largest term in power series (39):  Net sum: g s 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 76 s = , the ratio of largest to net sum is 34 5.072 10 × 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 compu-  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 40 s = , the round-off errors completely wipe out the numerical accuracy.

Density of Individual Injury Probability in the Crowd
In this subsection, we write out ( ) The density of individual injury probability Here f k is the number of terms to include in summation to make the truncation error smaller than 10 −16 . f k is estimated using (51)  ( ) ( ) The threshold ( 16 s = ) for switching from power series around 0 s = 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 se-

Numerical Results and Discussion
In this section, we study the numerical behaviors of the power series around  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 ( ) g s . For large s, we turn our attention to asymptotics (57) of ( ) g s . The numerical errors in using asymptotics (57) to calculate ( ) g s 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 ( ) g s 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 10 s > , the 12-group approximation is more accurate while for 10 s < , the 6-group approximation is more accurate. With a fixed number of groups of terms, the approximation is more accurate as s is increased.  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 ( ) g s is to use the power series around 0 s = 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 16 s = 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 ( ) g s that is accurate to 10 −8 for all values of s in With the switching mechanism developed above for calculating ( ) g s , we study the distribution density and cumulative distribution function of the individual injury probability. In Figure 5    by the precision of the number representation system. In IEEE double precision, we adopt a switching numerical procedure for calculating ( )   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.