A non-linear multiaxial fatigue damage model for the cervical intervertebral disc annulus

Abstract

A significant portion of the military population develops severe neck pain in the course of their duties. It has been hypothesized that neck pain is a consequence of accelerated degeneration of the intervertebral discs in the cervical spine, but more occupational and mechanistic-based tools and research are needed to positively confirm the link between neck pain and accelerated disc degeneration. Heavy head-supported mass including helmets and accessories worn by military personnel may subject the intervertebral discs of the cervical spine to complex cyclic loading profiles. In addition, some military operational travel which includes riding on high speed planing boats has also been reported to result in high magnitude cyclic loading on cervical spine discs. In this article, we present a methodology to computationally predict fatigue damage to cervical intervertebral discs over extended periods of time, by integrating kinematics-based biomechanical models with a continuum damage mechanics-based theory of disc degeneration. Through this computational approach, we can gain insights into the relationship between these military activities and possible accelerated fatigue degeneration of cervical intervertebral discs and provide a quantitative prediction tool for decade-long time ranges. The four significant improvements this computational framework adds to the area of modeling intervertebral disc degeneration are the following: (a) it addresses the non-linear nature of fatigue damage evolution, (b) it includes the effect of aging and damage recovery to accurately simulate biological phenomena, (c) it computes fatigue damage taking into account the multiaxial stress state in the disc, and (d) it correlates the computational damage parameter with established clinical grading systems for disc degeneration.

Keywords

Computational biomechanics cervical spine intervertebral disc annulus degeneration multiaxial fatigue

Introduction

The development of severe neck pain among personnel of the military is becoming a cause for concern.^1–3 Van den Oord and Veerle De Loose¹ conducted a study on military helicopter pilots and concluded that heavy helmets and associated head-mounted displays commonly used by these pilots may play a significant role in the development of neck pain. Moreover, this heavy head-supported equipment is believed to affect the biomechanics of the head–neck complex, giving rise to complex high amplitude cyclic loads on the cervical spine at the vertebra-disc level.^4,5 For a different military population, Ensign et al.² performed a survey of self-reported injuries among special boat operators who operate small boats in the open ocean and are thus exposed to high amplitude vibration forces on the head–neck complex. He concluded that an increased level of neck pain was prevalent among them. Furthermore, Cohen et al.³ performed an observational study on soldiers who were medically evacuated out of combat operations due to neck pain, concluding that these evacuations were associated with a low rate of return-to-duty.

It has been hypothesized that the biological cause for neck pain is the accelerated degeneration of the intervertebral discs of the cervical spine.^6–8 Podichetty⁶ argued that neck and back pain are manifestations of degenerative disc disorders of the spine and hypothesized that the degenerative cascade ultimately results in extensive structural defects in the discs and reduction of motion compliance of the spinal segments. These structural defects in the discs are manifested as fracture of the cartilaginous endplate, decrease in cell density, and formation of microstructural clefts and tears in the disc annulus material.⁹ Buckwalter⁷ postulated that one of the degenerative changes in cervical discs that contributes to spine stiffness and neck pain is fatigue failure of the disc annulus material. The main focus of this article is to computationally predict the evolution of fatigue damage in the annulus region of a cervical spine intervertebral disc, subjected to complex cyclic loading over long periods of time spanning several years. Previous studies by researchers have shown that the percentage of denatured collagen in the annulus increased with intervertebral disc degeneration.^7,9,10 Moreover, the clinical relevance of endplate and vertebral fractures as generators of pain or involvement in a degenerative process has not been established.¹¹ Hence, as an indicator of neck pain caused by disc degeneration, we attempt to predict the fatigue damage evolution specifically in the disc annulus.

Several studies have been carried out by researchers to investigate the fatigue failure properties of intervertebral discs, such as the number of load cycles to failure for various stress amplitudes and the variation in compressional stiffness of the disc as a result of fatigue. Yu et al.¹² conducted a study on the geometric and morphological changes of porcine intervertebral discs subjected to fatigue loading. Researchers have established that though porcine vertebrae are anatomically smaller than human vertebrae, they possess a similar ligamentous structure and facet joint orientation. Stiffness values are similar, and comparable injuries resulted from compressive and shear loads.¹³ Hence, porcine spines can be a useful model for studying human spinal injury. Using magnetic resonance imaging, they observed that fatigue damage in the discs manifested itself as a significant change in the shape and area of the nucleus pulposus region, bulging of the anterior and posterior longitudinal ligaments and dehydration of the annulus fibrosus region. It was concluded that disc degeneration and fatigue failure of the disc are strongly related. Snyder¹⁴ studied the biomechanical properties of functional spine units (FSUs) of the cervical spine subjected to different amplitudes of tension-compression fatigue loading. An FSU consists of two adjacent vertebrae, the intervertebral disc and all adjoining ligaments between the vertebrae. Using a dual ultrasound system, he recorded the variation in compressional stiffness of the FSUs caused by fatigue and concluded that discs that had lost their hydration and structural integrity as a result of fatigue loading had a lower motion compliance than healthy discs. Green et al.¹⁵ conducted fatigue failure tests on the annulus fibrosus region of human cadaveric lumbar intervertebral discs. They identified the failure point as the instant at which the annular fibers detached from the cartilaginous endplate of the underlying vertebra. Under fully reversed cyclic loading, they recorded the number of loading cycles to failure for a range of alternating stress amplitudes.

These experimental results have been used by a few researchers to develop computational models for fatigue damage prediction in intervertebral discs. In order to quantify the state of fatigue damage in the disc, the concept of structural damage within a material as postulated by Kachanov¹⁶ has been implemented. Kachanov proposed the use of an internal state variable of the material which represents the level of structural damage at a given point in the material. This numerical value of this state variable can range from 0 to 1, corresponding to completely undamaged and completely damaged material, respectively. Qasim et al.¹⁷ ran computational fatigue simulations on a finite element model of a lumbar spinal segment and predicted damage evolution in the disc annulus using continuum damage mechanics and Kachanov’s concept of the damage variable. The location of damage initiation in the disc was hypothesized to be at the integration point where the damage variable first crossed a threshold value of 0.9. Makwana et al.¹⁸ employed a finite element model of a cervical spinal segment to study the pattern of damage initiation and progression in the disc annulus. They modeled the disc as a poroelastic material, which enabled further quantification of damage in terms of water loss from the disc and change in disc height.

It is important to note that these past modeling efforts utilize a linear damage incrementation scheme given by the Palmgren–Miner hypothesis,¹⁹ which has been widely used throughout literature for convenience and simplicity. However, linear damage incrementation schemes have some major shortcomings. First, linear models do not capture the effect of the order of application of load cycles with different alternating stress amplitudes.²⁰ Moreover, linear models neglect cycle ratios that are below the fatigue limit. This can cause under-prediction of damage, as prior loading can result in lowering of the fatigue limit.²¹ These shortcomings of linear damage incrementation schemes can result in discrepancies of up to an order of magnitude between predicted and experimental fatigue life.²² In order to address these limitations of linear models, Marco and Starkey²³ proposed the first non-linear load-dependent damage evolution theory. As an extension of this concept of a non-linear fatigue damage model, Lemaitre and Chaboche²⁴ proposed an analytical differential equation for damage evolution of the form

\frac{δ D}{δ N} = f (D, σ_{m a x}, σ_{m e a n})

(1)

