ABSTRACT

We study a general framework for assessing the injury probability corresponding to an input dose quantity. In many applications, the true value of input dose may not be directly measurable. Instead, the input dose is estimated from measurable/controllable quantities via numerical simulations using assumed representative parameter values. We aim at developing a simple modeling framework for accommodating all uncertainties, including the discrepancy between the estimated input dose and the true input dose. We first interpret the widely used logistic dose-injury model as the result of dose propagation uncertainty from input dose to target dose at the active site for injury where the binary outcome is completely determined by the target dose. We specify the symmetric logistic dose-injury function using two shape parameters: the median injury dose and the 10 - 90 percentile width. We relate the two shape parameters of injury function to the mean and standard deviation of the dose propagation uncertainty. We find 1) a larger total uncertainty will spread more the dose-response function, increasing the 10 - 90 percentile width and 2) a systematic over-estimate of the input dose will shift the injury probability toward the right along the estimated input dose. This framework provides a way of revising an established injury model for a particular test population to predict the injury model for a new population with different distributions of parameters that affect the dose propagation and dose estimation. In addition to modeling dose propagation uncertainty, we propose a new 3-parameter model to include the skewness of injury function. The proposed 3-parameter function form is based on shifted log-normal distribution of dose propagation uncertainty and is approximately invariant when other uncertainties are added. The proposed 3-parameter function form provides a framework for extending skewed injury model from a test population to a target population in application.

Keywords:

Dose Injury Relation, Dose Propagation Uncertainty, Median Injury Dose, 10 - 90 Percentile Width, Skewness, Mapping Injury Model from One Population to Another

1. Introduction

In many injury assessment situations, injury status of a subject is simply characterized in the form of binary outcome. For example, in a study of skull fracture injury related to highway traffic safety [1] and in a study of rib fracture injury caused by blunt-impact non-lethal weapons [2] , in each situation subjects tested are classified as either fractured or not fractured. Mathematically, occurrences of binary injury outcomes are statistically described by injury probability (also called injury risk). Let

• $\left(v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_k\right)$ be a list of input factors that affect the injury outcome,

• I be the binary injury outcome (random variable), and

• p be the corresponding injury probability: p = Pr (I = “injured”)

Here the binary injury outcome I is a random variable even when all input factors $\left(v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_k\right)$ are given and fixed. One approach of building a simple and practical model for assessing the injury risk is to use a single metric x to capture the overall effects of all input variables $\left(v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_k\right)$ [3] [4] . Quantity x is called the input dose, serving as the single metric best predictor of the injury probability. Input dose x may be one of the input variables $\left(v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_k\right)$ or a combination of these input variables. Depending on the application situations, input dose x is also called the determinant of injury, the risk factor, the exposure level, or the predictor variable [3] [4] .

When the input dose x is directly controllable and measurable, an experimental data set consists of m entries, each containing a measured value of input dose and the corresponding binary injury outcome in an independent trial:

$\text{Data}=\left\{\left(x_j\mathrm{,}{I}_j\right)\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}m\right\}\mathrm{.}$ (1)

Injury models are constructed in the general form of injury probability vs input dose.

$p\left(x\right)=\text{injury}\text{probability}\text{at}\text{input}\text{dose}x$

In many application situations, however, the input dose is not directly measurable. For example, for bone fracture injuries, we may use the stress at the impact site as the input dose. But it is difficult to measure directly the stress at impact site. In a study of behind-armor blunt trauma (BABT) [5] and a study of human body response to blunt impacts using advanced total body model (ATBM) [6] , an estimated value of stress caused by the impact is calculated via computer simulations. The estimation is based on the measured mass and velocity of projectile and using representative median material properties of the projectile, the subject body and the armor. When the true input dose is not directly measurable, an experimental data set contains pairs of estimated input dose and the corresponding binary injury outcome:

$\text{Data}=\left\{\left(x_j^{\left(est\right)}\mathrm{,}{I}_j\right)\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}m\right\}\mathrm{.}$ (2)

In these situations, practical injury models are constructed in the form of injury probability vs the estimated input dose

$p\left(x^{\left(est\right)}\right)=\text{injury}\text{probability}\text{at}\text{estimated}\text{input}\text{dose}{x}^{\left(est\right)}$

The estimated input dose, in general, is different from the true input dose, and the discrepancy between the two is population dependent since the actual material properties of individual subjects are different from the selected representative material properties and are population dependent. In addition, the relation of injury probability vs true input dose is also population dependent because the material properties of subjects significantly affect the injury outcome even when the true input dose is fixed. For example, at a fixed impact force, the injury probability varies considerably among groups of different ages, among groups of different body types, body sizes and body compositions. The experimentally established relation of injury probability vs estimated input dose is heavily influenced by the particular population tested. As a result, applying the injury model established for one population, straightforwardly without modification, to assess the injury risk of a different population will inevitably lead to large errors. In many applications, however, we face exactly this task: we are given an injury model established on a particular test population and we need to predict the injury risk of a different population. For example, a data set for human forearm fracture was assembled in [7] from drop test results conducted on PMHS forearms from cadaver donors of average age 55. The purpose of assembling the data set, however, is to build an injury model for assessing the risk of forearm fracture among a population of live human subjects with an age distribution significantly different from that of cadaver donors. In this study, we develop a simple mathematical framework for this task. The key idea is based on interpreting the probabilistic injury model as the consequence of dose propagation uncertainty from input dose to target dose at the active site for injury where the binary outcome is uniquely determined by the target dose. The framework of dose propagation uncertainty makes it mathematically convenient to accommodate different uncertainties associated with different populations. The formulation developed provides a mechanism of mapping injury function from one population to another by simply updating the model parameters.

2. Mathematical Formulation

We first review the logistic model for binary outcomes [8] . Note that the injury probability p is not directly observable in experiments unless we repeat the experiment a large number of times at each fixed value of input dose. The logistic regression model was designed by working with the hidden injury probability p and considering the logit function of p, which is defined as the

logarithm of the injury odds: $\text{logit}\left(p\right)\equiv log\left(\frac{p}{1-p}\right)$. In the logistic model, $\text{logit}\left(p\right)$ is postulated to be a linear function of input dose x,

$\text{logit}\left(p\right)=\alpha \left(x-D_{\text{50}}\right)$ (3)

Writing probability p as a function of x, we obtain the logistic dose response relation

$p={\text{logit}}^{-1}\left(\alpha \left(x-D_{\text{50}}\right)\right)=\frac{1}{1+exp\left(-\alpha \left(x-D_{\text{50}}\right)\right)}\equiv f_{\text{logistic}}\left(x\right)$ (4)

We write the linear function in (3) as $\alpha \left(x-D_{\text{50}}\right)$ so that constant $D_{\text{50}}$ has the meaning of the median injury dose, at which the injury probability is 50%: $p\left(D_{50}\right)=50\%$ [9] . In general, we introduce $D_{\eta }$ to denote the dose value at which $p\left(D_{\eta }\right)=\eta \%$. For example, $D_{\text{90}}$ is the dose with $p\left(D_{90}\right)=90\%$. Coefficient $\alpha$ controls the steepness of transition (i.e., the sensitivity of injury probability with respect to dose change). We define the width of injury function as

$W\equiv D_{\text{90}}-D_{\text{10}}$

Conceptually, the width W is not the 10 - 90 percentile range of x since dose x is not a random output of an experiment; it is the controlled input. However, if we view the injury function as the cumulative distribution function (CDF) for x and draw random samples of x based on the CDF, then the width W is indeed the 10 - 90 percentile range of random samples drawn. For simplicity, we shall call W the 10 - 90 percentile width even though x is not a random variable. In the logistic model, the width W is inversely proportional to coefficient $\alpha$.

$W\equiv D_{90}-D_{10}=\frac{2\ln \left(9\right)}{\alpha }$ (5)

We point out that the steepness coefficient $\alpha$ exists only in the logistic model. In contrast, the width of injury function is universally defined and meaningful for all injury models. To facilitate the comparison of various models, we shall use the width (W) instead of coefficient $\alpha$ whenever it is appropriate to do so. The logistic model in terms of shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ has the expression.

$p=\frac{1}{1+exp\left(-\frac{2ln\left(9\right)}{W}\left(x-D_{\text{50}}\right)\right)}\equiv f_{\text{logistic}}\left(x\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}\mathrm{,}W\right)$ (6)

