Development and Validation of a Predictive Bone Fracture Risk Model for Astronauts

Abstract

There are still many unknowns in the physiological response of human beings to space, but compelling evidence indicates that accelerated bone loss will be a consequence of long-duration spaceflight. Lacking phenomenological data on fracture risk in space, we have developed a predictive tool based on biomechanical and bone loading models at any gravitational level of interest. The tool is a statistical model that forecasts fracture risk, bounds the associated uncertainties, and performs sensitivity analysis. In this paper, we focused on events that represent severe consequences for an exploration mission, specifically that of spinal fracture resulting from a routine task (lifting a heavy object up to 60 kg), or a spinal, femoral or wrist fracture due to an accidental fall or an intentional jump from 1 to 2 m. We validated the biomechanical and bone fracture models against terrestrial studies of ground reaction forces, skeletal loading, fracture risk, and fracture incidence. Finally, we predicted fracture risk associated with reference missions to the moon and Mars that represented crew activities on the surface. Fracture was much more likely on Mars due to compromised bone integrity. No statistically significant gender-dependent differences emerged. Wrist fracture was the most likely type of fracture, followed by spinal and hip fracture.

Appendices

Appendix A: Calculation of Fracture Probability from Fracture Risk Index

Davidson et al. 26 provided a means of linking estimates of fracture risk index, FRI, to fracture probability, p _fx, through the use of logistic regression to compare the binary condition of actual fractures to non-fractures, with reference to a set of appropriate controls and specified loading events. The logistic regression used in this study results in a sigmoidal curve represented by:

$$ p_{\text{fx}} = {\frac{1}{{\left( {1 + { \exp }\left( { - 1*\left( {{\text{FRI}} - \mu } \right)*\phi } \right)} \right)}}} $$

(A1)

where μ is the position factor of the curve (the value of FRI where the probability is 0.50) and ϕ is the slope factor (a measure of the steepness of the curve). Ideally, the development of relevant values of μ and ϕ would be derived from data of actual skeletal loading events, in which conditions and outcomes (fracture or no fracture), are well documented. Such data are needed to ensure that the FRI to probability relation is consistent with the methods used in the prediction of FRI. Unfortunately, data of this type are lacking in the literature. Digitizing and post-processing Davidson et al.’s data, derived from radial arm fractures in children caused by falls from playground equipment, reveal values of μ and ϕ in the range of 0.58 to 0.6 and 7.5 to 12, respectively. Although not directly applicable to the case of astronaut fractures, these values provide a rational first comparison in evaluating a translation function linking the FRI to estimated fracture probability for astronauts.

To develop the translation function for astronauts without direct loading and fracture data, we must make several assumptions concerning the sigmoidal function parameters, and determine an acceptable fracture threshold range. Several articles in the literature suggest that, at the proximal femur, the fracture threshold (the range when the probability of fracture is not negligible) occurs when the applied load, AL, equals fracture load, FL, within ±1σ of the bone strength17 ^, 58 or

$$ 0 < p_{\text{fx}} < 1\quad {\text{when}}\; {\text{FL}} - \sigma_{\text{FL}} < {\text{AL}} < {\text{FL}} + \sigma_{\text{FL}} $$

where p _fx is the probability of fracture (a number between 0 and 1) and σ is the standard deviation. For our purposes, this threshold would formally translate to:

$$ \begin{aligned} & 0 < p_{\text{fx}} < 1\quad {\text{when}}\; 1 - \sigma_{\text{FRI}} < {\text{FRI}} < 1 + \sigma_{\text{FRI}} \\ & \sigma_{{{\text{FRI}} = 1}} = \left[ {\left( {\sigma_{\text{AL}} {\frac{1}{\text{FL}}}} \right)^{2} + \left( {\sigma_{\text{FL}} {\frac{\text{AL}}{{2{\text{FL}}^{2} }}}} \right)^{2} } \right]^{1/2} \\ \end{aligned} $$

(A2)

The threshold of fracture has a mean value of FRI = 1, and the uncertainty of the threshold value is dependent on the standard deviations of AL and FL. Making use of the fact that at FRI = 1, AL = FL, and using the uncertainty of the loading condition and bone strengths, we estimated that 0.135 < σ_FRI=1 < 0.22. However, testing of the model using the parameter space imposed by the astronaut data resulted in a standard deviation of the FRI estimates ranging from 0.29 to 0.67. This is due, in part, to the inclusion of other parameter uncertainties and the generally higher mean values of FRI that are calculated. Therefore, to be inclusive of all available data and corresponding uncertainties, the formal estimate of σ at FRI = 1 was adjusted to range from 0.135 to 0.67 in all of the calculations presented in this paper.