where the increment in damage per cycle $δ D / δ N$ is a function of the current value of damage D and the maximum and mean stress of the load cycle. Researchers have successfully utilized this model to study damage evolution in biological tissue such as human bones and tendons.^25–27

Viscoelastic damage recovery of biological tissue has been studied and implemented in computational models for musculoskeletal tissues such as bones and tendons.^26,28 Daily recovery of the accumulated damage during rest periods has also been observed in intervertebral disc tissue,^7,29–31 though to a lesser degree, as compared to musculoskeletal tissue. Johannessen et al.³¹ elucidated that though the disc is an avascular tissue characterized by an absence of active biological processes that heal fatigue damage to the disc, a form of passive recovery can be seen in the disc. They observed that the compressive loads experienced by the disc were greatly reduced during periods of supine rest, and passive recovery occurs in the form of re-imbibition of fluids into the disc. This fluid flux results in the restoration of intradiscal hydrostatic pressure and disc height, which translates to restoration of the intervertebral disc mechanics prior to cyclic loading. Moreover, it has been observed that human intervertebral discs undergo multifactorial biochemical and morphologic degenerative changes during the process of aging, such as modification of proteins in the disc matrix and accumulation of degraded disc matrix molecules.^6,7 Therefore, it is important to address the effects of these biological aging processes and recovery mechanisms when modeling disc fatigue damage over long periods of time, spanning several years. These effects have not been explicitly captured in previous models of disc fatigue damage.

Previous efforts toward modeling fatigue damage to the disc consider only a uniaxial stress state in the disc material. In these models, an assumption is made that the only non-zero component of the stress tensor of the disc material is the normal compressive stress acting along the inferior–superior direction. However, the stress state in the disc matrix is multiaxial in reality, with non-zero shear stresses acting along the disc surface. Therefore, it is important to consider the effects of this multiaxial stress state while estimating fatigue damage to the disc. Ion et al.³² elucidate the use of an equivalent stress measure known as the signed Von-Mises equivalent stress for multiaxial fatigue calculation in a material. This signed Von-Mises stress measure can be used instead of just the compressive normal stress in order to estimate fatigue damage to the disc, taking into account the multiaxial stress state in the disc.

We propose to develop a novel computational model to predict damage in the annulus of cervical spine intervertebral discs, due to cyclic loading over extended periods of time. This model will incorporate the following four improvements over previous disc fatigue modeling efforts: (a) it addresses the non-linear nature of fatigue damage evolution, (b) it includes the effect of aging and damage recovery to accurately simulate biological phenomena, (c) it computes fatigue damage taking into account the multiaxial stress state in the disc, and (d) the computational damage parameter can be correlated with established clinical grading systems for disc degeneration.

Cyclic loading histories on the disc which are relevant to various military activities can be given as input to the damage model, in order to obtain an estimate of disc fatigue damage in military personnel over extended periods of time. In order to validate this model, the value of the damage variable reported by the model needs to be correlated to disc degeneration data from experiments performed on real soldiers over long periods of time.

Overview of the damage model framework

Intervertebral disc degeneration can arise from several different factors such as mechanical loading, diseases, and aging. Consider the following equation proposed by Nash,³³ describing the total damage to biological tissue

D_{tot} = D_{m} + D_{a} + D_{d} - H

(2)

where $D_{m}$ is the damage due to mechanical loading, $D_{a}$ is the damage due to aging phenomena, $D_{d}$ is the damage due to disease, and $H$ represents the recovery of damage due to healing processes in intervertebral discs. Moreover, Nash³³ postulated that the damage due to aging phenomena cannot be recovered, which leads to the conclusion that only the damage due to mechanical loading and disease can be recovered. Hansen³⁴ observed that intervertebral disc disease manifests itself as herniation of the nucleus pulposus of the disc through the annular fibers, extrusion of the nucleus material into the spinal canal, and protrusion of the annulus caused by shifting of the central nucleus material. The damage model presented in this article does not consider damage due to disease, because it is not likely that individuals affected by spine-related diseases will be allowed to serve in the military for a long period of time. Therefore, in this article, the term $H$ representing damage recovery has been incorporated into the expression for the mechanical damage $D_{m}$ . This leads to the following simplified version of equation (2), which describes the damage model in this article

D_{tot} = D_{m} + D_{a}

(3)

A brief overview of the computational framework adopted for the damage model presented in this article has been shown in Figure 1. The damage prediction process may be broadly classified into three categories: model development, model calibration, and model application. In model development, the first step involves determining the parameters in the equation for mechanical damage, using the disc fatigue failure data provided by Green et al.¹⁵ It must be noted that there is a lack of fatigue failure studies in the literature pertaining to the cervical disc annulus. Therefore, as a reasonable approximation, we adopted the lumbar annulus fatigue data of Green et al.¹⁵ to determine the mechanical damage parameters of the model. The next step involves determining the parameters in the equation for damage due to aging, using disc degeneration data made available by Thompson et al.³⁵ The sum of the mechanical damage and the damage due to aging represents the total damage to the intervertebral disc. In model calibration, the damage model is run for the case of neck loading due to normal walking without head-supported mass, using the determined model parameters. The total damage predicted by the model is compared to disc degeneration data over several decades for the general population from multiple sources in the literature.^36–38 This comparison is used to finally calibrate parameters of the mechanical damage equation that could not be determined from the disc fatigue failure data provided by Green et al.¹⁵ Finally, in model application, the calibrated damage model is used to predict damage evolution in soldiers due to neck loading arising from walking with head-supported mass and also from riding a fast boat. Each of these three broad categories will be discussed in detail in the following sections.

Figure 1.

Overview of the computational framework of the damage model: model development, model calibration, and model application.

Model development

Equation for mechanical damage evolution

It is well established that anatomically, the annulus fibrosus of the intervertebral disc is a layered structure, with each layer enforced by nearly aligned collagen fibers. From the point of view of continuum mechanics, the annulus fibrosus should be modeled as a composite with strongly anisotropic, non-linear material behavior. However, the damage model presented in this article is independent of such a constitutive material description of the annulus fibrosus, as will become clear in the following sections. The equation used to describe the mechanical damage evolution in the intervertebral disc annulus was adopted from the following non-linear continuous fatigue damage model proposed by Lemaitre and Chaboche²⁴ and Chaboche and Lesne³⁹

\frac{δ D_{m}}{δ N} = {D_{m}}^{α (t)} {[\frac{σ_{a}}{M_{0} (t) (1 - b (t) σ_{mean})}]}^{β (t)}

(4)

α (t) = 1 - a (t) \frac{σ_{a}}{σ_{u} - σ_{\max}}

(5)

where $D_{m}$ is the mechanical damage $(0 \leq D_{m} \leq 1)$ . $D_{m} = 0$ represents the undamaged state of the disc and $D_{m} = 1$ represents the damage in a completely degenerated disc. $δ D_{m} / δ N$ is the increment in damage per cycle of loading, $σ_{a}$ is the alternating stress, $σ_{mean}$ is the mean stress, $σ_{\max}$ is the maximum value of stress in the load cycle, and $σ_{u}$ is the ultimate tensile strength of the material. $α (t)$ is a function of $σ_{a}$ , $σ_{u}$ , and $σ_{\max}$ . The symbol 〈〉 is defined as $x = 0$ if $x < 0$ , and $x = x$ if $x \geq 0$ . $M_{0} (t)$ , $a (t)$ , $b (t)$ , and $β (t)$ are material parameters which can be determined from curve fitting using stress-life data from fatigue testing experiments on the disc annulus.¹⁵ In reality, the values of these material parameters may change over time (tissue re-modeling) depending on various factors such as the current damage state of the disc annulus and magnitude of the applied loads. However, there is a lack of relevant data in the literature to establish these relationships over a time period of several years. Therefore, as a simplification, these material parameters were assumed to be constant throughout the damage evolution process.