Logistic model is widely used as a phenomenological model for binary outcomes [3] [4] [10] . In this study, we interpret it and approximate it in the framework of dose propagation uncertainty from input dose to target dose. The key assumption in our interpretation is that there is an active site at which target dose Z uniquely determines the binary outcome I. Mathematically, target dose Z at the active site has the features described below:

Ÿ Binary outcome I is the indicator function of $Z>z^{(c)}$

$I=\{\begin{array}{ll}1,\hfill & \text{if}Z>z^{\left(c\right)}\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$ (7)

where $z^{\left(c\right)}$ is the critical threshold for target dose in transition from non-injury to injury. The transition is a discontinuous jump with respect to target dose Z at the active site. However, with respect to the input dose x that is away from the active site, the injury probability vs x generally is a smooth and gradual transition.

• The target dose Z is caused by the input dose x. While in most experiments the input dose x can be controlled, at least to some extent, the target dose Z is neither directly observable nor directly controllable.

• For a given input dose x, the corresponding target dose Z is a random variable, reflecting the uncertainty in the propagation from input dose to target dose.

We use an example to illustrate the propagation from input dose to target dose.

Example: Passing exam vs amount of study time

In this example, the input dose x is the amount of study time. Note that although the target dose Z is caused by the input dose x, quantities Z and x may have different physical dimensions. For passing an exam, the target dose Z is the effective fraction of actual exam contents correctly completed in the exam by the student. We use a flow chart to show a possible propagation from input dose to target dose.

x = the nominal amount of study time invested

®Z₁ = effective amount of study time

affected by the student’s attentiveness, effciency, and overall load

®Z₂ = amount of course contents learned

affected by the student’s prior preparation and ability of memorizing key items

®Z₃ = fraction of actual exam contents learned

affected by the exam scope and weighting of components in exam

®Z = effective fraction of actual exam contents correctly completed

affected by the student’s general health condition on exam day, and ability of working under time pressure and in presence of noise/disturbance (8)

Mathematically, we write the target dose explicitly as $Z\left(x\mathrm{,}\omega \right)$, emphasizing that Z is a random variable depending on the input dose x and depending on the random factor $\omega$ in the dose propagation. The probability that a given input dose x leads to injury is

$\text{Pr}\left(I=“\text{injured}”\right)=\text{Pr}\left(Z\left(x\mathrm{,}\omega \right)>z^{\left(c\right)}\right)$ (9)

We consider two models for uncertainty in dose propagation: 1) target dose $Z\left(x\mathrm{,}\omega \right)$ has a normal distribution; and 2) target dose $Z\left(x\mathrm{,}\omega \right)$ is expressed in terms of a normally distributed intermediate variable. For example, intermediate variable $Y\left(x\mathrm{,}\omega \right)\equiv ln\left(Z\left(x\mathrm{,}\omega \right)-x_0\right)$ has a normal distribution, and target dose $Z\left(x\mathrm{,}\omega \right)$ is a shifted log normal distribution, expressed in terms of intermediate variable $Y\left(x\mathrm{,}\omega \right)$ as $Z\left(x\mathrm{,}\omega \right)=exp\left(Y\left(x\mathrm{,}\omega \right)\right)+x_0$.

3. Logistic Dose-Injury Relation Interpreted as Normally Distributed Target Dose

We model the target dose as proportional to the sum of the input dose and an additive Gaussian noise.

$Z\left(x\mathrm{,}\omega \right)=r\times \left(x-\mu +\sigma \epsilon \right)$

where $\epsilon ~N\left(\mathrm{0,1}\right)$, a standard normal random variable. We scale target dose Z and the associated critical threshold $z^{\left(c\right)}$ to make $r=1$ by changing the physical unit for measuring z-values, or equivalently by changing the physical unit for measuring x-values. Thus, we set $r=1$ and proceed with

$Z\left(x\mathrm{,}\omega \right)=x-\mu +\sigma \epsilon$ (10)

In this section, we first examine the dose-response relation for normally distributed dose uncertainty, which is the probit model [11] . Then we discuss how to accommodate different uncertainties corresponding to different populations, including how to incorporate additional uncertainties into the dose-response relation.

3.1. Dose-Response Relation

The binary injury outcome is governed by the sign of random variable

$Z\left(x,\omega \right)-z^{\left(c\right)}=x-z^{\left(c\right)}-\mu +\sigma \epsilon$ (11)

The injury probability (p) corresponding to input dose x is

$p=\text{Pr}\left(\left(x-z^{\left(c\right)}-\mu +\sigma \epsilon \right)>0\right)=\text{Pr}\left(-\epsilon <\frac{x-z^{\left(c\right)}-\mu }{\sigma }\right)$

Recall that the cumulative distribution function (CDF) of standard normal is given by the error function, $\text{erf}\left(u\right)$, which is defined as

$\text{erf}\left(u\right)\equiv \frac{2}{\sqrt{\text{π}}}{\displaystyle {\int }_0^u}exp\left(-s^2\right)\text{d}s$

The dose response relation for normally distributed target dose Z has the expression:

$p=\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{x-z^{\left(c\right)}-\mu }{\sqrt{2}\sigma }\right)\equiv f_{\text{normal}}\left(x\right)$ (12)

We approximate dose-response relation (12) using the logistic function form (4) with tunable parameters $D_{\text{50}}$ and $\alpha$. First, we match the two functions at $p=50\%$ to obtain $D_{\text{50}}=z^{\left(c\right)}+\mu$. To simplify the search for optimal $\alpha$, we apply the transformation

$x_{\text{new}}=\frac{x_{\text{old}}-D_{\text{50}}}{\sigma }$

After the transformation, (4) and (12) as functions of $x_{\text{new}}$ have standard forms:

$f_{\text{logistic}}\left(x_{\text{old}}\right)\mapsto f_L\left(x\right)\equiv \frac{1}{1+exp\left(-{\alpha }^{\prime }x\right)}$ (13)

$f_{\text{normal}}\left(x_{\text{old}}\right)\mapsto f_N\left(x\right)\equiv \frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{x}{\sqrt{2}}\right)$ (14)

where the scaled coefficient ${\alpha }^{\prime }$ is related to $\alpha$ by ${\alpha }^{\prime }=\sigma \alpha$. For conciseness, we denote $x_{\text{new}}$ simply as x. The task of approximating (12) with (4) is reduced to finding an optimal value of ${\alpha }^{\prime }$ such that the distance between $f_N\left(x\right)$ and $f_L\left(x\right)$ is minimized. Using numerical optimization, we find that the best approximation is achieved at ${{\alpha }^{\prime }}_{\text{opt}}=1.701$.

Figure 1 compares functions (14) and (13) at ${{\alpha }^{\prime }}_{\text{opt}}=1.701$. It is clear that the two functions are very good approximations of each other. The maximum difference is bounded by 0.01 (i.e., difference in predicted injury probability is less than 1%). With that error tolerance, the logistic model and the normal distribution model can practically substitute each other. In other words, the widely used logistic model can be viewed as a very good approximation of the normal distribution model, which was derived based on normally distributed dose propagation uncertainty from input dose to target dose.

Models (13) and (14) are nevertheless mathematically different. When the data set of binary injury outcomes (I) is sufficiently large, eventually, the two models will be distinguishable. Let m be the number of samples in the data set. We look into the question of how large m needs to be in order to statistically distinguish the two models. We consider a collection of independent data sets, each of the form

$D=\left\{\left(x_j\mathrm{,}{I}_j\right)\mathrm{,}j=\mathrm{1,2,}\cdots \mathrm{,}m\right\}$

where $x_j$ is the input dose of the j-th experiment and $I_j$ the corresponding binary injury outcome. To test if the two models are statistically distinguishable, we generate data sets according to the normal distribution model $f_N\left(x\right)$ in (14). In all data sets, values of input dose $\left\{x_j\right\}$ are uniformly distributed in $\left[-\mathrm{3,3}\right]$, and for each input dose $x_j$ the corresponding binary injury outcome $I_j$ is sampled using injury probability $f_N\left(x_j\right)$.

Given data set D, the log-likelihood for a general probability function $f\left(x\right)$ is

