Shock attenuation in the human lumbar spine during walking and running

Human lumbar lordosis (LL) is a postural adaptation that facilitates normal spinal function and helps maintain balanced trunk orientation while minimizing the mechanical and metabolic demands of bipedalism (Gracovetsky and Iacono, 1987; Abitbol, 1988; Farfan, 1995; Lovejoy, 2005; Whitcome et al., 2007; Saha et al., 2007, 2008; Lovejoy and McCollum, 2010; Castillo et al., 2017). However, skeletal evidence reveals considerable variation in LL within modern human populations and among fossil hominins (Been et al., 2012, 2017a,b), suggesting an adaptive function for variations in lumbar curvature. It has been hypothesized that straighter lumbar spines (i.e. low degrees of LL) offer greater stability whereas more curved lumbar postures allow the lower spine to act like a ‘shock absorber’ during locomotion (Kapandji, 1974; Adams and Hutton, 1985; Rak, 1993; Kobayashi et al., 2008; Been et al., 2012, 2017a; Gómez-Olivencia et al., 2017). Some authors have also speculated that lordosis may allow for higher shock attenuation because spinal curves help dissipate energy associated with impacts as bending and rotational deformations rather than axial compression, increasing the effects of hysteresis within back tissues (Adams et al., 2006). Previous work has shown that dynamic motions of the lumbar spine in flexion and extension oscillate approximately 1–2 deg in amplitude per lumbar vertebral level during walking (Syczewska et al., 1999), perhaps indicating spring-like behavior that loads and unloads the spine during locomotion. However, no studies have experimentally tested how variations in posture affect the spine's ability to cope with impact-related accelerations in vivo. Here, we examined the relationship between postural variations in human LL and shock attenuation (SA) in the lumbar spine during walking and running.

During locomotion, the rapid deceleration of the body at foot contact results in an impact-related shock wave, which travels up from the ground through the body until it reaches the head (Voloshin et al., 1981, 1998; Voloshin and Wosk, 1982; Shorten and Winslow, 1992; Hamill et al., 1995; Derrick et al., 1998; Mercer et al., 2002; Edwards et al., 2012; James et al., 2014; Gruber et al., 2014; Giandolini et al., 2016). Ground reaction force (GRF) impact peaks can be substantial, with magnitudes between 0.6 and 1.0 times body weight (BW) during walking and 1.0–3.0 times BW during running (Nigg et al., 1995; Whittle, 1999). Impact loading rates can be as high as 500 BW s⁻¹ when running barefoot using a rearfoot strike pattern (Lieberman et al., 2010). Attenuation of these large, rapid impact peaks is crucial because the resulting impulse can disrupt the vestibulo-ocular reflex and potentially lead to injury and other pathology (Pozzo et al., 1990, 1991; Whittle, 1999; Daoud et al., 2012; Davis et al., 2016).

In addition to shock peak magnitude, the other component of lumbar shock transmission to consider is frequency range. Shock accelerations during locomotion are composed of three main frequency components that propagate differently through the body (Shorten and Winslow, 1992): low, mid and high. Whereas low frequencies below 10 Hz typically result from active motions rather than impact, and high frequencies around 30 Hz primarily reflect the resonant frequencies of the inertial device and the vibrations of underlying soft tissues, we focus here on the mid-frequency component between 10 and 30 Hz, which represents the frequency range related to impacts (Shorten and Winslow, 1992; Angeloni et al., 1994; Hamill et al., 1995; James et al., 2014). In particular, the mid-frequency component typically occurs between 10 and 20 Hz during barefoot and shod running (Bobbert et al., 1991; Shorten and Winslow, 1992; Derrick et al., 1998; Mercer et al., 2002; Hamill et al., 1995; Edwards et al., 2012; Gruber et al., 2014; Giandolini et al., 2016) but may be slightly higher at approximately 18–22 Hz for barefoot and shod walking (James et al., 2014).