In order to determine the relationships of these material parameters with time, we would need to obtain multiple sets of S-N curve failure data, by performing fatigue failure tests on the annulus fibrosus of discs of various ages and damage levels. Then, a separate curve fit would need to be done for each of these S-N curve failure datasets in order to establish the variation of the material parameters of the model over time. This can be pursued in future work.

The increment in mechanical damage per cycle in equation (4) depends on the current value of the damage variable. In other words, the mechanical damage variable $D_{m}$ is present on both sides of equation (4), and the exponent $α$ depends on the applied load. This makes equation (4) inseparable in terms of damage and the applied load,²⁴ giving rise to the non-linear nature of fatigue damage evolution under cyclic loading. This is an important feature captured by the model presented in this article, as it has been observed that changing the order of application of cyclic loads of different magnitudes has an effect on the fatigue damage caused to the material,²⁰ and separable equations are not capable of modeling these effects.⁴⁰

Integrating equation (4) as presented by Dattoma et al.⁴¹ produces the following set of equations

N_{f} = \frac{1}{1 - α} {[\frac{σ_{a}}{M_{0} (1 - b σ_{mean})}]}^{- β}

(6)

D_{m} = {(\frac{n}{N_{f}})}^{\frac{1}{1 - α}}

(7)

where $N_{f}$ is the number of load cycles required for complete failure, that is, for the value of the damage variable to increase to 1.

The equations for mechanical damage evolution presented above only apply to cyclic loads having constant stress amplitudes. However, the cyclic load amplitudes on the cervical intervertebral disc during one gait cycle of walking vary continuously with time as shown in Figure 5. Therefore, the equations for damage evolution presented in this section were modified for variable amplitude loading, as elucidated by Dattoma et al.⁴¹ These modified equations are given below

N_{i} = N_{fi} {(\frac{n_{i - 1}^{eff}}{N_{fi - 1}})}^{\frac{(1 - α_{i})}{(1 - α_{i - 1})}}

(8)

n_{i}^{eff} = n_{i} + N_{fi} {(\frac{n_{i - 1}^{eff}}{N_{fi - 1}})}^{\frac{(1 - α_{i})}{(1 - α_{i - 1})}}

(9)

D_{mi} = {(\frac{n_{i}^{eff}}{N_{fi}})}^{\frac{1}{1 - α_{i}}}

(10)

where $N_{fi}$ is the number of cycles to failure at iteration i, which is determined from $σ_{ai}$ , the alternating stress amplitude at iteration i; $n_{i}$ is the actual number of load cycles with alternating stress amplitude $σ_{ai}$ ; $N_{i}$ is the number of cycles with alternating stress amplitude $σ_{ai}$ that would need to be applied to a completely undamaged material, in order to cause the same damage as the damage due to all previous load iterations (1 to i – 1); $n_{i}^{eff}$ represents the effective number of cycles with alternating stress amplitude $σ_{ai}$ that would need to be applied to a completely undamaged material to cause the same damage as the actual damage at the end of load iteration i; and $D_{mi}$ is the effective mechanical damage at the end of load iteration i.

Multiaxial fatigue estimation

The stress state in the matrix of the disc annulus is multiaxial in reality, with non-zero shear stresses acting along the disc surface. Ion et al.³² elucidated the use of an equivalent stress measure known as the signed Von-Mises equivalent stress for multiaxial fatigue calculation in a material. This signed Von-Mises stress measure was used in the model presented in this article in order to estimate fatigue damage to the disc, taking into account the multiaxial stress state in the disc.

The loading configuration of the disc was idealized, assuming that the top surface of the disc is subjected to a compressive normal load in the inferior–superior direction and shear loads in the anterior–posterior and lateral directions. Under this loading configuration, we propose that the Cauchy stress tensor in the disc annulus material at any instant can be represented by

σ = [\begin{matrix} 0 & σ_{xy} & 0 \\ σ_{xy} & σ_{yy} & σ_{yz} \\ 0 & σ_{yz} & 0 \end{matrix}]

(11)

where x represents the anterior–posterior direction, y represents the inferior–superior direction, and z represents the lateral direction. The only non-zero normal stress is $σ_{yy}$ and the non-zero shear stresses are $σ_{xy}$ and $σ_{yz}$ . These assumptions were verified by means of a simple one-element static simulation in ABAQUS. The bottom surface of this element was fixed, and the element was subjected to the loading configuration described above. It was seen that the resulting Cauchy stress tensor in this element agrees with the form described by equation (11) and was uniform throughout the element.

The expression for the equivalent stress measure for multiaxial fatigue calculation known as the signed Von-Mises equivalent stress has been given below, as explained by Ion et al.³²

\begin{matrix} σ_{VM, s} = \frac{sign (σ_{xx} + σ_{yy} + σ_{zz})}{\sqrt{2}} \\ \sqrt{{(σ_{xx} - σ_{yy})}^{2} + {(σ_{yy} - σ_{zz})}^{2} + {(σ_{zz} - σ_{xx})}^{2} + 6 ({σ_{xy}}^{2} + {σ_{yz}}^{2} + {σ_{zx}}^{2})} \end{matrix}

(12)

where $σ_{VM, s}$ is the signed Von-Mises stress measure. For the particular loading configuration described above, equation (12) can be simplified to yield equation (13) shown below. From equation (13), we can see that the sign of the signed Von-Mises stress measure is given by the sign of $σ_{yy}$ , the normal stress in the y-direction

σ_{VM, s} = sign (σ_{yy}) \sqrt{{(σ_{yy})}^{2} + 3 ({σ_{xy}}^{2} + {σ_{yz}}^{2})}

(13)

In the damage model presented in this article, the signed Von-Mises stress measure shown in equation (13) was computed at every instant of the cyclic loading histories. This equivalent stress measure was then used to compute $σ_{a}$ , $σ_{mean}$ , and $σ_{\max}$ for each load iteration, in order to estimate the multiaxial fatigue damage evolution in the disc annulus.

Recovery of mechanical damage

Previous studies have shown that intervertebral discs can recover a considerable amount of the damage caused due to cyclic loading. This recovery typically takes place during a rest period, when the disc is not loaded.^31,42,43 We have proposed a novel method as an effort to capture this effect in our damage model. At the end of each day, we assume that a rest period occurs during sleep when the disc is unloaded and recovers some of the mechanical damage. We have introduced a recovery parameter r that reduces the mechanical damage as shown in equation (14)

D_{mi}^{r} = D_{mi - 1}^{r} + (D_{mi} - D_{mi - 1}^{r}) (1 - r)

(14)

where $D_{mi}$ is the mechanical damage at the end of day i without recovery and $D_{mi}^{r}$ is the mechanical damage at the end of day i after recovery. From equation (14), it can be seen that the recovery parameter r is the fraction of mechanical damage accumulated during one day that is recovered during the daily rest period. From equation (9), it is clear that the effective number of cycles $n_{i}^{eff}$ at each iteration is calculated using the number of load cycles to failure for the previous iteration, $N_{fi - 1}$ . Therefore, the value of $N_{fi}$ needs to be modified as follows to incorporate the effect of recovery in the next iteration