To be reasonably confident that the estimated range of μ and ϕ includes the “true” range of μ and ϕ, the following assumptions are made based on the guidance found in the literature:

• 1. The mean FRI threshold value is 1

• 2. At FRI = 1 + σ _FRI=1, p _fx = 0.95 represents the upper limit of the fracture threshold

• 3. At FRI = 1 − σ _FRI=1, p _fx = 0.05 represents the lower limit of the fracture threshold

In order to use these assumptions to generate a range of possible ϕ values, the sigmoid equation is first solved for ϕ. This gives:

$$ \phi = - 1*{\frac{{{ \ln }\left( {{\frac{1}{{p_{\text{fx}} }}} - 1} \right)}}{{\left( {{\text{FRI}} \pm \sigma_{{{\text{FRI}} = 1}} - \mu } \right)}}} $$

(A3)

Under our assumptions that the threshold value for FRI = 1, and assuming μ = 1, then the range for ϕ in the threshold region is found to be dependent only on σ _FRI=1. For the estimated range of σ _FRI=1, the range of estimated ϕ is found to be 4.4 < θ < 22. This approach produces a range on ϕ that is inclusive of the values estimated by Davidson et al. 26

We can again turn to the sigmoid equation to assist in estimating the range of values for μ. Rearranging the sigmoid equation, μ can be calculated as:

$$ \mu = {\frac{{\ln \left( {{\frac{1}{{p_{\text{fx}} }}} - 1} \right)}}{\phi }} + \left( {{\text{FRI}} \pm \sigma_{{{\text{FRI}} = 1}} } \right) $$

(A4)

Assuming that the mean FRI = 1 and using the previous calculation, the mean value for ϕ is evaluated as 13.2. Applying the change in p _fx and σ _FRI=1 over the combination of their respective ranges, the range for μ is found to be 0.55 < μ < 1.45, which is inclusive of the Davidson et al. 26 estimates. The limits on the family of sigmoid curves, as described over the range of μ and ϕ, are illustrated in Fig. A.1. Although additional data would reduce the uncertainty in these calculations, the approach described here reasonably bounds the available knowledge base, makes use of expert opinion in the literature, and provides a consistent means of representing the uncertainty in estimating fracture risk probability.

Figure A.1

Sigmoid curve family representing the estimated fracture threshold translation function for transforming FRI to probability of fracture. Solid line represents mean value of μ and ϕ, while upper and lower fracture thresholds are shown in dotted lines at specified combinations of (μ, ϕ)

To implement the sigmoid translation function for FRI to fracture probability in the probabilistic model for fracture:

• 1. Uniform probability distributions were established over the range of μ and ϕ.

• 2. At each trial when an FRI value was calculated in the model, the μ and ϕ distributions were randomly sampled.

• 3. Using the randomly sampled combined values of μ and ϕ, a translation function using the general form of the sigmoid equation was used to estimate the probability of fracture from the FRI value.

Uniform probability distributions are sometimes referred to as “maximum uncertainty” distributions since no value is more probable to occur than any other value within the limits of the distribution. We chose this type of distribution for estimating μ and ϕ due to the lack of knowledge regarding the shape of a representative probability density function and to formally represent the large (epistemic) uncertainty in this area of the model.

Appendix B: Equations for Biomechanical Model of the Hip and Lumbar Spine

The biomechanical models for dynamic loading shown in Figs. 3b–3d lead to sets of first-order differential equations with unknown time-dependent displacements, x, and velocities, $\dot{x}.$ For the case of the hip, the equation of motion is:

$$ m_{\text{H}} \ddot{x}_{\text{H}} + b_{\text{H}} \dot{x}_{\text{H}} + k_{\text{H}} x_{\text{H}} = m_{\text{F}} g $$

(B1)

where the subscript H refers to the hip, m is the mass of the body, $\ddot{x}$ is the acceleration, k is the spring constant, and b is the damping coefficient. In all cases, the initial displacement is set to zero, and the impact velocity is set to$\sqrt {2gh} ,$ where h is the distance fallen. In this case, h is defined as the distance between the ground and the hip.

The equations can be more compactly written as matrix equations, in which the unknown displacement x _H is set to x ₁ and the unknown velocity, $\dot{x}_{1}$ to x ₂