Impact-related shocks can be attenuated passively via soft tissues, footwear and ground substrate compliance, or actively via muscles that do negative work or modify gait kinematics (Paul et al., 1978; McMahon et al., 1987; Derrick et al., 1998; Whittle, 1999; Edwards et al., 2012; Butler et al., 2003; Addison and Lieberman, 2015). For example, several studies have explored how footwear and foot strike patterns affect impact shock transmission through the musculoskeletal system (Ogon et al., 2001; Divert et al., 2005; Lieberman et al., 2010; Kulmala et al., 2013; Boyer et al., 2014; Gruber et al., 2014; Giandolini et al., 2016). However, few studies have examined impact-related shocks in the spinal column. Existing studies of spinal accelerations primarily have relied on measurements from accelerometers mounted below and above the spine (e.g. lower limb and head), and these studies have mostly focused on time-domain rather than frequency-domain accelerations (e.g. Voloshin and Wosk, 1982; Ogon et al., 2001; Delgado et al., 2013).

The purpose of this study was to investigate how variations in LL affect impact-related shock attenuation in the lumbar spine. We tested the main hypothesis that individuals with higher LL show greater spinal SA during walking and running compared with individuals with straighter lumbar postures, as speculated by previous authors (Kapandji, 1974; Adams and Hutton, 1985; Rak, 1993; Adams et al., 2006; Kobayashi et al., 2008; Been et al., 2012, 2017a; Gómez-Olivencia et al., 2017). As LL is not static and given that lumbar motions potentially affect spinal SA (Syczewska et al., 1999; Adams et al., 2006), this study also examined the relationship between dynamic changes in LL and SA. Finally, we investigated the role of intervertebral discs in spinal SA. Discs are often considered the primary passive ‘shock absorbers’ of the spine (Voloshin et al., 1981, 1998; Voloshin and Wosk, 1982; Alexander, 1997; Adams et al., 2006), yet sagittal disc shape accounts for much of the variation in lumbar curvature among adults (Shefi et al., 2013). Thus, when accounting for covariation between intervertebral discs and lumbar curvature, we predict that LL remains strongly associated with spinal SA.

Study participants were recruited from the greater Boston area for a series of three experiments conducted on different days over the course of 4 months (see Castillo et al., 2017). Only young adults 18–35 years old were recruited to minimize the potential for age-related degenerative changes in spinal posture (Schwab et al., 2006). A general health questionnaire was administered prior to experiments, and participants were excluded if they reported a history of back pain, sciatica, scoliosis, major illness or injury in the last 3 months that could compromise gait. Twenty-seven participants (14 male, 13 female) completed the experiment reported here. Participants sampled a range of heights and body masses (Table 1). The study was approved by the Committee on the Use of Human Subjects at Harvard University. Written informed consent was given prior to participation. All experimental protocols took place between 15:00 h and 18:00 h to control for the circadian effects of viscoelastic creep in spinal tissues (Strickland and Shearin, 1972; Whitehouse et al., 1974; Lampl, 1992; Botsford et al., 1994; Voss and Bailey, 1997; Tillmann and Clayton, 2001).

Before participation in experiments, individuals were scanned using magnetic resonance imaging (MRI) conducted at the Center for Brain Science Neuroimaging Facility at Harvard University using a Siemens TIM Trio (3-T) scanner. Imaging protocols are described in Castillo et al. (2017). Briefly, a standard spine-array and large-flex coils scanned participants in a neutral, supine position (legs extended and arms resting at their sides). A midsagittal, single-slice ‘localizer’ scan (repetition time TR=8.6 ms, echo time TE=4 ms; 7 mm thickness, 1.7 mm pixel⁻¹) was used to orient individuals in the scanner. From the midsagittal scan, the maximum disc height (MDH) was measured at the dorsoventral center of the disc as the greatest cranio-caudal distance along the axis of the lumbar spine for six intervertebral levels (T12/L1 through L5/S1) using ImageJ (National Institutes of Health, Bethesda, MD, USA). Assuming isometric scaling of sagittal disc thickness with body mass, average disc height was standardized for size by calculating relative maximum disc height (relMDH) as MDH cubed divided by body mass.