N_{fi} = \frac{n_{i}^{eff}}{{(D_{mi}^{r})}^{\frac{1}{(1 - α_{i})}}}

(15)

Curve fitting to determine mechanical damage parameters

In order to implement all the equations developed so far to estimate mechanical damage evolution in the disc annulus material, the parameters $a$ , $b$ , $β$ , and $M_{0}$ in equations (5) and (6) need to be determined for this particular material. These parameters were determined using the ultimate tensile strength and stress-life data for the disc annulus obtained from the fatigue failure experiments of Green et al.¹⁵ Having also conducted tensile failure tests on vertical slices of lumbar intervertebral discs, Green et al.¹⁵ reported that the ultimate tensile strength of the anterior disc annulus is 3.9 MPa, whereas that of the posterior disc annulus is 8.6 MPa. It must be noted that though our study concerns the cervical disc annulus, this is the most relevant data we could obtain to estimate the ultimate tensile strength of the cervical disc annulus. As we are considering a uniform stress distribution throughout the disc annulus, an average value of 6.25 MPa was chosen for the ultimate tensile strength used in our damage model.

A non-linear curve fitting tool (lsqnonlin) in MATLAB (version R2015a) was used to determine the above-mentioned mechanical damage parameters in equations (5) and (6) by fitting equations (5) and (6) to these experimental data. The plot obtained for this curve fitting process has been shown in Figure 2. The coefficient of determination for this model fit was determined to be 0.7638, and the original variance of the experimental data about the mean was determined to be 1.134 MPa. This signifies that the model fit accounts for 76.38% (0.866 MPa) of the original experimental variance. From Figure 2, it can be seen that there are a few experimental data points that deviate from the non-linear curve fit. However, these outliers could not be avoided as this is the best possible fit we could obtain through the curve fitting process outlined above.

Figure 2.

Fitting of the stress–life relationship in our model to the stress-life data obtained from the disc fatigue failure experiments of Green et al.¹⁵

It must be noted that the experimental stress-life data used for calibrating equations (5) and (6) were only for fully reversed loading, with a mean stress of zero. Therefore, the value of the mean stress parameter b in equation (6) could not be determined from this curve fitting process. The procedure followed to determine b has been outlined in the “Calibration of model parameters r and b” section, in which model calibration is discussed.

Damage due to aging

Strehler⁴⁴ proposed that the damage evolution in a mature biological structure due to aging phenomena can be considered a linear function of time. The damage model presented in this article adopts the following linear equation for damage evolution due to aging phenomena

D_{a} = C_{a} \times age

(16)

where $D_{a}$ is the damage due to aging, age is the age of the individual being evaluated, and $C_{a}$ is a proportionality constant.

To determine the aging proportionality constant $C_{a}$ in equation (16), disc degeneration data for the general population recorded by Thompson et al.³⁵ were utilized. These data were presented as average grades of degeneration in the intervertebral discs from spines of donors with ages ranging from 15 to 90 years. The average grades ranging from Grade I to Grade V in the grading system given by Thompson et al.³⁵ were scaled using linear interpolation, corresponding to the values of the total damage variable $D_{tot}$ ranging from 0 to 1.

According to equation (3), the total damage variable obtained in this manner represents the total damage to the intervertebral disc due to mechanical loading $(D_{m})$ and aging phenomena $(D_{a})$ . Since the data provided by Thompson et al.³⁵ were for a general population, $D_{m}$ was estimated by applying the loading case of walking without head-supported mass to our calibrated damage model (developed in the “Model calibration” section). This mechanical damage $D_{m}$ was subtracted from the total damage $D_{tot}$ obtained from the disc degeneration data, in order to obtain an estimate for the damage due to aging. A least squares method was used to fit a straight line to this estimate while constraining the line to pass through the origin, in order to obtain the aging proportionality constant $C_{a}$ . Figure 3 shows a plot of this straight line fit with a 95% confidence interval.

Figure 3.

To determine the aging constant $C_{a}$ , a straight line was fit to the estimates of damage due to aging obtained from the disc degeneration data provided by Thompson et al.³⁵

Model calibration

In the “Model development” section, we determined the parameters $a$ , $β,$ $M_{0}$ , and $C_{a}$ . Further calibration of the model is required to determine the mean stress parameter b and the recovery parameter r. The following section elucidates the procedure followed to determine these parameters.

Cyclic loading input for the damage model

The full-body musculoskeletal model shown in Figure 4 was used to simulate the forces in all three directions (inferior–superior, anterior–posterior, and lateral) on the cervical spine joints. The following three loading cases were simulated for the purpose of providing the appropriate loading input to the damage model presented in this article: neck loading during one gait cycle due to normal walking, neck loading during one gait cycle due to walking while wearing an Advanced Combat Helmet (ACH) weighing 1.43 kg, and neck loading during a single impact period of a fast boat ride. The musculoskeletal neck model used for this simulation consisted of eight cervical vertebral joints, 24 degrees of freedom of the head and neck, and 84 muscle fascicles.^45,46 The present study used an improved version⁴⁷ of the neck model developed by Vasavada et al.,⁴⁶ which was also validated in other studies.^48,49 This neck model was further modified and integrated with a full-body model⁴⁵ shown in Figure 4. Maximum neck muscle force was revised according to Chancey et al.⁵⁰ and Nelson et al.⁵¹ The segment mass inertia properties of the cervical spine and head were adopted from De Jager⁵² and Van der Horst.⁵³ A medium-sized ACH was assembled onto the head and neck model using the mass and inertia properties listed in Table 1. In order to mimic the neck loading due to riding a fast boat, measured lumbar accelerations of a dummy on the boat were applied to the pelvis and torso of the same musculoskeletal model, without the ACH.

Figure 4.

The detailed neck and full-body musculoskeletal model (left) with predicted muscle contraction activation level (right) at a point of the walking gait cycle (middle). The predicted muscle contraction forces contribute to compressive forces in the intervertebral discs.

Table 1.

Mass and inertia properties of a medium-sized ACH.

	Mass (kg)	COM (cm)		MI (kg cm²)
	Mass (kg)	x	y	I_x	I_y	I_z
Head	4.69	1.0	3.5	181.0	173.0	236.0
ACH	1.43	–0.7	8.9	130.0	70.0	130.0
Head + ACH	6.12	0.4	4.8	357.2	272.0	421.2

ACH: Advanced Combat Helmet: COM: center of mass; MI: moment of inertia; HACS: head anatomical coordinate system.

ACH data were obtained from LaFiandra et al.⁵⁴ The head-supported mass COM position is presented relative to the HACS associated with Frankfurt planes. The mass is distributed symmetrically with respect to the lateral z-direction and the MI is expressed relative to the corresponding COM. Here, x is the anterior–posterior direction, y is the inferior–superior direction, and z is the lateral direction.

These simulations were performed using in-house developed computational software, CoBi-Dyn.⁵⁵ This software simulates full-body dynamics using the reduced (or internal) coordinate approach and the multibody dynamics formulation based on spline joint.⁵⁶ The use of CoBi-Dyn for simulations is well accepted in the biomechanics literature.^45,57–59 The simulations were aimed at estimating the intervertebral joint forces during normal walking, walking with the medium-sized ACH, and while riding a fast boat.