$\mathcal{L}\left(f(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)\equiv \frac{1}{m}{\displaystyle \sum _{j=1}^m}\left(I_{j}log\left(f\left(x_j\right)\right)+\left(1-I_j\right)log\left(1-f\left(x_j\right)\right)\right)$ (15)

We use log-likelihood (15) to compare models $f_N\left(x\right)$ and $f_L\left(x\right)$. Since $f_N\left(x\right)$ is the exact probability model for the data set while $f_L\left(x\right)$ is a slightly incorrect model, the difference in log-likelihood $\mathcal{L}\left(f_N(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)-\mathcal{L}\left(f_L(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)$ is expected to be positive. However, due to randomness of data sets, the difference in log-likelihood between two models fluctuates from one date set to another. We examine the sample distribution of differences in log-likelihood based on $N=100000$ independent data sets. Figure 2 plots the histograms of $\mathcal{L}\left(f_N(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)-\mathcal{L}\left(f_L(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)$ for various values of m, the size of each data set.

To clarify, here N is the number of data sets used in each histogram and m is the number of binary outcomes in each data set. In Figure 2, each sample of difference in log-likelihood requires one data set. That is why we use $N=100000$ independent data sets to plot each histogram.

Suppose we use the sign of $\mathcal{L}\left(f_N(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)-\mathcal{L}\left(f_L(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)$ to classify data sets as the normal distribution model (positive sign) or as the logistic model (negative sign). All data sets examined in Figure 2 are generated based on the normal distribution model. Thus, data sets with $\mathcal{L}\left(f_N(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)-\mathcal{L}\left(f_L(\cdot )\mathrm{<mo>\; |\; </mo>}D\right)<0$ will be falsely identified as the logistic model (false negative). In Figure 2, all counts to the left of the dashed black line in each histogram correspond to false negative identification. For data sets of $m=500$ samples each (top left panel), the false negative rate is 25.26%. For $m=1000$ (top right panel), the false negative rate decreases to 19.44%. When the sample size is increased to $m=2000$ (bottom left panel), the false negative rate falls to 12.49%. Finally, when the sample size is doubled again to $m=4000$ (bottom right panel), the false negative rate drops down to 5.63%. Based on the simulation results, we see that to reduce the false negative rate to less than 20%, for example, we need to work with data sets, each consisting of $m=1000$ samples. This is above the typical sample size of data sets for injury models. Thus, in real applications, the normal distribution model (14) and logistic model (13) are practically the same unless we work with injury data sets of very large sample size.

We go back to the pre-transformation logistic model, function (4) specified by steepness coefficient $\alpha$, and function (6) specified by width W. The corresponding optimal values for $\alpha$ and for W are respectively

$\begin{array}{l}{\alpha }_{\text{opt}}=\frac{{{\alpha }^{\prime }}_{\text{opt}}}{\sigma }=\frac{1.701}{\sigma }\\ W_{\text{opt}}=\frac{2ln\left(9\right)}{{\alpha }_{\text{opt}}}=\frac{2ln\left(9\right)}{1.701}\sigma =2.584\sigma \end{array}$ (16)

Since the 10 - 90 percentile width is well defined for all injury functions, we choose to specify the logistic model using width W instead of coefficient $\alpha$. We conclude that normal distribution model (12) based on dose propagation uncertainty is practically equivalent to logistic model (6) with shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ given by

$W=2.584\sigma ,D_{50}=z^{\left(c\right)}+\mu$ (17)

(17) describes the best approximation to the normal distribution model (12) from the logistic model family (6). The best approximation is obtained numerically by minimizing the distance between the two functions (Figure 1). Alternatively, a straightforward approximation can be written out by simply matching the widths of two injury functions. The width of normal distribution model is given by the inverse error function

$W_{\text{normal}}=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)\sigma =2.563\sigma$

Notice that the two widths, the width of normal distribution model $W_{\text{normal}}$ and the width of its best logistic model approximation $W_{\text{opt}}$, are indeed very close to each other. We will use these two interchangeably.

Similar to the situation of logistic model, the normal distribution model is also completely specified by the shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$. It has the form

$p=\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{2{\text{erf}}^{-1}\left(0.8\right)\cdot \left(x-D_{\text{50}}\right)}{W}\right)\equiv f_{\text{normal}}\left(x\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}\mathrm{,}W\right)$ (18)

where shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ are related to parameters of dose propagation uncertainty in (17). It should be pointed out that in general, the target dose Z is hidden, not observable or controllable; none of parameters $z^{\left(c\right)}$, $\mu$ or $\sigma$ is directly observable. These are internal quantities in the mathematical model, explaining why the injury probability follows the normal distribution model (12). In an idealized situation, the input dose x should be a controllable/measurable variable, and shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ may be determined from experimental measurements. In realistic applications, however, the true input dose x may not be directly measurable, which we will discuss in next subsection. At the end of this subsection, we summarize the normal distribution model for dose propagation uncertainty, and its connection to the widely used logistic model.

Summary of the injury model based on dose propagation uncertainty

• We select the physical unit for measuring the target dose Z such that in the absence of dose propagation uncertainty, target dose Z is the same as input dose x:

${Z|}_{\left(\text{zerouncertainty}\right)}=x\mathrm{.}$

• In the normal distribution model, the difference between target dose and input dose is an additive Gaussian noise:

$Z\left(x\mathrm{,}\omega \right)-x=-\mu +\sigma N\left(\mathrm{0,1}\right)\mathrm{.}$

• The binary injury outcome is completely determined by the condition $Z\left(x,\omega \right)>z^{\left(c\right)}$ where $z^{\left(c\right)}$ is the critical threshold for target dose Z.

• The probability of injury caused by the input dose x is described by the CDF of normal distribution. Practically the injury probability is very well approximated by the widely used logistic dose-response relation.

• As given in (17), the median injury dose of injury function is the critical threshold for the target dose, shifted by the bias in the dose propagation:

$D_{50}=z^{\left(c\right)}+\mu ,$

and the width of injury function is proportional to the uncertainty in dose propagation (standard deviation of the Gaussian noise):

$W=2.584\sigma$

The larger the uncertainty, the more spread out the injury function is.

• In terms of shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$, the logistic model is expressed in (6); the normal distribution model is given in (18).

Next, we study how to incorporate additional uncertainties in the framework of dose-response relation, and how to model a new population with different uncertainty.

3.2. Effects of Additional Uncertainties

In the previous subsection, we interpreted the dose-response relation as a consequence of dose propagation uncertainty. In this subsection we study how to incorporate additional uncertainties by changing the shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ in logistic model (6) or in normal distribution model (18).

We start by considering a homogeneous population consisting of statistically identical subjects, which means quantities $z^{\left(c\right)}$, $\mu$ and $\sigma$ are fixed and stay the same for all subjects in the population. In a homogeneous population, the dose propagation uncertainty is statistically the same for all subjects. Its effect is already reflected in the dose response relation specified by shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$, which are related to internal parameters $\left(z^{\left(c\right)}\mathrm{,}\mu \mathrm{,}\sigma \right)$ in (17). In particular, the width W is proportional to the standard deviation of uncertainty. If there is no uncertainty present in the dose propagation, the dose-response relation would be a sharp transition (a step function).

Now we consider a more realistic situation: a heterogeneous population consisting of subjects with variable critical threshold $z^{\left(c\right)}$, denoted here in the new setting as $Z^{\left(c\right)}$, following the convention of using uppercase letters for random variables. In addition to the uncertainty in $Z^{\left(c\right)}$, the input dose x may not be directly measurable. In some situations, the input dose x is not directly measured; instead, input dose x is derived from a controllable/measurable variable y. In these situations, the value of input dose x is calculated via computer simulations from measurable quantities using idealized representative properties of subjects, such as the 50-percentile properties of the general population [5] [6] . We use the example below to illustrate the situation of controllable variable y vs true input dose $X^{\left(true\right)}$ vs estimated input dose $X^{\left(est\right)}$. Consider the experiment in which we test the shatter resistance of a product by dropping it from a specified height. In this example, the various quantities in the model are described as follows:

• The height y is the controllable/measurable variable.

• The estimated input dose $X^{\left(est\right)}$ is the impact force calculated in a computer simulation from height y using the representative median properties, such as the weight of the product, the aerodynamic properties, the mechanical properties of the product and the ground surface, and the orientation angle of the product at impact.

• The true input dose $X^{\left(true\right)}$ is the actual impact force, which in general is different from the estimated input dose $X^{\left(est\right)}$. The difference $\left(X^{\left(true\right)}-X^{\left(est\right)}\right)$ depends on how much the true properties deviate from the selected representative properties. The distribution of difference varies from one population to another.

• The target dose $Z\left(X^{\left(true\right)}\mathrm{,}\omega \right)$ is the maximum stress at the most vulnerable part of the product.

The bottom line is that the true input dose $X^{\left(true\right)}$ is a random variable when the controllable variable y is specified. We model the difference $\left(X^{\left(true\right)}-X^{\left(est\right)}\right)$, the dose propagation uncertainty $\left(Z\left(X^{\left(true\right)}\mathrm{,}\omega \right)-X^{\left(true\right)}\right)$, and the critical threshold $Z^{\left(c\right)}$ as additive Gaussian noises. Mathematically, we formulate the problem as

$\left(Z\left(X^{\left(true\right)}\mathrm{,}\omega \right)-X^{\left(true\right)}\right)=-{\mu }_1+{\sigma }_1{\epsilon }_1$ (19)

$\left(X^{\left(true\right)}-X^{\left(est\right)}\right)=-{\mu }_2+{\sigma }_2{\epsilon }_2$ (20)

$Z^{\left(c\right)}=z^{\left(c\right)}+{\mu }_3+{\sigma }_3{\epsilon }_3$ (21)

where $\left\{{\epsilon }_1\mathrm{,}{\epsilon }_2\mathrm{,}{\epsilon }_3\right\}$ are i.i.d. samples of $N\left(\mathrm{0,1}\right)$. The binary injury outcome is governed by the sign of random variable

$Z\left(X^{\left(true\right)},\omega \right)-Z^{\left(c\right)}=X^{\left(est\right)}-z^{\left(c\right)}-\left({\mu }_1+{\mu }_2+{\mu }_3\right)+\sqrt{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}\epsilon$ (22)

At a given value of $X^{\left(est\right)}$, random variable $\left(Z\left(X^{\left(true\right)}\mathrm{,}\omega \right)-Z^{\left(c\right)}\right)$ has the same mathematical form as random variable $\left(Z\left(x\mathrm{,}\omega \right)-z^{\left(c\right)}\right)$ in (11). As a result, the injury probability vs the estimated input dose has the expression

$\begin{array}{c}{p|}_{X^{\left(est\right)}=x}=\text{Pr}\left(\left[x-z^{\left(c\right)}-\left({\mu }_1+{\mu }_2+{\mu }_3\right)+\sqrt{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}\epsilon \right]>0\right)\\ =\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{x-z^{\left(c\right)}-\left({\mu }_1+{\mu }_2+{\mu }_3\right)}{\sqrt{2}\sqrt{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}}\right)\end{array}$ (23)

Injury function (23) has the same form as (12). Thus, ${p|}_{X^{\left(est\right)}=x}$ is described by the normal distribution model with shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ given as follows.

${p|}_{X^{\left(est\right)}=x}=f_{\text{normal}}\left(x\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}\mathrm{,}W\right)$ (24)

$W=2.584\sqrt{{\sigma }_1^2+{\sigma }_2^2+{\sigma }_3^2}$

$D_{50}=z^{\left(c\right)}+\left({\mu }_1+{\mu }_2+{\mu }_3\right)$

In a well controlled lab setting, the true input dose $X^{\left(true\right)}$ is measurable. For example, in experiments of male forearm fracture [12] , a cylinder of specified mass is dropped from a specified height along a vertical track onto the PMHS forearm sample. Both the forearm sample and the cylinder are connected to accelerometers, allowing accurate measurements of the dynamic impactor load and the support loads. In addition, in situ strain gauges are used to record time series of strains at various locations during the loading. In this idealized setting, there is no measurement error in $X^{\left(true\right)}$. The injury probability as a function of the true input dose, ${p|}_{X^{\left(true\right)}=x}$, can be determined from the observed binary injury outcomes vs measured values of $X^{\left(true\right)}$. Injury function ${p|}_{X^{\left(true\right)}=x}$ follows the normal distribution model with shape parameters $\left(D_{\text{50}}^{\left(0\right)}\mathrm{,}{W}^{\left(0\right)}\right)$ given below.

${p|}_{X^{\left(true\right)}=x}=f_{\text{normal}}\left(x\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}^{\left(0\right)}\mathrm{,}{W}^{\left(0\right)}\right)$ (25)

$W^{\left(0\right)}=2.584\sqrt{{\sigma }_1^2+{\sigma }_3^2}$

$D_{\text{50}}^{\left(0\right)}=z^{\left(c\right)}+\left({\mu }_1+{\mu }_3\right)$

With this formulation, we can map back and forth between injury functions ${p|}_{X^{\left(true\right)}=x}$ and ${p|}_{X^{\left(ext\right)}=x}$. We can also revise the injury function ${p|}_{X^{\left(ext\right)}=x}$ measured on one population to construct the injury function for a different population. We now discuss these two problems.

Problem 1:

Suppose we are given an injury model ${p|}_{X^{\left(true\right)}=x}$, specified by shape parameters $\left(D_{\text{50}}^{\left(0\right)}\mathrm{,}{W}^{\left(0\right)}\right)$. The given injury function is for an idealized setting where the true input dose is directly measured. Our goal is to extend the given injury function ${p|}_{X^{\left(true\right)}=x}$ to predict the injury probability, ${p|}_{X^{\left(ext\right)}=x}$, as a function of estimated input dose for the same population when the true input dose is not measurable.

Solution:

Injury function ${p|}_{X^{\left(true\right)}=x}$ is specified by shape parameters $\left(D_{\text{50}}^{\left(0\right)}\mathrm{,}{W}^{\left(0\right)}\right)$ given in (25) while injury function ${p|}_{X^{\left(ext\right)}=x}$ is specified by shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ given in (24). Combining (25) with (24), we write $\left(D_{\text{50}}\mathrm{,}W\right)$ as an update on $\left(D_{\text{50}}^{\left(0\right)}\mathrm{,}{W}^{\left(0\right)}\right)$.

$\begin{array}{l}{W}^2=W_0^2+{\left(2.584\right)}^2{\sigma }_2^2\\ D_{50}=D_{50}^{\left(0\right)}+{\mu }_2\end{array}$ (26)

Problem 2:

Suppose we are given an injury model ${p|}_{X^{\left(ext\right)}=x}$, specified by shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$. The given injury function is established based on measurements of a heterogeneous population, labeled population 1. Population 1 is characterized by uncertainties in the input dose estimation and in the critical threshold, as described in (20) and (21)

$\left(X^{\left(true\right)}-X^{\left(est\right)}\right)=N\left(-{\mu }_2\mathrm{,}{\sigma }_2^2\right)$

$Z^{\left(c\right)}=z^{\left(c\right)}+N\left({\mu }_3,{\sigma }_3^2\right)$

Now consider a different heterogeneous population, labeled population 2, with uncertainties described by

$\left(X^{\left(true\right)}-X^{\left(est\right)}\right)=N\left(-{\tilde{\mu }}_2,{\tilde{\sigma }}_2^2\right)$

$Z^{\left(c\right)}=z^{\left(c\right)}+N\left({\tilde{\mu }}_3,{\tilde{\sigma }}_3^2\right)$

Here we assume that the propagation uncertainty from true input dose to target dose $\left(Z\left(X^{\left(true\right)}\mathrm{,}\omega \right)-X^{\left(true\right)}\right)=N\left(-{\mu }_1\mathrm{,}{\sigma }_1^2\right)$ is statistically the same for the two populations. Our goal is to predict the injury function ${\left({p|}_{X^{\left(ext\right)}=x}\right)}_2$ for population 2 based on the given injury function ${p|}_{X^{\left(ext\right)}=x}$ for population 1.

Solution:

Injury function ${p|}_{X^{\left(ext\right)}=x}$ for population 1 is specified by shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$ while injury function ${\mathrm{<mo\; stretchy="false">\; (\; </mo>}{p|}_{X^{\mathrm{<mo\; stretchy="false">\; (\; </mo>}est\mathrm{<mo\; stretchy="false">\; )\; </mo>}}\mathrm{<mo>\; =\; </mo>}x}\mathrm{<mo\; stretchy="false">\; )\; </mo>}}_2$ for population 2 is specified by shape parameters $\left({\tilde{D}}_{\text{50}}\mathrm{,}\tilde{W}\right)$. We write $\left({\tilde{D}}_{\text{50}}\mathrm{,}\tilde{W}\right)$ as an update on $\left(D_{\text{50}}\mathrm{,}W\right)$ to take into account the differences in uncertainties between the two populations.

$\begin{array}{l}{\tilde{W}}^2=W^2+{\left(2.584\right)}^2\left({\tilde{\sigma }}_2^2-{\sigma }_2^2+{\tilde{\sigma }}_3^2-{\sigma }_3^2\right)\\ {\tilde{D}}_{50}=D_{50}+\left({\tilde{\mu }}_2-{\mu }_2\right)+\left({\tilde{\mu }}_3-{\mu }_3\right)\end{array}$ (27)

4. Dose-Injury Function for Target Doze of Log-Normal Distribution

For the discussion below, we adopt the normal-distribution model as the base formulation, switching away from the logistic model. There are several reasons behind the switching.

• The normal-distribution model is based on 1) viewing the binary injury outcome as completely determined by the target dose at the active site, 2) explaining the randomness in injury outcome as the consequence of uncertainty in dose propagation from input dose to target dose, and 3) modeling the dose propagation uncertainty as an additive Gaussian noise. This interpretation is both theoretically and operationally appealing.