Two tri-axial piezoelectric accelerometers (Endevco model 35A, San Juan Capistrano, CA, USA) were secured to tiny rectangular pieces of aluminium (21.5×14.0×0.5 mm) using cyanoacrylate. The total mass of each sensor including the aluminium mount and lead wires was 2.7 g. Accelerometers were affixed firmly to the skin overlying the T12/L1 and L5/S1 vertebral levels using adhesive tape. These anatomical locations were determined by manually palpating and counting the underlying vertebral spinous processes to locate the intervertebral levels of interest, assuming the five pre-sacral vertebrae to be lumbar. The T12/L1 and L5/S1 accelerometer locations were chosen to isolate SA within the lumbar vertebral column, measuring incoming shocks from the ground at the lumbosacral joint (i.e. ‘low-back’) and outgoing impact shocks at the thoracolumbar joint (i.e. ‘mid-back’). The vertical (z) axes of accelerometers were aligned with the spine's craniocaudal axis, and the transverse (y) axes were aligned with the body's dorsoventral axis. Sensors were powered and amplified by an Isotron signal conditioner (Endevco model 2793) that passed signal sampled at 1000 Hz to a common analog-to-digital converter board, which synchronized with kinematic data (see below). To reduce soft tissue oscillations, skin laterally adjacent to each accelerometer was ‘pre-loaded’ by manually stretching and taping the skin perpendicular to the spine's craniocaudal axis using kinesiology (KT™) tape. We chose to pre-load the skin with tape and use low-mass accelerometers to increase the stiffness of the attachment between the sensor and the skin. Using low-mass accelerometers and pre-loading the skin have been shown to be effective methods for reducing motion artifact due to vibrations of the inertial sensor and skin (Saha and Lakes, 1977; Ziegert and Lewis, 1979; Nokes et al., 1984; Trujillo and Busby, 1990; Ogon et al., 2001; Forner-Cordero et al., 2008).

For visualization of whole-body movements, markers were affixed to the left and right calcaneal tuberosities, first and fifth metatarsal heads, medial and lateral maleoli, medial and lateral femoral epicondyles, greater trochanters, anterior and posterior superior iliac spines, iliac crests, acromion processes, medial and lateral humeral epicondyles, ulnar and radial styloid processes, the sternal notch, the spinous process of C7, and the frontal eminences of the forehead. Small reflective markers were also affixed to the seven spinous processes approximating the T12 to S1 vertebral levels. 3D kinematic data were captured at 200 Hz using an 8-camera infrared motion-capture system (Oqus 1 Series, Gothenburg, Sweden) and Qualysis tracking software (v.2.10).

The degree of LL was measured as the sagittal central angle (θ) using external spinal markers (Fig. 1). The size-standardized central angle method has been shown to measure lumbar curvature reliably, correlating strongly with standard radiological techniques for quantifying LL including Cobb angle (Castillo et al., 2017). The central angle, which also allows for dynamic measures of LL during in vivo experimental studies, is calculated as the ratio of the curved arc length of the lumbar spine divided by its radius of curvature. Mean standing LL (Standing LL_mean) was computed by finding the central angle of the best-fit, least-squares circle passing through the seven points representing markers on the T12–S1 spinous processes averaged over the 2 min static standing trial. To measure lordosis dynamically (Dynamic LL), the same procedure was accomplished by finding the central angle of the best-fit circle through these seven lumbar spinal markers for each kinematic frame of the trial. The mean central angle over the 30 s trial was used to quantify Dynamic LL_mean. Because dynamic changes in LL during locomotion oscillate in an approximately sinusoidal pattern, the average amplitude of lumbar lordosis angular displacement (LL_amp) was found by calculating the root mean square amplitude of Dynamic LL for an entire trial. The average velocity of lumbar lordosis angular displacement (LL_vel) was found by calculating the average first-order derivative of change in Dynamic LL with respect to time (Fig. 2).

To test whether variations in LL influence lumbar shock attenuation, ordinary least-squares (OLS) regressions were separately conducted between SA and Standing LL_mean and Dynamic LL_mean during walking and running. Focusing on shocks in the sagittal plane, resultant mean shock attenuation (MSA_r) was used as the dependent variable in these models, as previous work has suggested resultant values may be of greater importance for shock assessment than individual vertical or transverse components alone (Giandolini et al., 2016). As negative values of MSA_r indicate higher lumbar shock attenuation, the correlation between lordosis and MSA_r was predicted to be negative.

Dynamic changes in LL were analyzed to explore whether differences in MSA_r are explained by lumbar spinal motion during locomotion. Only running trials were examined because MSA_r was found to be uncorrelated with LL during walking. Assuming the lumbar spine behaves like a Euler–Bernoulli beam containing elastic and viscous elements, we assumed that the overall lumbar spine's viscoelastic response to shock vibrations would follow a Kelvin–Voigt generalization (Herrmann, 2008). Thus, the spine's elastic response to impact shocks was predicted to be proportional to LL_amp, while its damping response was predicted to be proportional to LL_vel. To account for potential covariation between curvature displacement amplitude and rate, a multiple regression tested the effects of LL_vel and LL_amp as independent variables against running MSA_r as the dependent variable. Dynamic LL_mean was also included as a model covariate to account for potential covariation between lumbar posture and spinal motion parameters.