From these simulations, the neck loading in all three directions (x, y, and z) at the C2C3 intervertebral joint was extracted for the following cases: one gait cycle of normal walking, one gait cycle of walking while wearing an ACH weighing 1.43 kg, and a single impact period of a fast boat ride. From these extracted load histories, the histories of the effective stress measure known as the signed Von-Mises stress were computed according to equation (13). The stresses $σ_{yy}$ , $σ_{xy}$ , and $σ_{yz}$ were computed from the extracted loads by considering the cross-sectional area of the C2C3 intervertebral disc to be 308 mm².¹⁸ This disc area of 308 mm² is in agreement with the average disc area for the C2C3 disc as measured by Moroney et al.⁶⁰ As a first approximation to the damage model presented in this article, a uniform stress distribution was assumed in the disc material.

Figure 5 depicts the computed signed Von-Mises stress histories at the C2C3 intervertebral joint, during one gait cycle of walking with and without the ACH. This walking gait cycle starts at a heel strike (0% at 0.583 s) and ends when this heel strikes the ground again (100% at 1.817 s). It was concluded that the intervertebral joint compressive loading was increased significantly while walking with head-supported mass as compared to normal walking.

Figure 5.

Signed Von-Mises stress history at the C2C3 joint of the cervical spine during one gait cycle of walking with and without head-supported mass.

Figure 6 depicts the computed signed Von-Mises stress history at the C2C3 intervertebral joint, during one short cycle of a fast boat ride. This cycle includes a single impact period and the subsequent period dominated by hydrodynamic lift and buoyancy forces.⁶¹ It must be noted that the stress histories depicted in Figures 5 and 6 represent the cumulative result of the inertial effect (due to motion) and muscle contraction induced joint compression.

Figure 6.

Signed Von-Mises stress history at the C2C3 joint of the cervical spine during a single impact period on a fast boat.

It must be noted that the C2C3 disc was chosen as the subject of our study because it is the intervertebral disc in closest proximity to the point of load application (head-supported mass/helmet). This would lead to a higher magnitude of stresses in the C2C3 disc as compared to lower cervical spine discs. Therefore, the damage predicted by the model would be highest at the C2C3 disc, giving a more conservative damage estimate. Moreover, Bogduk and Marsland⁶² observed a higher incidence of neck pain arising from damage to the C2C3 FSU, as compared to lower levels of the cervical spine. For these reasons, the C2C3 disc was chosen as the subject of our study.

Rainflow cycle counting

In order to implement the equations for variable amplitude loading developed in the “Equation for mechanical damage evolution” section, the rainflow cycle counting method⁶³ was used to extract stress cycles from the complex loading histories shown in Figures 5 and 6. Further elaboration of the algorithm used by the rainflow cycle counting method can be found in Irvine⁶⁴ and Downing and Socie,⁶⁵ and some examples of its usage in fatigue analysis can be found in previous works.^66–68 This method can be used to obtain the mean stresses and stress amplitudes of each extracted cycle from the complex loading histories shown in Figures 5 and 6. These mean stresses and stress amplitudes were then fed as input to the damage model presented in this article.

In order to estimate the damage evolution in the disc on a timescale of years, we need to estimate the average number of times per day that the loading histories shown in Figures 5 and 6 are applied to the damage model. For the case of normal walking without head-supported mass, one day was represented by 9000 steps on the basis of results from a meta-analysis by Bohannon.⁶⁹ For the case of walking with head-supported mass, a soldier was assumed to walk 12 miles a day on average, based on the typical length of a loaded march required to gain the Expert Infantryman Badge. Assuming there are 1900 steps in 1 mile,⁷⁰ 12 miles per day translates to 22,800 steps per day. It must be noted that the loading over one gait cycle shown in Figure 5 signifies two complete steps. For the case of loading due to riding a fast boat, it was assumed that a typical training session in a fast boat is 15 min long.⁷¹ Therefore, in this case, the 2.85-s long loading history shown in Figure 6 was applied 315 times per day to the damage model.

Calibration of model parameters r and b

As discussed in the “Curve fitting to determine mechanical damage parameters” section, the experimental stress-life data¹⁵ used for fitting equations (5) and (6) were only for fully reversed loading, with a mean stress of zero. Therefore, the value of the mean stress parameter b in equation (6) could not be determined from this curve fitting process. Moreover, due to a lack of quantitative data in the literature regarding the amount of mechanical damage recovery in the disc, it was not possible to obtain an estimate for the recovery parameter r. The procedure followed to estimate these model parameters r and b has been outlined in the following paragraphs.

The damage model was run for the case of normal walking without head-supported mass. Since the values of r and b have not been determined yet, the model was run for a range of values of r and b. The estimates of damage evolution predicted by the model were compared against disc degeneration data for cervical and lumbar intervertebral discs from the generic population, from multiple sources in the literature.^36–38 It must be noted that both cervical and lumbar disc degeneration data were used for this comparison, due to the scarcity of data in the literature exclusive to cervical discs. These data represent the number of discs at various stages of degeneration in different age groups. Each stage of degeneration is represented by a grade which was scaled corresponding to the range of values of the total damage variable $D_{tot}$ , as outlined in the “Damage due to aging” section. By comparing these scaled data to the model results for a range of values of r and b, the values of r and b for which the model results best agreed with the scaled data were determined. A parametric study investigating the effect of the parameters r and b on the model results has been presented in the “Results and discussion” section.

Now, having determined all the parameters required by the damage model, we can apply the model to predict damage evolution in intervertebral discs over long periods of time due to walking with head-supported mass and due to riding a fast boat. This model application process has been discussed in the “Results and discussion” section.

Results and discussion

Model parameters

The value of the ultimate tensile strength $σ_{u} = 6.25 MPa$ for the disc annulus material was used in the damage model presented in this article, as outlined in the “Curve fitting to determine mechanical damage parameters” section. Table 2 contains the values of the model parameters $a$ , $β$ , and $M_{0}$ determined from the non-linear curve fitting process presented in the “Curve fitting to determine mechanical damage parameters” section and the value of the aging proportionality constant $C_{a}$ which was estimated as explained in the “Damage due to aging” section.

Table 2.

Damage model parameters determined from curve fitting processes.

Parameter	a	β	M₀ (MPa)	C_a (year⁻¹)
Value	9.8865	2.9737	17.8154	0.0030

As explained in the “Calibration of model parameters r and b” section, the recovery parameter r and the mean stress parameter b could not be determined from curve fitting processes. The “Calibration of model parameters r and b” section also outlines the procedure followed to estimate the model parameters r and b. After a parametric study, it was seen that a representative value of r = 0.5 and a value of b = 1.5 produced a good fit between the damage model prediction and the disc degeneration data.^36–38Figure 7 shows a plot comparing the disc degeneration data and the damage evolution predicted by the model over a period of 80 years, using the calibrated values of r = 0.5 and b = 1.5, for the loading case of normal walking without head-supported mass. The coefficient of determination for this fit between the damage model results and the disc degeneration data was determined to be 0.7729, and the original variance of the experimental data about the mean was determined to be 0.06. This signifies that the model fit accounts for 77.29% (0.046) of the original experimental variance.

Figure 7.

Comparison of damage model prediction with disc degeneration data,^36–38 for the loading case of normal walking without head-supported mass. The calibrated values of r = 0.5 and b = 1.5 were used.

Damage due to walking with head-supported mass