$$ \left[ {\begin{array}{*{20}c} {\dot{x}_{1} } \\ {\dot{x}_{2} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 0 & 1 \\ { - b/m} & { - k/m} \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {x_{1} } \\ {x_{2} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} 0 \\ g \\ \end{array} } \right] $$

(B2)

The initial conditions are zero displacement and impact velocity as set by the local gravitational acceleration:

$$ \left[ {\begin{array}{*{20}c} {x_{1} (0)} \\ {x_{2} (0)} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 0 \\ {\sqrt {2gh} } \\ \end{array} } \right] $$

(B3)

Once the system of equations is solved, the ground reaction force, F _GR, is found by summing the damping and spring forces that are applied between the body and the ground:

$$ F_{\text{GR}} = bx_{2} + kx_{1} $$

(B4)

For the femoral fracture model, F _GR multiplied by attenuation factors for the fall orientation and active response is the applied load to the bone. In the biomechanical models for the lumbar spine and the wrist, the body is represented as a series of linked masses connected with linear springs and dampers. For the lumbar spine:

$$ \begin{aligned} & m_{\text{F}} \ddot{x}_{\text{F}} + \dot{x}_{\text{F}} b_{\text{G}} + x_{\text{F}} k_{\text{G}} + (x_{\text{F}} - x_{\text{PL}} )k_{\text{L}} = m_{\text{F}} g \\ & m_{\text{PL}} \ddot{x}_{\text{PL}} + (\dot{x}_{\text{PL}} - \dot{x}_{\text{HAT}} )b_{\text{LS}} + (x_{\text{PL}} - x_{\text{HAT}} )k_{\text{LS}} + (x_{\text{PL}} - x_{\text{F}} )k_{\text{L}} = m_{\text{PL}} g \\ & m_{\text{HAT}} \ddot{x}_{\text{HAT}} + (\dot{x}_{\text{HAT}} - \dot{x}_{\text{PL}} )b_{\text{LS}} + (x_{\text{HAT}} - x_{\text{PL}} )k_{\text{LS}} = m_{\text{HAT}} g \\ \end{aligned} $$

(B5)

When we perform a similar re-definition of the matrix variables by:

$$ \begin{aligned} x_{1} = & x_{\text{F}} \\ x_{2} = & \dot{x}_{1} \\ x_{3} = & x_{\text{PL}} \\ x_{4} = & \dot{x}_{3} \\ x_{5} = & x_{\text{HAT}} \\ x_{6} = & \dot{x}_{5} \\ \end{aligned} $$

(B6)

then the following matrix equation can be found:

$$ \left[ {\begin{array}{*{20}c} {\dot{x}_{1} } \\ {\dot{x}{}_{2}} \\ {\dot{x}_{3} } \\ {\dot{x}_{4} } \\ {\dot{x}_{5} } \\ {\dot{x}_{6} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 0 & 1 & 0 & 0 & 0 & 0 \\ { - (k_{\text{G}} + k_{\text{L}} )/m_{\text{F}} } & { - b_{\text{G}} /m_{\text{F}} } & {k_{\text{L}} /m_{\text{F}} } & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ {k_{\text{L}} /m_{\text{PL}} } & 0 & { - (k_{\text{L}} + k_{\text{LS}} )/m_{\text{PL}} } & { - b_{\text{LS}} /m_{\text{PL}} } & {k_{\text{LS}} /m_{\text{PL}} } & {b_{\text{LS}} /m_{\text{PL}} } \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & {k_{\text{LS}} /m_{\text{HAT}} } & {b_{\text{LS}} /m_{\text{HAT}} } & { - k_{\text{LS}} /m_{\text{HAT}} } & { - b_{\text{LS}} /m_{\text{HAT}} } \\ \end{array} } \right] \cdot \left[ {\begin{array}{*{20}c} {x_{1} } \\ {x_{2} } \\ {x_{3} } \\ {x_{4} } \\ {x_{5} } \\ {x_{6} } \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} 0 \\ g \\ 0 \\ g \\ 0 \\ g \\ \end{array} } \right] $$

(B7)

with initial conditions:

$$ \left[ {\begin{array}{*{20}c} {x_{1} (0)} \\ {x{}_{2}(0)} \\ {x_{3} (0)} \\ {x_{4} (0)} \\ {x_{5} (0)} \\ {x_{6} (0)} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 0 \\ {\sqrt {2gh} } \\ 0 \\ {\sqrt {2gh} } \\ 0 \\ {\sqrt {2gh} } \\ \end{array} } \right] $$

(B8)

Similar to the hip model, the ground reaction force, F _GR, and the force on the lumbar spine, F _LS, is given by:

$$ F_{\text{GR}} = b_{\text{G}} x_{2} + k_{\text{G}} x_{1} $$

$$ F_{\text{LS}} = b_{\text{LS}} (x_{6} - x_{4} ) + k_{\text{LS}} (x{}_{5} - x_{3} ) $$

The force on the lumbar spine is equivalent to the applied load to the spine, which was then used in the spinal fracture model.