As discs are often considered passive ‘shock absorbers’ of the spine, we also investigated the extent to which intervertebral discs explain variations in SA. An OLS regression was used to find the bivariate association between relMDH and running MSA_r. To test whether LL is a predictor of lumbar SA after controlling for covariation with disc height, the proportion of variance explained in the OLS regression was compared with a multiple regression containing running MSA_r as the dependent variable and relMDH and Standing LL_mean as independent variables.

All data processing and analyses were performed in R v3.3.2 (https://www.R-project.org/). We examined assumptions of normality in variable distributions using a Shapiro–Wilk test. Standing LL_mean and Dynamic LL_mean were found to be log-normally distributed and thus log-transformed. Paired t-tests compared differences between walking versus running gait variables. Bivariate associations between continuous variables were tested using OLS regression and Pearson's correlations, and R² values quantified the proportion of variance explained by regression models. In multiple regression models, standardized β values are reported to compare between effect sizes.

Lumbar posture was similar between resting and dynamic trials. Standing LL_mean and Dynamic LL_mean were strongly correlated with each other during walking (r=0.81, P<0.0001) and running (r=0.75, P<0.0001). The slope of the regressions for walking (slope=0.58, 95% CI=0.50–0.76) and running (slope=0.56, 95% CI=0.35–0.77) were substantially less than 1, indicating that participants' lumbar postures were straighter dynamically than when standing (Fig. 1). This effect was most pronounced among individuals with Standing LL_mean greater than 30 deg. Standing LL_mean (40.5±15.9 deg) was 46% greater than walking Dynamic LL_mean (P<0.0001) and 31% greater than running Dynamic LL_mean (P=0.0002). Comparing lumbar motion between gaits, both walking and running Dynamic LL showed a repeating pattern of oscillation with each step (Fig. 2A,B). Dynamic LL_mean was 11% greater during running than during walking (P=0.01) (Table 2). LL_amp during running was 44% greater than that during walking (P=0.0001), and LL_vel during running was 57% greater than that during walking (P=0.0005) (Table 2).

The power spectra of low- and mid-back accelerometer signals clearly separated mid-frequency impact ranges from high-frequency resonance ranges (Fig. 3). The impact frequency component ranged from 16 to 27 Hz with a peak at about 23 Hz for walking (Fig. 3A), and from 10 to 20 Hz with a peak at 16 Hz for running (Fig. 3B). The higher resonance frequency component was similar for the two gaits, ranging between 25 and 35 Hz with peaks at 32 Hz for walking and 30 Hz for running. Comparing between gaits, mean resultant peak power during running was approximately 50 times greater at the mid-back and more than 100 times greater at the low-back compared with that during walking (P<0.0001) (Table 2). Within gaits, mean resultant peak power was 46% greater at the mid- versus low-back accelerometer during walking (P=0.03), but the opposite pattern occurred during running, with peak power at the low-back being 47% greater than that at the mid-back (P=0.02). As Table 2 shows, differences in resultant peak power were primarily driven by differences in transverse rather than vertical power during walking and vertical power during running. Walking transverse peak power was twice as high as vertical power at the mid-back (P<0.0001), but there were no differences between vertical and transverse peak power at the low-back. In contrast, there were no differences in transverse peak power at the mid- and low-back during running, but running vertical peak power was 28 times larger than transverse peak power at the low-back (P<0.0001) and 16.5 times larger than transverse peak power at the mid-back (P<0.0001) (Table 2).

MSA_r was −0.77 dB during running compared with 2.10 dB during walking (P=0.003), but there were no differences between gaits for individual components of transverse or vertical shock attenuation (Table 2). There was no relationship between LL and shock attenuation during walking (Fig. 4). Walking MSA_r was uncorrelated with Standing LL_mean (r=0.25, P=0.21) and Dynamic LL_mean (r=0.15, P=0.47). However, running MSA_r correlated negatively with both Standing LL_mean (r=−0.55, P=0.004) and Dynamic LL_mean (r=−0.50, P=0.009). As negative MSA_r values indicate attenuation, OLS regressions demonstrated that a 1% increase in standing lordosis was associated with a 9.8% increase in shock attenuation, and a 1% increase in dynamic lordosis was associated with a 10% increase in attenuation. Translating these effects from decibels to signal power ratios, the OLS model predicts over a 64% reduction in low- versus mid-back signal power during running for participants with the highest Dynamic LL_mean sampled in this study. For participants with the lowest Dynamic LL_mean sampled in this study, the OLS regression predicted a 42% gain in signal power between low- and mid-back during running.

LL_amp and LL_vel had opposite associations with MSA_r during running (Fig. 5). The multiple regression model (R²=0.41) revealed that LL_amp had a negative association with running MSA_r (β=−0.40, P=0.05), suggesting that lumbar spines that underwent larger lordosis angular displacements during running had greater lumbar shock attenuation (Table 3). In addition, LL_vel had a positive association with running MSA_r (β=0.40, P=0.05), suggesting that faster rates of lordosis angular displacement were related to reduced lumbar shock attenuation. The model covariate, Dynamic LL_mean (not shown), was also a strong predictor of running MSA_r (β=−0.53, P=0.005), indicating that lumbar curvature remained strongly associated with greater shock attenuation after controlling for the effects of spinal movement.

relMDH was strongly associated with shock attenuation during running. relMDH showed a negative correlation with MSA_r (r=−0.46, P=0.02), indicating that thicker discs attenuated greater amounts of impact shock (Fig. 6A). After controlling for the effects of lumbar posture, relMDH had a negative effect on MSA_r (β=−0.38, P=0.02) (Fig. 6B), but Standing LL_mean had an even stronger negative association with MSA_r (β=−0.49, P=0.005) (Fig. 6C). As a result, standing lordosis was a stronger predictor of shock attenuation than intervertebral disc height (Table 4).

The purpose of this study was to investigate how variations in lumbar lordosis affect in vivo shock attenuation through the lower spine during locomotion in healthy humans. Lordosis was quantified in natural standing posture as well as dynamically during barefoot walking and running, and lumbar SA was measured using small accelerometers taped to the skin overlying the thoracolumbar and lumbosacral joints. Our main finding was a strong association between LL and lumbar SA during running (but not walking), which explains approximately 30% of the variation in resultant SA during running. These results suggest that for every 1% increase in LL there is a 10% increase in lumbar SA, supporting the main hypothesis of this study for running but not walking.

Resultant SA within the lumbar spine during running (−0.8±3.1 dB) was much less intense than levels of resultant SA measured previously using accelerometers attached to the tibia and sacrum (−4.0±3.1 dB) (Giandolini et al., 2016). This difference is likely explained by the fact that the majority of passive and active mechanisms for SA are associated with the lower limb (e.g. footwear, lower limb compliance, substrate stiffness), so the magnitude of shock is already mostly attenuated by the time it reaches the vertebral column (Paul et al., 1978; Whittle, 1999; Derrick et al., 1998; Butler et al., 2003; Edwards et al., 2012; James et al., 2014; Addison and Lieberman, 2015). Lower limb SA probably also explains differences in lumbar shock accelerations between walking and running. Resultant peak power measured at the low-back accelerometer during walking was two orders of magnitude smaller than that during running (Table 2), suggesting that most of the impact shock generated during walking is attenuated before reaching the spine, and only very low levels of spinal shock are left to attenuate during walking compared with running.

Examining the relationship between dynamic changes in lordosis and shock attenuation, we found support for the idea that the lumbar spine acts like a viscoelastic system to attenuate impact-related shocks during running. Lumbar spinal motion during walking and running showed a mostly double-peaked pattern over the gait cycle, and the total change in lordosis angular displacement reported here is consistent with segmental spinal displacements reported in previous studies (Syczewska et al., 1999; Levine et al., 2007; Dimitriadis et al., 2011). However, motion of the lumbar spine revealed that less dynamically compliant lumbar spines (i.e. smaller amplitudes of lordosis angular displacement) are associated with lower levels of SA during running. Additionally, damping behavior was observed from the association between SA and lower rates of lordosis angular displacement. Thus, lower back stiffness and damping properties may co-vary with LL to modulate lumbar SA (Fig. 5). However, which tissues account for this behavior and how the tissues respond to lumbar shocks remain unknown. Numerous elastic (e.g. tendons) and viscous (e.g. discs) structures in the lumbar spine can passively respond to impact shocks, but the capacity for spinal tissues to store large amounts of strain energy and dissipate heat effectively in structures like discs (which are mostly avascular) is likely minimal (Alexander, 1997; Adams et al., 2006). Thus, we speculate that trunk muscle activity is the primary factor influencing dynamic lumbar stiffness and damping behavior. Further research is needed to understand how muscle activation affects lumbar motion to influence spinal shock attenuation during human gait.

Intervertebral discs account for the majority of variation in lumbar curvature in adults (Shefi et al., 2013), but another key finding of this study is that the central relative thickness of the discs explained only 20% of the variation in resultant SA during running (Fig. 6). After controlling for the effects of LL and disc thickness, variations in lumbar curvature had a much stronger influence on SA than the discs, which together explained more than twice the variation (44%) in lumbar SA compared with disc thickness alone. Contrary to the view that intervertebral discs are the primary ‘shock absorbers’ of the lumbar spine (Voloshin et al., 1981, 1998; Voloshin and Wosk, 1982; Alexander, 1997; Adams et al., 2006), these findings suggest that sagittal spinal shape may be more important for attenuating impact accelerations than discs.

Because the hypotheses tested in this study require reliable measurements of SA in the lumbar spine, we examined accelerations within the frequency domain using PSD analysis, which revealed distinct power peaks centered at about 23 Hz for walking and 16 Hz for running (Fig. 3), both within the ranges of impact frequencies reported previously (James et al., 2014; Bobbert et al., 1991; Shorten and Winslow, 1992; Derrick et al., 1998; Mercer et al., 2002; Hamill et al., 1995; Edwards et al., 2012; Gruber et al., 2014; Giandolini et al., 2016). We also found the resonance frequency of the inertial sensor to be centered at approximately 31 Hz for both walking and running, close to the mean natural vibrating frequency of 30 Hz reported by Kitazaki and Griffin (1995), who also used low-mass accelerometers attached to the skin overlying lumbar spinous processes. Together, these results increase the confidence in the methods used to measure lumbar shock transmission in this study.

This study further highlights the importance of analyzing both transverse and vertical components of impact shock during locomotion. As previous authors have concluded, transverse and resultant accelerations provide useful information for shock assessment compared with analyses of vertical accelerations alone (Giandolini et al., 2016). For example, the dimensional contributions of peak power explained differences in SA between gaits and the relationship between SA and LL. The main difference between walking and running accelerations was the elevated transverse peak power during walking. To reiterate, walking transverse peak power was twice as high as vertical power at the mid-back, causing an average gain in signal, rather than attenuation – which is not unusual in studies of impact shock and most commonly attributed to differential body segment motion, vibrations near resonance frequency or joint kinematics during stance phase (Shorten and Winslow, 1992; Hamill et al., 1995; Gruber et al., 2014). Our finding of high transverse shock accelerations during walking also may be due to overall greater antero-posterior GRF impulses reported for barefoot walking. Nilsson and Thorstensson (1989) examined GRFs according to foot strike patterns and showed that individuals walking barefoot at 1.5 m s⁻¹ had antero-posterior braking impulses that were roughly two-thirds greater than those of individuals using a rearfoot strike pattern while running at 3 m s⁻¹. As previous studies of shock attenuation during barefoot walking have mostly considered vertical shock accelerations (e.g. James et al., 2014), it is unknown whether elevated transverse shocks found in this study are representative of barefoot walking overall.

During running, peak power in the vertical axis was 16–28 times the peak power in the transverse axis. This finding suggests that the correlation between LL and SA is primarily driven by vertical rather than transverse forces. Though the mechanisms underlying dimensional differences in running acceleration peaks and their relationship to lumbar spinal shape are unclear, one explanation may be that vertical peak power is related to gait kinematics. Greater pelvic tilt (and presumably higher LL) is known to be associated with longer strides and stiffer lower limb postures (Levine and Whittle, 1996; Vialle et al., 2005; Franz et al., 2009; Hamill et al., 2009). Thus, it is possible that increased lower extremity stiffness causes higher levels of vertical forces through the spine during running. However, future work is necessary to better understand how lumbopelvic motions interact and to what extent they influence spinal shock attenuation.

There were several notable limitations of this study. First, only healthy young adults (mostly students) were recruited, thus sampling a narrow range of variation in physical activity levels, age and other factors that may influence lumbar SA. Second, this study relied on external measures of LL. Although the central angle method used here strongly correlates with radiological standards (e.g. Cobb LA), this method is less precise than radiological imaging methods (for discussion, see Castillo et al., 2017). Third, these experiments were conducted on barefoot participants. Although participants were instructed to use a rearfoot strike pattern during all trials, we did not conduct a rigorous kinematic analysis of foot strike pattern. There is extensive literature demonstrating the effects of footwear and strike pattern on impact during locomotion, which could have influenced the results in this study (Boyer et al., 2014; Divert et al., 2005; Lieberman et al., 2010; Kulmala et al., 2013; Giandolini et al., 2016). However, a post hoc qualitative examination of force plate data collected in this study found that 25 out of 27 participants had clearly visible and consistent impact transients in their vertical GRFs. As forefoot strike patterns are associated with reduced impact transients (Lieberman et al., 2010), most of the study participants are assumed to have used rearfoot strikes. We are therefore confident that foot strike patterns did not have a strong effect on the relationship between LL and lumbar SA measured in this study. Fourth, we analyzed only one speed for each gait. Although speed may influence SA, our goal was to study only moderate gait speeds during which we could establish baseline measures of lumbar SA and spinal motion. Future studies would benefit from examining the relationship between LL and lumbar SA across speeds. A fifth limitation of this study is that it focused on analyses of intervertebral disc thickness rather than other aspects of disc morphology. We relied on maximum disc height in the middle of the disc because we assumed this metric would be associated with the mechanical response of the nucleus puloposus, which is often considered to be the source of intervertebral shock absorption. But other aspects of disc morphology, such as disc wedging, may be biomechanically relevant, and their role in spinal SA should be examined further. Finally, this study did not examine other features of active SA such as lower limb kinematics or trunk muscle activation. Because the focus of this experiment was to understand the relationship between lumbar posture (and changes in posture) on spinal shock transmission, only variables related to spinal motion and accelerations were examined. However, follow-up research should consider testing how lower limb joint stiffness and trunk muscle activation may affect the patterns of SA found here.

Overall, this study demonstrates that LL has a functional influence on shock transmission through the human axial skeleton during high-impact dynamic activities such as running, further underscoring the important role of the lumbar spine in transferring energy between the upper and lower body during gait (Gracovetsky and Iacono, 1987; Syczewska et al., 1999; Grasso et al., 2000). But lordosis is only one structural component in a highly integrated biomechanical system involving many parameters (e.g. limb kinematics, pelvic tilt, thoracic kyphosis, etc.), all with complex and interacting effects. Although this study was not able to directly test the fundamental cause of higher levels of SA among individuals with greater lordosis, we speculate that correlated kinematic patterns of LL, stride length, limb compliance and foot strike may explain differences in shock transmission. Individuals with reduced LL tend to use slower walking speeds and shortened strides, and individuals with higher LL tend to use longer strides when running (Grasso et al., 2000; Sarwahi et al., 2002; Hirose et al., 2004; Franz et al., 2009). Longer strides are also associated with greater reliance on rearfoot strike patterns and increased leg stiffness, both characteristics hypothesized to increase impact forces and injury risk (Lieberman et al., 2015; Davis et al., 2016). These factors may provide clues for understanding the link between variations in spinal posture and shock transmission in the human spine.

The relationship between SA and LL found here may also help explain some of the observed differences in lumbar curvature among modern athletes (Been and Kalichman, 2014). Whereas sprinters, long-distance runners and soccer players tend to have more curved lumbar spines, body builders and swimmers – athletes who do not experience repetitive impacts but employ high levels of isometric contractions in the upper body – tend to have much straighter lumbar spines (Uetake and Ohtsuki, 1993; Wodecki et al., 2002). Evidence from this study suggests that LL is higher among athletes for whom running and dynamic impacts are key, but athletes who do not experience high levels of impact-related loading tend to have much straighter lumbar spines, possibly for stability. Though many other physical characteristics and fitness variables co-vary with LL differences and spinal function (see Been and Kalichman, 2014), this study contributes further support to the notion that LL differences have a considerable effect on spinal function.

Given the benefits of greater lumbar curvature for SA, a final question to consider is why there are such high levels of variability among modern humans and between hominin groups. Fossil estimates of lumbar curvature suggest that most hominins had degrees of lordosis within the modern human range (Been et al., 2012). However, Neanderthals had much straighter lumbar spines at the extreme low range of modern variation, suggesting an adaptive function for lordosis variation (Gómez-Olivencia et al., 2017; Been et al., 2017a,b). Although differences in upper body size and shape may be one source of lordosis variability (Castillo et al., 2017), another hypothesis is that novel environmental conditions explain higher levels of LL variability today. As we recently argued, reduced physical activity levels and novel behaviors since the Industrial Revolution – such as sleeping on soft mattresses and prolonged sitting in chairs throughout the day – may have led to abnormally low patterns of spinal loading and weaker, less stable back tissues (Castillo and Lieberman, 2015). But testing this hypothesis requires detailed comparative studies of populations around the world including hunter–gatherers and non-industrial societies that vary in activity levels and spinal loading behaviors (e.g. habitually carrying heavy loads, using harder sleeping surfaces, etc.).

Another hypothesis is that the evolution of human lumbar spinal curvature represents tradeoffs between competing selection pressures. As other authors have shown, the relative strength of the trunk muscles supporting the spine has been shown to be an important factor underlying lumbar curvature variations (Kim et al., 2006; Elsayed et al., 2018). However, the strength of the hypaxial versus epaxial muscles shows a complex tradeoff with sagittal lumbar flexibility to drive variations in lower back curvature (Castillo et al., 2017), possibly representing an underlying adaptive constraint on lumbar spinal posture. Another important tradeoff may be between injury risk and effective shock dissipation. On the one hand, spinal attenuation of impact-related shocks may be beneficial for head stabilization, preservation of vision during locomotion, and dynamic stability (Pozzo et al., 1990, 1991; Whittle, 1999). On the other hand, high levels of LL and SA, especially via passive means, may contribute to tissue strain and injury (Hamill et al., 1995). Clinical evidence suggests a link between LL and spinal pathology (Berlemann et al., 1999; Umehara et al., 2000; Kumar et al., 2001; Rajnics et al., 2002; Labelle et al., 2005; Barrey et al., 2007; Chen and Wei, 2009). Furthermore, the results of this study show that increased SA is associated with increased dynamic lumbar compliance, possibly suggesting higher amounts of strain in back tissues. Our results demonstrate that the passive effects of intervertebral disc height account for about 20% of SA, further suggesting potentially high strain on discs during dynamic changes in lumbar posture. This may be important, for instance, among individuals with back pain who show reduced SA capacity (Voloshin and Wosk, 1982) and a reduction in disc height following prolonged dynamic activities such as running (Dimitriadis et al., 2011). Thus, curved lumbar spines may allow for greater SA at the cost of increased risk of injury while straighter lumbar spines, which may be more stable and less injury prone, are less able to contend with shock forces generated during running.

Altogether, we speculate that this may suggest the less lordotic Neanderthal lumbar spines may have been better adapted for stiffness and stability at the expense of a reduced capacity for SA during dynamic activities such as running. Lordosis was likely an early bipedal adaptation for balancing the mass of the bipedal upper body over the lower limb, but why earlier hominins such as australopiths evolved degrees of lordosis close to the modern human average remains unclear (Been et al., 2012). Although australopiths were unlikely to be running long distances to the same degree as later Homo (Bramble and Lieberman, 2004), another possible explanation is that modern human-like lordosis was useful for climbing, a behavior that may have been retained among early hominins as evidenced by their thorax and shoulder morphology (Stern and Susman, 1983). Much like modern rock climbers today, who tend to have more curved lumbar spines (Förster et al., 2009), perhaps the increased lower back curvature among australopiths was useful for climbing and the need for sufficient lordosis-related spinal compliance (Castillo et al., 2017). More research is needed to test these hypotheses to interpret LL variability in the hominin fossil record.

We thank Connie Hsu, Ross Mair, Stephanie McMains and Tammy Moran for help with MRI scanning and analysis. We also thank David Williams, Madhusudhan Venkadesan, Anna Warrener, Erik Otárola-Castillo, Brian Addison, Heather Dingwall and Andrew Yegian for helpful discussions related to study design and analyses. Special thanks to David Pilbeam, Andrew Biewener, Guoan Li and two anonymous reviewers for their comments, which have greatly improved the manuscript.