The fully calibrated model can now be used to predict fatigue damage evolution in the intervertebral discs of military personnel due to walking with heavy head-supported mass. The loading history over one gait cycle of walking with a medium-sized ACH as shown in Figure 5 was applied to the damage model. As explained in the “Rainflow cycle counting” section, a soldier was assumed to walk 22,800 steps (12 miles) per day, where 2 complete steps represent one gait cycle. Figure 8 shows the damage evolution predicted by the model for this case, over a period of 30 years.

Figure 8.

Predicted damage evolution in cervical intervertebral discs of military personnel due to walking with a medium-sized ACH, over a period of 30 years.

From Figure 8, it can be seen that the total damage predicted by the model due to walking with the ACH reaches a value of 1 (complete damage) in a period of 26 years. This is in stark contrast to the predicted damage evolution for the loading case of walking without head-supported mass. The damage evolution in discs of the generic population due to normal walking without head-supported mass is depicted in Figure 7, which shows that the total damage predicted by the model only reaches a value of 0.13 at 26 years. Although the lifestyle of a typical soldier is highly simplified in this damage model, these results indicate that heavy head-supported equipment worn by military personnel may accelerate the rate of intervertebral disc degeneration.

Damage due to riding a fast boat

The damage model was also used to predict fatigue damage evolution in the intervertebral discs of military personnel due to riding a fast boat. The loading history over the 2.85-s long single impact period on a fast boat as shown in Figure 6 was applied to the damage model. As explained in the “Rainflow cycle counting” section, it was assumed that the typical training duration in a fast boat is 15 min per day.⁷¹Figure 9 shows the damage evolution predicted by the model for this case, over a period of 13 years.

Figure 9.

Predicted damage evolution in cervical intervertebral discs of military personnel due to riding a fast boat, over a period of 13 years.

From Figure 9, it can be seen that the total damage predicted by the model due to riding a fast boat reaches a value of 1 (complete damage) in a period of 12.5 years. This is in stark contrast to the predicted damage evolution for both the loading cases of walking without and with head-supported mass. Clearly, the damage evolution predicted by the model for both these cases (represented by Figures 7 and 8, respectively) at 12.5 years is significantly lower than the complete damage $(D_{tot} = 1)$ predicted at 12.5 years for the loading case of riding a fast boat. These results indicate that the high-frequency cyclic loads on the cervical spine experienced by military personnel while riding a fast boat may accelerate cervical disc degeneration at a rate significantly higher than that due to walking with heavy head-supported mass. It must be noted that the x-axes of Figures 8 and 9 represent the number of years elapsed after a person has joined the military. It may also be noted that in both Figures 8 and 9, the initial value of the damage due to aging is a non-zero value. This is based on the assumption that a person joining the military would typically be around 20 years of age, resulting in an initial value of 0.06 (C_a × 20) for the damage due to aging.

Effect of r and b on the predicted damage evolution

As explained in the “Calibration of model parameters r and b” section, the recovery parameter r and the mean stress parameter b could not be determined from curve fitting processes. They were estimated by the model calibration process presented in the “Calibration of model parameters r and b” section. Therefore, it is important to investigate the effect of these parameters on the damage predicted by the model. Since both the parameters r and b only affect the evolution of mechanical damage, this parametric study only investigates the effect of r and b on the mechanical damage $D_{m}$ .

First, a parametric study of the recovery parameter r was performed in which the results of the damage model were compared for various values of r ranging from 0 to 0.9. The loading case of walking with the medium-sized ACH was applied to the damage model to perform this study. Figure 10 shows mechanical damage evolution curves over a period of 40 years for various values of r. From Figure 10, it can be observed that varying the recovery parameter has a considerable effect on the results. For a change of r from 0.5 to 0.7 (40% increase), the value of the predicted mechanical damage at 10 years changed from 0.23 to 0.11 (52% decrease).

Figure 10.

Comparison of the mechanical damage evolution predicted by the model for different values of r.

Then, a parametric study of the mean stress parameter b was performed in which the results of the damage model were compared for various values of b ranging from 0 to 2. The loading case of walking with the medium-sized ACH was applied to the damage model to perform this study. Figure 11 shows mechanical damage evolution curves over a period of 40 years for various values of b. From Figure 11, it can be observed that varying the mean stress parameter has a considerable effect on the results. For a change of b from 1.5% to 2 (33% increase), the value of the predicted mechanical damage at 10 years changed from 0.23% to 0.05 (78% decrease).

Figure 11.

Comparison of the mechanical damage evolution predicted by the model for different values of b.

As outlined in “Calibration of model parameters r and b” and “Model Parameters” sections, calibrated values of r = 0.5 and b = 1.5 were adopted through the process of model calibration. Further experiments are needed to quantify the recovery parameter and the mean stress parameter in order to improve the accuracy of the damage model presented in this article.

Limitations and future work

The problem discussed in this article poses several difficult challenges, ranging from obtaining cervical spine loading data from the field of combat, to quantifying biological factors such as recovery and aging in a computational model.

One limitation of this article is that it only considers damage due to walking, whereas soldiers may even run for several miles during extensive physical training. It is very difficult to predict the physical activity pattern of soldiers during combat operations, and an attempt has been made to quantify this to some extent in this article. Moreover, we were unable to acquire relevant data pertaining to degeneration of the intervertebral discs of military personnel, which would be required to validate the damage model for the military population. In order to obtain such data, we would need to conduct a longitudinal study, which tracks soldier occupation and keeps count of steps using a tracking device. Then, we would need to conduct histopathological studies on the discs upon soldier death, in order to characterize the damage. If such military disc degeneration data could be obtained, the damage model results for the military population can be validated in order to improve the usability of the damage model.

In order to perform this validation, we would have to obtain separate intervertebral disc degeneration datasets spanning several years for the two categories of military personnel studied in this article (soldiers walking with heavy head-supported mass and soldiers riding a fast boat). Then, in a manner similar to the method adopted in the “Calibration of model parameters r and b” section, the damage model could be validated for the military population by comparing the model predictions with the appropriate dataset for these two categories of military personnel. This can be pursued in future work on the damage model.

The values of the recovery parameter r and the mean stress parameter b could not be determined from curve fitting processes, due to the reasons explained in the “Calibration of model parameters r and b” section. Calibrated values of r = 0.5 and b = 1.5 were adopted through the process of model calibration. Further experiments are needed to obtain disc fatigue failure data for a cyclic loading condition with non-zero mean stress and to quantify the recovery parameter r. This would lead to more reliable estimates of r and b, thus improving the accuracy and reliability of the damage model.

Some of the experimental data used for curve fitting and model calibration were for the lumbar spine, due to the scarcity of cervical spine data available in the literature. Furthermore, the experimental data available were for tensile testing of the disc annulus, whereas the disc is subjected to predominantly compressive stresses in the loading scenarios considered in this article. Future experiments could be aimed at generating cervical spine data for compressive stresses, which could be used to improve the accuracy and reliability of the damage model.

In this article, we have presented a damage model that predicts damage evolution over long periods of time for a given stress history. We have considered a uniform stress distribution throughout the disc annulus, as a first approximation of the damage model. Our next step is aimed at developing a finite element model of a cervical spine segment, which would enable us to better understand the non-uniform stress distribution pattern in the disc annulus. For this, the annulus fibrosus would be modeled as an anisotropic composite, resulting in a non-uniform stress distribution in the disc. The effective signed Von-Mises stress histories of each element from this finite element model can then be fed as input to the computational damage model, which would enable us to investigate the pattern of damage initiation and progression in the intervertebral disc annulus. The location of maximum stress in the annulus would be the most probable location of damage initiation.