• Mathematically, the injury function form of normal-distribution model is exactly invariant when additional normally distributed noise/uncertainty is incorporated into the model.

• We will study dose-injury models based on normally distributed intermediate variable. Mathematically, such an injury model is conveniently treated as a transformation of the normal-distribution model since the target doze is expressed as a function of the normally distributed intermediate variable.

• As we demonstrated in the previous section, the logistic model is practically equivalent to the normal-distribution model with the same shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$.

We first recall the function form of the normal-distribution model. In terms of internal variables $\left(\sigma \mathrm{,}\mu \mathrm{,}{z}^{\left(c\right)}\right)$, it is given by (12). In terms of shape parameters $\left(D_{\text{50}}\mathrm{,}W\right)$, it is expressed in (18). Geometric quantities $D_{\text{50}}$, $D_{\text{10}}$, $D_{\text{90}}$ and W of the injury function are related to internal variables $\left(\sigma \mathrm{,}\mu \mathrm{,}{z}^{\left(c\right)}\right)$ as

$\begin{array}{l}{D}_{\text{50}}=z^{\left(c\right)}+\mu \\ D_{\text{10}}=D_{\text{50}}-{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \\ D_{\text{90}}=D_{\text{50}}+{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \\ W\equiv D_{\text{90}}-D_{\text{10}}=2\sqrt{2}{\text{erf}}^{-1}\left(0.8\right)\sigma \end{array}$ (28)

Because of the symmetry of error function $\text{erf}\left(z\right)$, the normal distribution model (18) is symmetric around the median injury dose $D_{\text{50}}$ :

$\text{Symmetry}\mathrm{<mo>\; :\; </mo>}{f}_{\text{normal}}\left(D_{\text{50}}+x\right)-\frac{1}{2}=\frac{1}{2}-f_{\text{normal}}\left(D_{\text{50}}-x\right)$ (29)

We now study a skewed injury function that breaks this symmetry. Consider the situation where the target dose $Z\left(x\mathrm{,}\omega \right)$ has a log-normal distribution

$Z\left(x\mathrm{,}\omega \right)=x\cdot exp\left(-\mu +\sigma \epsilon \right)$

Again $\epsilon ~N\left(\mathrm{0,1}\right)$ is a standard normal random variable. In this case, $ln\left(Z\left(x\mathrm{,}\omega \right)\right)$ and $ln\left(x\right)$ are simply related by an additive Gaussian noise.

$ln\left(Z\left(x\mathrm{,}\omega \right)\right)=ln\left(x\right)-\mu +\sigma \epsilon$ (30)

If we use $ln\left(x\right)$ and $ln\left(Z\left(x\mathrm{,}\omega \right)\right)$ to measure, respectively, the input dose and the target dose, then the injury probability vs $ln\left(x\right)$ follows the same function form as (12) with $\left(x\mathrm{<mo>\; ;\; </mo>}{z}^{\left(c\right)}\right)$ replaced by $\left(ln\left(x\right)\mathrm{<mo>\; ;\; </mo>}ln\left(z^{\left(c\right)}\right)\right)$ :

$p\left(x\right)=f_{\text{normal}}\left(ln\left(x\right)\mathrm{<mo>;</mo>}ln\left(z^{\left(c\right)}\right)\right)=\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{ln\left(x\right)-ln\left(z^{\left(c\right)}\right)-\mu }{\sqrt{2}\sigma }\right)$ (31)

We examine the injury probability as a function of the original input dose x. The purpose is to investigate 1) under what condition the injury probability vs x can be approximated by the symmetric normal-distribution model, and 2) when the normal distribution approximation is invalid, what additional parameter we need to introduce to describe the injury function for the original input dose x.