Conclusion

This article presents a structured framework for the prediction of fatigue damage in the cervical intervertebral disc annulus, due to long-term cyclic loading experienced by military personnel during combat operations such as walking with heavy head-supported equipment and riding a fast boat. The damage model proposed in this article considers damage evolution based on different causes of damage to the disc, and this is a novel contribution to the field of intervertebral disc degeneration. The methodology used in this article quantifies mechanical damage due to cyclic loading with stress cycles of variable amplitude, based on a non-linear damage evolution model. The mechanical damage is computed based on the effective stress measure for multiaxial fatigue estimation known as the signed Von-Mises stress, which takes into account the multiaxial stress state in the disc annulus. Equations are developed to predict the damage due to aging and damage recovery during daily rest periods. The model parameters in these equations were obtained from curve fitting and model calibration, using experimental data from the literature. Then, the fully calibrated damage model was implemented to predict damage evolution in the discs of military personnel, for the loading scenarios of walking with heavy head-supported equipment and riding a fast boat. The model predictions indicate that the complex cyclic loads on the cervical spine typically experienced by military personnel during combat operations may significantly accelerate the rate of intervertebral disc degeneration. However, it must be kept in mind that these damage model predictions for the military population could not be validated, due to a lack of relevant intervertebral disc degeneration data in the literature. If such military data could be obtained, the damage model can be validated for the military population, which would increase the reliability of the damage model.

Long-term damage prediction in military personnel is a highly complicated problem with several challenges in obtaining experimental data for model calibration and model validation. Therefore, this article also identifies several gaps in the relevant experimental work and suggests various directions for future work on this model.

Footnotes

Handling Editor: Liyuan Sheng

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was partially supported by the US Army under grant W81XWH-14-C-003. This work was also supported in part through instrumentation funded by the National Science Foundation through grant OCI0821527.

ORCID iD

Adhitya Subramani

References

Van den Oord

Veerle De Loose

. Neck pain in military helicopter pilots: prevalence and associated factors. Mil Med 2010; 175: 55–60.

Ensign

Hodgdonl

Prusaczyk

et al . A survey of self-reported injuries among special boat operators. San Diego, CA: Naval Health Research Center, 2000, http://www.dtic.mil/docs/citations/ADA421234

Cohen

Kapoor

Nguyen

et al . Neck pain during combat operations: an epidemiological study analyzing clinical and prognostic factors. Spine 2010; 35: 758–763.

Mathys

Ferguson

. Simulation of the effects of different pilot helmets on neck loading during air combat. J Biomech 2012; 45: 2362–2367.

Manoogian

Kennedy

Duma

. A literature review of musculoskeletal injuries to the human neck and the effects of head-supported mass worn by soldier. Fort Rucker, AL: Army Aeromedical Research Lab, 2005, http://www.dtic.mil/docs/citations/ADA441139

Podichetty

. The aging spine: the role of inflammatory mediators in intervertebral disc degeneration. Cell Mol Biol 2007; 53: 4–18.

Buckwalter

. Aging and degeneration of the human intervertebral disc. Spine 1995; 20: 1307–1314.

Peng

Hao

Hou

et al . Possible pathogenesis of painful intervertebral disc degeneration. Spine 2006; 31: 560–566.

Adams

Roughley

. What is intervertebral disc degeneration, and what causes it?Spine 2006; 31: 2151–2161.

10.

Iatridis

MacLean

Ryan

. Mechanical damage to the intervertebral disc annulus fibrosus subjected to tensile loading. J Biomech 2005; 38: 557–565.

11.

Hansson

Keller

Spengler

. Mechanical behavior of the human lumbar spine. II. Fatigue strength during dynamic compressive loading. J Orthop Res 1987; 5: 479–487.

12.

C-Y

Tsai

K-H

W-P

et al . Geometric and morphological changes of the intervertebral disc under fatigue testing. Clin Biomech 2003; 18: S3–S9.

13.

Yingling

Callaghan

McGill

. The porcine cervical spine as a model of the human lumbar spine: an anatomical, geometric, and functional comparison. Clin Spine Surg 1999; 12: 415–423.

14.

Snyder

. Development of ultrasound to measure in-vivo dynamic cervical spine intervertebral disc mechanics. Boston, MA: Beth Israel Deaconess Medical Center, 2015, http://www.dtic.mil/docs/citations/ADA613867

15.

Green

Adams

Dolan

. Tensile properties of the annulus fibrosus. II. Ultimate tensile strength and fatigue life. Eur Spine J 1993; 2: 209–214.

16.

Kachanov

. Rupture time under creep conditions. Int J Fract 1999; 97: 11–18.

17.

Qasim

Natarajan

et al . Initiation and progression of mechanical damage in the intervertebral disc under cyclic loading using continuum damage mechanics methodology: a finite element study. J Biomech 2012; 45: 1934–1940.

18.

Makwana

Krishna

Yuan

et al . Towards a micromechanical model of intervertebral disc degeneration under cyclic loading. In: International mechanical engineering congress and exposition, Montreal, QC, Canada, 14–20 November 2014, pp.V003T03A012-1–V003T03A012-7. New York: ASME.

19.

Miner

. Cumulative damage in fatigue. J Appl Mech 1945; 12: A159–A164.

20.

Manson

Freche

Ensign

. Application of a double linear damage rule to cumulative fatigue. In: Fatigue crack propagation. ASTM International, 1967, http://www.astm.org/DIGITAL_LIBRARY/STP/PAGES/STP47237S.htm

21.

Manson

Nachtigall

Ensign

et al . Further investigation of a relation for cumulative fatigue damage in bending. J Eng Ind 1965; 87: 25–35.

22.

Fatemi

Yang

. Cumulative fatigue damage and life prediction theories: a survey of the state of the art for homogeneous materials. Int J Fatigue 1998; 20: 9–34.

23.

Marco

Starkey

. A concept of fatigue damage. T ASME 1954; 76: 627–632.

24.

Lemaitre

Chaboche

J-L

. Mechanics of solid materials. Cambridge University Press, 1994, https://books.google.com/books?hl=en&lr=&id=YfHTNHfjqXoC&oi=fnd&pg=PR11&dq=Lemaitre+and+Chaboche+1994&ots=kvRSFlY1uy&sig=wX5WDBI5GbereYj6kFiF04Wg9SE

25.

Pidaparti

Wang

Burr

. Modeling fatigue damage evolution in bone. Biomed Mater Eng 2001; 11: 69–78.

26.

Schechtman

Bader

. In vitro fatigue of human tendons. J Biomech 1997; 30: 829–835.

27.

Zioupos

Wang

Currey

. Experimental and theoretical quantification of the development of damage in fatigue tests of bone and antler. J Biomech 1996; 29: 989–1002.

28.

Fratzl

Weinkamer

. Hierarchical structure and repair of bone: deformation, remodelling, healing. In: Self healing materials. Springer, 2007. pp.323–335, https://link.springer.com/chapter/10.1007/978-1-4020-6250-6_15

29.

Bradford

Oegema

Jr Cooper

et al . Chymopapain, chemonucleolysis, and nucleus pulposus regeneration. A biochemical and biomechanical study. Spine 1984; 9: 135–147.

30.

O’Connell

Jacobs

Sen

et al . Axial creep loading and unloaded recovery of the human intervertebral disc and the effect of degeneration. J Mech Behav Biomed Mater 2011; 4: 933–942.

31.

Johannessen

Vresilovic

Wright

et al . Intervertebral disc mechanics are restored following cyclic loading and unloaded recovery. Ann Biomed Eng 2004; 32: 70–76.

32.

Ion

Lorand

Mircea

et al . The equivalent stress concept in multiaxial fatigue. J Eng Stud Res 2011; 17: 53–62.

33.

Nash

. Fatigue of self-healing structure: a generalized theory of fatigue failure. New York: ASME, 1967, p.68.

34.

Hansen

H-J

. A pathologic-anatomical study on disc degeneration in dog, with special reference to the so-called enchondrosis intervertebralis. Acta Orthop Scand 1952; 23: 1–130.

35.

Thompson

Pearce

Schechter

et al . Preliminary evaluation of a scheme for grading the gross morphology of the human intervertebral disc. Spine 1990; 15: 411–415.

36.

Lawrence

. Disc degeneration. Its frequency and relationship to symptoms. Ann Rheum Dis 1969; 28: 121–138.

37.

Miller

JAA

Schmatz

Schultz

. Lumbar disc degeneration: correlation with age, sex, and spine level in 600 autopsy specimens. Spine 1988; 13: 173–178.

38.

Antoniou

Steffen

Nelson

et al . The human lumbar intervertebral disc: evidence for changes in the biosynthesis and denaturation of the extracellular matrix with growth, maturation, ageing, and degeneration. J Clin Invest 1996; 98: 996–1003.

39.

Chaboche

Lesne

. A non-linear continuous fatigue damage model. Fatigue Fract Eng Mater Struct 1988; 11: 1–17.

40.

Ostergren

Krempl

. A uniaxial damage accumulation law for time-varying loading including creep-fatigue interaction. ASME T J Press Vessel Technol 1979; 101: 118–124.

41.

Dattoma

Giancane

Nobile

et al . Fatigue life prediction under variable loading based on a new non-linear continuum damage mechanics model. Int J Fatigue 2006; 28: 89–95.

42.

Lotz

. Animal models of intervertebral disc degeneration: lessons learned. Spine 2004; 29: 2742–2750.

43.

Virgin

. Experimental investigations into the physical properties of the intervertebral disc. J Bone Joint Surg Br 1951; 33: 607–611.

44.

Strehler

. The biology of aging. Washington, DC: American Institute of Biological Science, 1960, pp.170–175.

45.

Zhou

Whitley

Przekwas

. A musculoskeletal fatigue model for prediction of aviator neck manoeuvring loadings. Int J Hum Factors Model Simul 2014; 4: 191–219.

46.

Vasavada

Delp

. Influence of muscle morphometry and moment arms on the moment-generating capacity of human neck muscles. Spine 1998; 23: 412–422.

47.

Zheng

. Sex differences in human neck musculoskeletal biomechanics and modeling. Washington State University, 2011, http://search.proquest.com/openview/2c2cec56accb8359a9459b155e7f2c12/1?pq-origsite=gscholar&cbl=18750&diss=y

48.

Nevins

Zheng

Vasavada

. Inter-individual variation in vertebral kinematics affects predictions of neck musculoskeletal models. J Biomech 2014; 47: 3288–3294.

49.

Suderman

. The effect of curvature on neck muscle moment arms. Washington State University, 2012, http://search.proquest.com/openview/8b30c6036afd4c3a472b30f994f2a60a/1?pq-origsite=gscholar&cbl=18750&diss=y

50.

Chancey

Nightingale

Van Ee

et al . Improved estimation of human neck tensile tolerance: reducing the range of reported tolerance using anthropometrically correct muscles and optimized physiologic initial conditions. Stapp Car Crash J 2003; 47: 135–153.

51.

Nelson

Pandy

Myers

et al . Variation of neck muscle strength along the human cervical spine. Stapp Car Crash J 2004; 48: 397–417.

52.

De Jager

. Mathematical head-neck models for acceleration impacts. PhD Thesis, Eindhoven University of Technology, University of Limburg and TNO Crash-Safety Research Center, Eindhoven, 1996.

53.

Van der Horst

. Human head neck response in frontal, lateral and rear end impact loading: modelling and validation, 2002, http://repository.tue.nl/554047

54.

LaFiandra

Harman

Cornelius

et al . The effects of the personal armor system for ground troops (PASGT) and the Advanced Combat Helmet (ACH) with and without PVS-14 night vision goggles (NVG) on neck biomechanics during dismounted soldier movements. Natick, MA: Military Performance Div, Army Research Inst of Environmental Medicine, 2007, http://www.dtic.mil/docs/citations/ADA470312

55.

CFDRC. Cobi-dyn, 2016, http://www.cfdrc.com/CoBi-SUITE-10-47.html

56.

Lee

S-H

Terzopoulos

. Spline joints for multibody dynamics. ACM T Graphic 2008; 27: 22.

57.

Zhou

Przekwas

. A fast and robust whole-body control algorithm for running. Int J Hum Factors Model Simul 2011; 2: 127–148.

58.

Dura-Bernal

Zhou

Neymotin

et al . Cortical spiking network interfaced with virtual musculoskeletal arm and robotic arm. Front Neurorobot 2015; 9: 13.

59.

Whitley

Roos

Zhou

. Comparison of gender specific and anthropometrically scaled musculoskeletal model predictions using the Sorensen test. In: Cassenti

(ed.) Advances in human factors in simulation and modeling. Cham: Springer, 2017, pp.469–477.

60.

Moroney

Schultz

Miller

et al . Load-displacement properties of lower cervical spine motion segments. J Biomech 1988; 21: 769–779.

61.

Riley

Coats

Murphy

. Acceleration response mode decomposition for quantifying wave impact load in high-speed planing craft. Norfolk, VA: Carderock Div, Naval Surface Warfare Center, 2014, http://www.dtic.mil/docs/citations/ADA621230

62.

Bogduk

Marsland

. The cervical zygapophysial joints as a source of neck pain. Spine 1988; 13: 610–617.

63.

Matsuishi

Endo

. Fatigue of metals subjected to varying stress. Jpn Soc Mech Eng 1968; 68: 37–40.

64.

Irvine

. Rainflow cycle counting in fatigue analysis. Revis Vib 2010, http://www.vibrationdata.com/tutorials2/rainflow_counting.pdf

65.

Downing

Socie

. Simple rainflow counting algorithms. Int J Fatigue 1982; 4: 31–40.

66.

Frendahl

Rychlik

. Rainflow analysis: Markov method. Int J Fatigue 1993; 15: 265–272.

67.

Amzallag

Gerey

Robert

et al . Standardization of the rainflow counting method for fatigue analysis. Int J Fatigue 1994; 16: 287–293.

68.

Bishop

Sherratt

. Finite element based fatigue calculations. NAFEMS, 2000. The International Association for the Engineering Analysis Community, Farnham.

69.

Bohannon

. Number of pedometer-assessed steps taken per day by adults: a descriptive meta-analysis. Phys Ther 2007; 87: 1642–1650.

70.

Welk

Differding

Thompson

et al . The utility of the Digi-walker step counter to assess daily physical activity patterns. Med Sci Sports Exerc 2000; 32: S481–S488.

71.

Kearns

. Analysis and mitigation of mechanical shock effects on high speed planing boats. Massachusetts Institute of Technology, 2001, https://dspace.mit.edu/bitstream/handle/1721.1/8235/50200001-MIT.pdf?sequence=2