Since the injury probability vs $ln\left(x\right)$ follows the normal distribution model (12), we use results (28) for (12) to write out $D_{\text{50}}^{\left(ln\right)}$, $D_{\text{10}}^{\left(ln\right)}$ and $D_{90}^{\left(ln\right)}$ for quantity $ln\left(x\right)$.

$\begin{array}{l}{D}_{\text{50}}^{\left(\ln \right)}=ln\left(z^{\left(c\right)}\right)+\mu \\ D_{\text{10}}^{\left(\ln \right)}=D_{\text{50}}^{\left(\ln \right)}-{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \\ D_{\text{90}}^{\left(\ln \right)}=D_{\text{50}}^{\left(\ln \right)}+{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \end{array}$ (32)

The corresponding $D_{\text{50}}$, $D_{\text{10}}$ and $D_{\text{90}}$ for quantity x are

$\begin{array}{l}{D}_{\text{50}}=exp\left(D_{\text{50}}^{\left(\ln \right)}\right)=z^{\left(c\right)}\cdot {\text{e}}^{\mu }\\ D_{\text{10}}=exp\left(D_{\text{10}}^{\left(\ln \right)}\right)=D_{\text{50}}\cdot exp\left(-{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\\ D_{\text{90}}=exp\left(D_{\text{90}}^{\left(\ln \right)}\right)=D_{\text{50}}\cdot exp\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\end{array}$ (33)

In this case, it is clear that $\left(D_{90}-D_{50}\right)>\left(D_{50}-D_{10}\right)$. The injury probability vs quantity x is not exactly symmetric around $D_{\text{50}}$. We introduce a measure of skewness to represent the asymmetry of injury probability vs quantity x.

$\gamma \equiv \ln \left(\frac{D_{90}-D_{50}}{D_{50}-D_{10}}\right)$ (34)

Specifically, $\gamma$ defined above measures the skewness of interval $\left[D_{\text{10}}\mathrm{,}{D}_{\text{90}}\right]$ around $D_{\text{50}}$.

• When $\gamma =0$, interval $\left[D_{\text{10}}\mathrm{,}{D}_{\text{90}}\right]$ is symmetric around $D_{\text{50}}$.

• When $\gamma >0$, we have $\left(D_{90}-D_{50}\right)>\left(D_{50}-D_{10}\right)$, which implies that the upper half (above $D_{\text{50}}$ ) of injury function is flatter than the lower half (below $D_{\text{50}}$ ).

• When $\gamma <0$, we have $\left(D_{90}-D_{50}\right)<\left(D_{50}-D_{10}\right)$, and that the upper half of injury function is steeper than the lower half.

Skewness $\gamma$ is an indicator of how well the injury function for x can be approximated by the symmetric normal distribution model. For a target dose of log-normal distribution, the skewness is $\gamma ={\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma$. When $\sigma$ is small, the skewness $\gamma \approx 0$, and the injury function is nearly symmetric around $D_{\text{50}}$. When ${\sigma }^2>0$, the skewness $\gamma$ is positive, and in (31) the injury probability as a function of x is not symmetric. In this case, the injury function is characterized by three shape parameters: $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$.

$\begin{array}{l}{D}_{\text{50}}=z^{\left(c\right)}\cdot {\text{e}}^{\mu }\\ W\equiv D_{\text{90}}-D_{\text{10}}=D_{\text{50}}\cdot 2sinh\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\\ \gamma \equiv ln\left(\frac{D_{\text{90}}-D_{\text{50}}}{D_{\text{50}}-D_{\text{10}}}\right)\mathrm{<mo>\; =\; </mo>}{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \end{array}$ (35)

Notice that even though expressions of $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ in (35) contain three variables $\left(\mu \mathrm{,}\sigma \mathrm{,}{z}^{\left(c\right)}\right)$, two variables $\mu$ and $z^{\left(c\right)}$ appear only as a combination $\left(z^{\left(c\right)}\cdot {\text{e}}^{\mu }\right)$ in $D_{\text{50}}$. Mathematically, the three shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ are completely specified by $\left(D_{\text{50}}\mathrm{,}\sigma \right)$, and thus, have only two degrees of freedom. As a result, the three shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ cannot be set independently of each other. For example, in (35) when $\gamma$ is small, the width W will be small unless the median dose $D_{\text{50}}$ is large. Formulation (35), based on target dose of log-normal distribution (30), cannot accommodate any negative skewness ( $\gamma <0$ ). It cannot even accommodate the simple symmetric case of $\gamma =0$ with finite $W>0$ and $D_{50}<+\infty$. We like to revise the formulation and construct an injury model in which the three shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ can be set independently of each other.

5. A Dose-Injury Model with Skewness Based on a Normally Distributed Intermediate Variable

We construct a model that accommodates the median injury dose ( $D_{\text{50}}$ ), the width (W) and the skewness ( $\gamma$ ) as 3 independent parameters. In previous section, we studied the formulation based on target dose of log-normal distribution, in which the skewness is always positive and the 3 shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ are not independent of each other. A log-normal random variable can be viewed as the exponential of normal random variable. To accommodate negative skewness and to make $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ independent of each other, we extend the formulation to the case of target dose being a more general function of normal random variable.

We consider the situation where the dose propagation uncertainty is an additive Gaussian noise in quantity $ln\left|x-x_0\right|$ with $x_0$ as a new tunable parameter. The target dose $Z\left(x\mathrm{,}\omega \right)$ and the input dose x are related by

$Z\left(x,\omega \right)-x_0=\left(x-x_0\right)\exp \left(-\mu +\sigma \epsilon \right)$

In this setting, $\left(Z\left(x\mathrm{,}\omega \right)-x_0\right)$ has the same sign as $\left(x-x_0\right)$. The domain of x is divided by $x_0$ into two regions: $x>x_0$ and $x<x_0$. Only the region containing the critical threshold $z^{\left(c\right)}$ will be relevant for the injury model. The other region of x produces target dose $Z\left(x\mathrm{,}\omega \right)$ always above or always below $z^{\left(c\right)}$. For example, when $x_0>z^{\left(c\right)}$, only the region $x<x_0$ is relevant for the injury model; the region $x>x_0$ leads to target doze $Z\left(x,\omega \right)>x_0>z^{\left(c\right)}$ and thus, leads to an injury probability of 100%. We discuss separately the case of $x_0>z^{\left(c\right)}$ and the case of $x_0<z^{\left(c\right)}$.

5.1. Case 1: $x_0<z^{(c)}$

In this case, the region $x<x_0$ yields target dose $Z\left(x,\omega \right)<x_0<z^{\left(c\right)}$ and an injury probability of 0%. We focus on the region $x>x_0$, the relevant region for the injury model. The logarithm of shifted target dose $ln\left(Z\left(x\mathrm{,}\omega \right)-x_0\right)$ and logarithm of shifted input dose $ln\left(x-x_0\right)$ are related by an additive Gaussian noise.

$\ln \left(Z\left(x,\omega \right)-x_0\right)=ln\left(x-x_0\right)-\mu +\sigma \epsilon$ (36)

where $\epsilon ~N\left(\mathrm{0,1}\right)$. We apply the shift $x_{new}=x_{old}-x_0$ on all dose quantities (including $D_{\text{50}}$ and $z^{\left(c\right)}$ ). After the shift, problem (36) above is exactly the same as problem (30) in the previous section. It follows that the injury probability has the same function form as (12) with $\left(x\mathrm{<mo>\; ;\; </mo>}{z}^{\left(c\right)}\right)$ replaced by $\left(ln\left(x-x_0\right)\mathrm{<mo>\; ;\; </mo>}ln\left(z^{\left(c\right)}-x_0\right)\right)$

$\begin{array}{c}p\left(x\right)=f_{\text{normal}}\left(ln\left(x-x_0\right)\mathrm{<mo>\; ;\; </mo>}ln\left(z^{\left(c\right)}-x_0\right)\right)\text{for}{x}_0<z^{\left(c\right)}\\ =\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{ln\left(x-x_0\right)-ln\left(z^{\left(c\right)}-x_0\right)-\mu }{\sqrt{2}\sigma }\right)\end{array}$ (37)

Based on results (33) and (35), we write out $\left(D_{\text{50}}\mathrm{,}\gamma \mathrm{,}W\right)$ for injury function (37).

$\begin{array}{l}{D}_{\text{50}}-x_0=\left(z^{\left(c\right)}-x_0\right)\cdot {\text{e}}^{\mu }\\ D_{\text{10}}-x_0=\left(D_{\text{50}}-x_0\right)\cdot exp\left(-{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\\ D_{\text{90}}-x_0=\left(D_{\text{50}}-x_0\right)\cdot exp\left(+{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\\ W\equiv D_{\text{90}}-D_{\text{10}}=\left(D_{\text{50}}-x_0\right)\cdot 2sinh\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\\ \gamma \equiv ln\left(\frac{D_{\text{90}}-D_{\text{50}}}{D_{\text{50}}-D_{\text{10}}}\right)={\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \end{array}$ (38)

Note that both $z^{\left(c\right)}$ and $D_{\text{50}}$ are on the right side of $x_0$ in the case of $x_0<z^{\left(c\right)}$. As we will see, $D_{\text{50}}$ and $z^{\left(c\right)}$ are always on the same side of $x_0$. With Formulas (38) for the case of $x_0<z^{\left(c\right)}$, we can accommodate shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ with positive skewness $\gamma >0$. Specifically, at any fixed $\mu$, for each given set of $\left(D_{50},W,\gamma >0\right)$ there is a unique corresponding set of $\left(x_0\mathrm{,}\sigma \mathrm{,}{z}^{\left(c\right)}\right)$.

$\begin{array}{l}\sigma =\frac{\gamma }{{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}}\\ \left(D_{\text{50}}-x_0\right)=\frac{W}{2sinh\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)}\mathrm{,}{x}_0=D_{\text{50}}-\left(D_{\text{50}}-x_0\right)\mathrm{.}\\ z^{\left(c\right)}=x_0+\left(D_{\text{50}}-x_0\right){\text{e}}^{-\mu }\end{array}$ (39)

This works for any positive skewness $\gamma >0$, corresponding to the situation where the injury probability has a flatter rise above the median injury dose $D_{\text{50}}$ than below it.

To accommodate negative skewness $\gamma <0$, however, we need $x_0>z_{\left(c\right)}$.

5.2. Case 2: $x_0>z^{(c)}$

In this case, we focus on the region $x<x_0$ since the region $x>x_0$ yields target dose $Z\left(x,\omega \right)>x_0>z^{\left(c\right)}$ and an injury probability of 100%. The target dose and input dose are related by

$-ln\left(x_0-Z\left(x\mathrm{,}\omega \right)\right)=-ln\left(x_0-x\right)+\mu +\sigma \epsilon$ (40)

where $\epsilon ~N\left(\mathrm{0,1}\right)$. Here we consider quantity $-ln\left(x_0-Z\left(x\mathrm{,}\omega \right)\right)$ with the negative sign because it is an increasing function of $Z\left(x\mathrm{,}\omega \right)$. Injury occurs when the target dose is above the critical threshold: $Z\left(x,\omega \right)>z^{\left(c\right)}$, which translates to

$\text{Injury}\text{probability}=\text{Pr}\left(-ln\left(x_0-Z\left(x\mathrm{,}\omega \right)\right)>-ln\left(x_0-z^{\left(c\right)}\right)\right)$

The injury probability has the expression

$\begin{array}{c}p\left(x\right)=f_{\text{normal}}\left(-ln\left(x_0-x\right)\mathrm{<mo>\; ;\; </mo>}-ln\left(x_0-z^{\left(c\right)}\right)\right)\text{for}{x}_0>z^{\left(c\right)}\\ =\frac{1}{2}+\frac{1}{2}\text{erf}\left(\frac{-ln\left(x_0-x\right)+ln\left(x_0-z^{\left(c\right)}\right)+\mu }{\sqrt{2}\sigma }\right)\end{array}$ (41)

Notice that (31) with quantities denoted by (' ) and (41) are connected by transformation

$\left(x-x_0\right)=\frac{-1}{x^{\prime }},\left(z^{\left(c\right)}-x_0\right)=\frac{-1}{z^{\left(c\right)}{}^{\prime }},\mu =-{\mu }^{\prime }$

We use results (33) and (35) for injury function (31) to write out $\left(D_{\text{50}}\mathrm{,}\gamma \mathrm{,}W\right)$ for (41).

$\begin{array}{l}{D}_{\text{50}}-x_0=-\left(x_0-z^{\left(c\right)}\right)\cdot {\text{e}}^{\mu }<0\\ D_{\text{10}}-x_0=\left(D_{\text{50}}-x_0\right)\cdot exp\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)<0\\ D_{\text{90}}-x_0=\left(D_{\text{50}}-x_0\right)\cdot exp\left(-{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)<0\\ W\equiv D_{\text{90}}-D_{\text{10}}=-\left(D_{\text{50}}-x_0\right)\cdot 2sinh\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)>0\\ \gamma \equiv ln\left(\frac{D_{\text{90}}-D_{\text{50}}}{D_{\text{50}}-D_{\text{10}}}\right)=-{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma <0\end{array}$ (42)

In the case of $x_0>z^{\left(c\right)}$, both $z^{\left(c\right)}$ and $D_{\text{50}}$ are on the left side of $x_0$. With Formulas (42) for the case of $x_0>z^{\left(c\right)}$, we can accommodate shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ with negative skewness $\gamma <0$. Specifically, at any fixed $\mu$, for each given set of $\left(D_{50},W,\gamma <0\right)$ there is a unique corresponding set of $\left(x_0\mathrm{,}\sigma \mathrm{,}{z}^{\left(c\right)}\right)$.

$\begin{array}{l}\sigma =\frac{\left(-\gamma \right)}{{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}}\\ \left(D_{\text{50}}-x_0\right)=\frac{-W}{2sinh\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)}\mathrm{,}{x}_0=D_{\text{50}}-\left(D_{\text{50}}-x_0\right)\mathrm{.}\\ z^{\left(c\right)}=x_0+\left(D_{\text{50}}-x_0\right){\text{e}}^{-\mu }\end{array}$ (43)

This works for $\gamma <0$, which indicates that the injury probability has a steeper rise above the median injury dose $D_{\text{50}}$ than below it.

Next we combine the results of $x_0<z^{\left(c\right)}$ and $x_0>z^{\left(c\right)}$ to derive a unified formulation for accommodating shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ regardless of the sign of $\gamma$.

5.3. A Unified Formulation for All Values of Skewness

In the previous sub-section, we studied models based on target dose of shifted log normal distribution with shift as a parameter. We now synthesize the results obtained to develop a unified formulation of injury function in which the 3 shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ can be specified independently.

First, we show that at any fixed value of $\mu$, there is one-to-one correspondence between $\left(x_0\mathrm{,}\sigma \mathrm{,}{z}^{\left(c\right)}\right)$ and $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$. For any given set of shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ regardless of the sign of $\gamma$, we combine results (39) and (43) to write out the corresponding $\left(x_0\mathrm{,}\sigma \mathrm{,}{z}^{\left(c\right)}\right)$.

$\begin{array}{l}\sigma =\frac{\left|\gamma \right|}{{\text{erf}}^{-1}\left(0.8\right)\sqrt{2}}\\ x_0=D_{\text{50}}-\frac{W}{2sinh\left(\gamma \right)}\\ z^{\left(c\right)}=x_0+\left(D_{\text{50}}-x_0\right){\text{e}}^{-\mu }\end{array}$ (44)

Conversely, for any given set of $\left(x_0\mathrm{,}\sigma \mathrm{,}{z}^{\left(c\right)}\right)$, we combine results (38) and (42) to write out the corresponding shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$.

$\begin{array}{l}{D}_{\text{50}}=x_0+\left(z^{\left(c\right)}-x_0\right)\cdot {\text{e}}^{\mu }\\ W=\left|D_{\text{50}}-x_0\right|\cdot 2sinh\left({\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \right)\\ \gamma =\text{sign}\left(z^{\left(c\right)}-x_0\right)\cdot {\text{erf}}^{-1}\left(0.8\right)\sqrt{2}\sigma \end{array}$ (45)

Again, $D_{\text{50}}$ and $z^{\left(c\right)}$ are always on the same side of $x_0$. Next we combine (37) for $x_0<z^{\left(c\right)}$ and (41) for $x_0>z^{\left(c\right)}$ to write out a unified injury probability vs x.

$p\left(x\right)=\frac{1}{2}+\frac{1}{2}\text{erf}\text{}\left(\frac{\text{sign}\left(z^{\left(c\right)}-x_0\right)}{\sqrt{2}\sigma }ln\left(\frac{x-x_0}{\left(z^{\left(c\right)}-x_0\right){\text{e}}^{\mu }}\right)\right)$ (46)

To specify the unified injury function in terms of shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$, we express all quantities in (46) using only $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ and x.

$\frac{\text{sign}\left(z^{\left(c\right)}-x_0\right)}{\sqrt{2}\sigma }=\frac{{\text{erf}}^{-1}\left(0.8\right)}{\gamma }$

$\left(z^{\left(c\right)}-x_0\right){\text{e}}^{\mu }=\left(D_{\text{50}}-x_0\right)=\frac{W}{2sinh(\gamma )}$

$x-x_0=\left(x-D_{\text{50}}\right)+\left(D_{\text{50}}-x_0\right)=\left(x-D_{\text{50}}\right)+\frac{W}{2sinh(\gamma )}$

With these expressions, we write the unified injury function as

$p\left(x\right)=\frac{1}{2}+\frac{1}{2}\text{erf}\text{}\left(\frac{{\text{erf}}^{-1}\left(0.8\right)}{\gamma }ln\left(2sinh\left(\gamma \right)\frac{\left(x-D_{\text{50}}\right)}{W}+1\right)\right)\equiv G\left(x\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ (47)

In injury model (47), the 3 shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ can be specified independently of each other. In particular, for small skewness $\gamma \ll 1$, expanding (47) in terms of $\gamma$ reduces it to the symmetric normal-distribution model (18)

$G\left(x\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)=\frac{1}{2}+\frac{1}{2}\text{erf}\text{}(\frac{2{\text{erf}}^{-1}\left(0.8\right)\cdot \left(x-D_{\text{50}}\right)}{W}+O(\gamma ))$

Figure 3 illustrates several injury functions of the form (47), respectively, for positive, zero and negative skewness. All injury functions shown have the same width $W=5$. In the left panel of Figure 3, injury functions are aligned at $D_{50}=16$. This alignment demonstrates that for $\gamma >0$ the left half of injury function is steeper than the right half; for $\gamma <0$ the left half of injury function is flatter than the right half; and for $\gamma =0$ the injury function is symmetric. In the right panel, injury functions are shifted to be aligned at $D_{10}=13.5$ and thus also aligned at $D_{90}=D_{10}+W=18.5$ because they all have the same width

$W=5$. With $D_{10}=13.5$ fixed, the median dose $D_{\text{50}}$ varies with skewness $\gamma$ from $D_{50}=14.1$ at $\gamma =2$, to $D_{50}=16$ at $\gamma =0$, and to $D_{50}=17.9$ at $\gamma =-2$. The alignment of interval $\left[D_{\text{10}}\mathrm{,}{D}_{\text{90}}\right]$ highlights that as $\gamma$ increases from negative to zero to positive, the injury function becomes more concave down.

6. Effect of Input Dose Uncertainty on the Injury Function with Skewness

We study the effect of input dose estimation uncertainty on the dose-injury function with skewness. We use the term “composite injury function” to denote the injury model after the input dose uncertainty has been incorporated into the model. In general, the composite injury function will be somewhat different from the 3-parameter function form (47) we derived in the previous section. We calculate the three shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ of the composite injury function. Then we explore approximating the composite injury function using function form (47). We examine the difference between the composite injury function and model (47) with the same shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$. If the approximation error is small, then the 3-parameter function form (47) is approximately invariant with respect to input dose uncertainty, and it serves as an adequate framework for accommodating uncertainty in estimating the input dose. Furthermore, framework (47) provides a mechanism of mapping the injury function for one particular dose propagation uncertainty to that for a different uncertainty. Using this mechanism, we can construct an injury model for a target population in application, based on measured injury data for a test population in experiments.

We start with a function of injury probability vs true input dose that is exactly of form (47) specified by 3 shape parameters $\left(D_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$ :

$p_0\left(x^{\left(true\right)}\right)=G\left(x^{\left(true\right)}\mathrm{<mo>\; ;\; </mo>}{D}_{\text{50}}\mathrm{,}W\mathrm{,}\gamma \right)$

We consider the situation where the true input dose $x^{\left(true\right)}$ is not measurable. Instead, an estimated input dose, x, is obtained as an approximation for $x^{\left(true\right)}$. We assume

• the difference $\left(x^{\left(true\right)}-x\right)$ is a normal random variable, and

• the difference $\left(x^{\left(true\right)}-x\right)$ is independent of x.

We assess the injury probability as a function of the estimated input dose x. For each fixed value of x, the corresponding $x^{\left(true\right)}$ is a normal random variable: $x^{\left(true\right)}=x+\mu +\sigma \epsilon$ where $\epsilon ~N\left(\mathrm{0,1}\right)$. The composite injury function, $p_{\sigma }\left(x\right)$, representing the injury probability at estimated input dose x, is a Gaussian weighted average of $p_0\left(x^{\left(true\right)}\right)$ :

$\begin{array}{c}{p}_{\sigma }\left(x\right)=E\left(G\left(x+\mu +\sigma \epsilon ;D_{50},W,\gamma \right)\right)\\ ={\displaystyle {\int }_{-\infty }^{\infty }}G\left(x+\mu +s;D_{50},W,\gamma \right)\frac{1}{\sqrt{2\text{π}{\sigma }^2}}\exp \left(\frac{-s^2}{2{\sigma }^2}\right)\text{d}s\end{array}$ (48)

When injury function $G(\cdot )$ has non-zero skewness, the Gaussian weighted average of $G(\cdot )$ on the right hand side of (48) does not have a simple analytical expression. We use numerical integration to calculate the composite injury function $p_{\sigma }\left(x\right)$ and calculate its shape parameters $\left(D_{\text{50}}\left(\sigma \right)\mathrm{,}W\left(\sigma \right)\mathrm{,}\gamma \left(\sigma \right)\right)$. We examine numerically if $p_{\sigma }\left(x\right)$ is still well described by function form (47) with $p_{\sigma }\left(x\right)$ ’s shape parameters $\left(D_{\text{50}}\left(\sigma \right)\mathrm{,}W\left(\sigma \right)\mathrm{,}\gamma \left(\sigma \right)\right)$.

In our numerical study, $p_0\mathrm{<mo\; stretchy="false">\; (\; </mo>}{x}^{\mathrm{<mo\; stretchy="false">\; (\; </mo>}true\mathrm{<mo\; stretchy="false">\; )\; </mo>}}\mathrm{<mo\; stretchy="false">\; )\; </mo>}$, the injury probability vs the true input dose before input dose uncertainty is incorporated, has function form (47) and is specified by shape parameters $D_{50}=16$, $W=5$, and $\gamma =0.8$. We consider input dose uncertainty of normal distribution with $\mu =0$ (mean) and various values of $\sigma$ (standard deviation). The composite injury function, $p_{\sigma }\left(x\right)$, contains the effect of input dose uncertainty, showing injury probability vs estimated input dose x. Figure 4 examines the composite injury function $p_{\sigma }\left(x\right)$ for $\sigma$ between 0 and 3.

The left panel of Figure 4 shows the injury probability vs the estimated input dose x, respectively, for $\sigma =0,1,2,3$. The most pronounced effect of input dose uncertainty is to spread out the injury function and increase the width. We examine the trend of shape parameters $\left(D_{\text{50}}\left(\sigma \right)\mathrm{,}W\left(\sigma \right)\mathrm{,}\gamma \left(\sigma \right)\right)$ when the input dose uncertainty $\sigma$ is added and increased. The right panel shows $\left(D_{\text{50}}\left(\sigma \right)\mathrm{,}W\left(\sigma \right)\mathrm{,}\gamma \left(\sigma \right)\right)$ vs $\sigma$, of the composite injury function. As the input dose uncertainty $\sigma$ increases, both the median injury dose $D_{50}\left(\sigma \right)$ and the width $W\left(\sigma \right)$ increase monotonically, with W increasing more prominently than $D_{50}$. At the same time, when $\sigma$ increases, the asymmetry of injury function is smoothed out by the Gaussian noise and as a result, the skewness $\gamma \left(\sigma \right)$ decreases. The change in median injury dose $D_{50}\left(\sigma \right)$ is attributed to the presence of skewness: the median injury dose increases (moves toward the right) when an injury function with positive skewness is smoothed out by a Gaussian noise. Conversely, the median injury dose decreases (moves toward the

left) when an injury function with negative skewness is smoothed out. The movement of the median injury dose is caused by smoothing an asymmetric function (see Figure 3 for the general shape of injury functions with positive, zero, and negative skewness). For an injury function of zero skewness, $D_{50}\left(\sigma \right)$ is invariant with respect to $\sigma$ when the injury function is smoothed out by a Gaussian noise.

Next we examine whether or not the composite injury functions for ${\sigma }^2>0$ shown in Figure 4 are still approximately described by model (47). Figure 5 compares the composite injury function and the approximation using function form (47) with shape parameters $\left(D_{\text{50}}\left(\sigma \right)\mathrm{,}W\left(\sigma \right)\mathrm{,}\gamma \left(\sigma \right)\right)$ of the composite injury function. The left panel of Figure 5 compares the composite injury function for $\sigma =1$ and its approximation. The two functions are barely distinguishable from each other. To quantitatively examine the error of approximation, in the right panel we plot the difference between the composite injury function and its approximation. For all values of $\sigma$ examined, the maximum error in approximation is less than 0.01 (1%). The results demonstrate that function form (47) specified by 3 independent shape parameters $\left(D_{50}\mathrm{,}W\mathrm{,}\gamma \right)$ is an adequate model for quantitatively describing general injury functions with skewness.

With the framework of function form (47) and mapping transformation (48), we can filter out the effect of input dose uncertainty in measured injury data. Suppose we are given a measured injury function, $p_{{\sigma }_1}\left(x\right)$, of form (47) for a particular population with input dose uncertainty ${\sigma }_1$. We use transformation (48) to map it back to $p_0\left(x^{\left(true\right)}\right)$, the injury function for the case of zero input dose uncertainty ( $\sigma =0$ ). From there, we can apply the mapping transformation again to predict the injury model for another population with input dose uncertainty ${\sigma }_2$. There is no simple analytical expression for the mapping

transformation. Both the forward and backward mappings need to be implemented numerically. The detailed numerical procedure will be discussed in a subsequent study.

7. Concluding Remarks

We considered injury models in the framework of dose propagation uncertainty. The mathematical formulation is based on that the binary injury outcome is completely determined by the target dose at the active site and the critical threshold. The randomness in the occurrence of injury at a given input dose is attributed to the dose propagation uncertainty from input dose to target dose. The normal distribution model describes the situation where the dose propagation uncertainty is normally distributed. We interpreted the widely used logistic model as a good approximation to the normal distribution model, and thus, interpreted it approximately as a consequence of normally distributed dose propagation uncertainty. In many applications, the input dose is not directly measurable. Instead, an estimated input dose is calculated via computer simulations from measured quantities using representative median parameter values of the general population. In many practical situations, injury models are constructed in the form of injury probability vs estimated input dose. The discrepancy between the estimated input dose and the true input dose can be viewed as an uncertainty in the input dose. With the interpretation of dose propagation uncertainty, the input dose uncertainty is conveniently incorporated into the injury model. The framework of dose propagation uncertainty provides a mechanism of extending an injury function established on a test population to predict the injury model for a different population in application. Both the logistic model and the normal distribution model are specified by two shape parameters: the median injury dose and the 10 - 90 percentile width. The mapping between the injury functions of two populations has a simple analytical form of updating the two shape parameters. Both the logistic model and the normal distribution model are symmetric around the median injury dose and have no skewness. To accommodate injury functions with skewness, we studied dose propagation uncertainties of shifted log normal distribution with shift as a parameter. Based on the shifted log normal model, we developed a function form for injury probability vs input dose that is specified by three shape parameters: median injury dose, the width, and the skewness. The proposed function form allows the three shape parameters to be set independent of each other. In particular, the proposed function form is capable of accommodating arbitrary skewness, positive or negative. In addition, we showed numerically that the proposed 3-parameter function form is approximately invariant with respect to additions or changes in input dose uncertainty. Therefore, the 3-parameter function form serves as a broad framework for modeling input dose uncertainty and modeling injury function skewness at the same time. This broad framework allows us to map injury function with skewness from a test population to a different population in applications.

Disclaimer and Acknowledgements

The authors thank C. Kramer and J. Swallow of Institute for Defense Analysis (IDA) for bringing the problem to their attention, and thank the Joint Non-Lethal Weapons Directorate of U.S. 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 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. (2018) Dose-Injury Relation as a Model for Uncertainty Propagation from Input Dose to Target Dose. American Journal of Operations Research, 8, 360-385. https://doi.org/10.4236/ajor.2018.85021

References