Load carrying with flexible bamboo poles: optimization of a coupled oscillator system

Backpacks with non-compliant straps are commonly used in western nations when carrying moderate to heavy loads. Elastic suspension systems – such as spring-loaded backpacks (Rome et al., 2005, 2006) or flexible poles (Castillo et al., 2014) – can reduce energetic expenditure of load carriage. While spring-loaded backpacks are relatively new, individuals in Asia have been using flexible bamboo poles to transport loads for centuries. However, the potential benefits of carrying an oscillating load are still under debate. Some studies have shown that these systems can increase energy expenditure (Foissac et al., 2009; Martin and Li, 2018) or have little effect (Kram, 1991) compared with carrying loads in a conventional backpack. These discrepancies may be explained by key differences between the various studies – for example, gait type [walking (Rome et al., 2005, 2006; Castillo et al., 2014; Foissac et al., 2009; Martin and Li, 2018) versus running (Kram, 1991)], directionality of load oscillations [vertical (Rome et al., 2005, 2006; Castillo et al., 2014; Foissac et al., 2009; Kram, 1991) versus medial–lateral (Martin and Li, 2018)] or load suspension mechanism [spring-loaded backpack (Rome et al., 2005, 2006; Foissac et al., 2009; Martin and Li, 2018) versus compliant poles (Castillo et al., 2014; Kram, 1991)].

Dynamics models of compliant loading have demonstrated the effect of spring constant and damping on the energetic cost of human walking (Ackerman and Seipel, 2014; Li et al., 2016a). Such models assume leg length changes over stance (prescribed as a sinusoidal path) and parameter values (e.g. frequency, amplitude) scaled to empirical data but neglect an individual's capacity to adjust gait and optimize interactions with the load. A human walking with an oscillating load is essentially a coupled oscillator system, with one oscillator passive (the flexing load) and the other active (the person). An individual may experience and adapt to interactive effects, particularly if there is an advantage to doing so. Such interactions are observed in various systems, e.g. while walking on a swaying pedestrian bridge, individuals spontaneously synchronize their steps with one another and entrain to the bridge's oscillations (Dallard et al., 2001). We hypothesize entrainment exploits energy-saving opportunities available when gait is coordinated with substrate motion. People also modify behaviour to wearable devices; infants bounce in an elastic harness and learn to match frequency to the system's resonance (Goldfield et al., 1993). Individuals wearing a knee exoskeleton optimize step frequency over a range of available energetic cost options (i.e. cost ‘landscape’) when naturally preferred frequencies are penalized with controlled resistance from the mechanical device (Selinger et al., 2015). These examples describe unusual circumstances that nonetheless illicit coherent gait adaptations over time. Modelling approaches that rely on fitting data from standard locomotion tasks (e.g. undisturbed treadmill walking) often fall short in predicting, and thus explaining, non-standard locomotion tasks (e.g. carrying load oscillating on a compliant bamboo pole).

Trajectory optimization is an alternative that has made accurate predictions about energy-minimizing gait solutions under a variety of standard and non-standard circumstances: walking/running at various step length–velocity combinations (Srinivasan and Ruina, 2006), uphill/downhill locomotion (Hasaneini et al., 2013), navigating a shaking platform (Joshi and Srinivasan, 2015), walking with an oscillating impulse applied to the body (Schroeder and Bertram, 2018), and even manoeuvres associated with the urban sport of parkour (Croft et al., 2018). Here, we used trajectory optimization to determine energy-minimizing solutions for humans carrying loads suspended from flexible poles of various stiffness and load. We then compared these results with empirical data from Vietnamese farmworkers highly experienced in this manner of load carriage. Perhaps these individuals are sensitized to functional advantages of this carrying mode and/or are capable of subtle adjustments that exploit those opportunities. Experience can profoundly influence the energetics of unusual load carrying. For example, women in some African tribes can carry loads on their heads for substantially less energetic cost than novice westerners carrying the same load in a backpack or on their heads (Maloiy et al., 1986; Cavagna et al., 2002). When inexperienced individuals are introduced to novel cost landscapes, they sometimes require guided ‘exploration’ where a range of gait options are introduced before they can spontaneously locate the optimal solution (Selinger et al., 2015). When carrying loaded bamboo poles, unguided exploration through natural variation probably occurs with extensive practice – particularly when exposed to various loads, poles and terrains. Past studies are largely limited by their exclusive use of naive, inexperienced participants, which may help explain inconsistent findings. Thus, we studied Vietnamese farmworkers with years of experience carrying substantial loads on bamboo poles.

Fourteen participants (eight males, six females) were recruited for the study after a local translator conducted a brief pre-participation interview to ensure that basic inclusion criteria were met (e.g. no recent injuries affecting gait and at least 5 years of experience, although many participants claimed lifelong experience). Of the 14 participants, seven used a pole daily, three weekly, one monthly, and three seasonally. At the time of data collection, the mean (±s.d.) age was 45±17 years, body mass was 50.2±5.9 kg and height was 1.56±0.07 m. All participants provided informed consent to participate through a qualified interpreter, and these studies were approved by ethics review boards of Thái Nguyên University of Medicine and Pharmacy, Edith Cowan University (Human Ethics Review Board 15249) and the University of Calgary (REB16-0910).

Participants were asked to walk along a 1 m-wide path with a steady gait while carrying one of two pole types: (i) a rigid, season-dried bamboo pole with a full circular diameter (i.e. it was not split to make the characteristic flattened carrying pole, so remained rigid) or (ii) their own personal compliant bamboo pole (pre-fabricated and well used) – property measurements of these poles were presented previously (Schroeder et al., 2018). Participants walked along a 20 m path at a steady preferred velocity and cadence while carrying loads (0%, 30% and 50% body weight, BW). For zero load, subjects carried their own pole unloaded on their shoulder in a casual manner. Otherwise, sandbags were used to load the pole to a proportion of the participant's body weight. The sandbags were suspended from the pole in carrying baskets on wire frames (baskets and frames purchased from a local market in the region) and were evenly distributed between the two ends of the pole in the standard manner employed in the region. The order of conditions (pole type and load) was randomized with five trial repetitions each. Participants were allowed practice trials to acclimate to each condition as needed. All instructions were communicated through a local guide and translator.

Wireless inertial measurement units (Xsens Technologies B.V., Enschede, The Netherlands) captured acceleration and angular displacement in the vertical, lateral and fore–aft directions (acquisition rate=60 Hz). Three sensors were placed on subjects (above each ankle and at the lower lumbar region of the back) and three along the pole's length (directly above the shoulder contact point and at each end near the attachments of the loaded baskets). The ankle sensors were used to determine step frequency by calculating the time difference between acceleration peaks, and the remaining sensors were integrated to characterize fluctuations in centre of mass (CoM) and pole–load velocity and displacement over each step. All sensor signals were nulled between each trial to minimize drift error.

Although reaction forces between the pole and shoulder were not measured directly, an approximation was calculated from kinematics. Compliant pole spring force (F_s) was calculated via the product of the pole spring constant (k_p) and the vertical displacement between the body CoM (y_c) and the load (y_L): reaction force of the compliant pole R_cmp≈F_s=k_p(y_c−y_L). Damping forces were neglected given the system's distinctly underdamped nature (Schroeder et al., 2018). For the rigid pole, reaction force was calculated as: ⁠, where m_L is the load mass, ÿ_c is vertical acceleration of the body CoM, and g is gravitational acceleration (9.81 m s⁻²), assuming motion of the load and body are tightly coupled.

A GoPro Hero 4 camera (GoPro Inc., San Mateo, CA, USA) recorded video in the sagittal and coronal planes (frame rate=120 Hz), although only the sagittal data were analysed. The sagittal camera was placed perpendicular to and approximately 4 m away from the pathway to reduce parallax and lens distortion. Video calibrations were performed between participants by filming two markers at known vertical and horizontal distances in the sagittal plane. Average forward velocity was determined from the videos by converting pixels travelled to metres travelled and dividing by the time duration. A rough indication of leg length was determined by approximating a point at the subject's greater trochanter and measuring the vertical distance to the ground while the subject stood still. This method was used instead of a measuring tape to avoid potential cultural sensitivities as well as any misunderstanding due to language barriers. All video digitization was performed in MATLAB (MathWorks, Natick, MA, USA) using custom DLTdv5 software (Hedrick, 2008).

where f_s is the participant's absolute step frequency. Relative step frequency has a profound influence on the magnitude and phase of the load's oscillation relative to the individual carrying the pole. The response is particularly transient near resonance (i.e. f_r=1), where relative magnitude spikes and phase shifts dramatically from 0° (in phase) to 180° (out of phase, Fig. 1; comparable to analysis in Castillo et al., 2014).

These foot contact positions define the initial position of the body CoM to begin at the middle of double stance (i.e. dual contact phase) at t=0s.

The FRS term is used as a smoothing factor to penalize extremely impulsive forces that are physiologically unrealistic and numerically challenging. Although the use of this cost is somewhat arbitrary, there is emerging evidence that metabolic consumption increases with frequency of muscle activation (and force rate) even as mechanical work remains constant (Doke and Kuo, 2007). Other gait optimization models have similarly incorporated force-rate terms in their cost functions (Rebula and Kuo, 2015; Handford and Srinivasan, 2018). Here, the scaling constant was adjusted until CoM acceleration peaks matched those from participants carrying zero load and was then held constant throughout all optimizations (Fig. 2B).

where v is average forward velocity. Note that the CoT is normalized to body mass instead of total system mass (body CoM plus load) to distinguish the effect of added cost due to increasing loads.

During the optimization, all variables were non-dimensionalized with the parameters L_max, g and m_t; however, all optimization outputs were re-dimensionalized with the appropriate parameter values, as needed. Subject parameter values such as load mass, body mass, maximum allowable leg length and average forward velocity were used to represent an average study participant as well as individual participants. A summary of these parameter values is displayed in Table S1.

The trajectory optimization procedure was implemented in MATLAB using a sparse non-linear optimizer program (SNOPT; Gill et al., 2005) in conjunction with GPOPS-II (Patterson and Rao, 2014) for problem discretization and setup. In order to procure robust solutions, a two-part optimization regime was used (Schroeder and Bertram, 2018). The first part implemented 15 random initial guesses to test for global optimality, and the second perturbed the prevailing optimum 15 times with random noise in order to fine tune the solution's local optimality. This procedure was employed to predict CoT of multiple optimal solutions sweeping a wide range of relevant parameter values. Specifically, step frequency (f_s) and damped resonant frequency (via spring constant) of the pole–load (f_DR) were systematically varied to probe the cost landscape as a function of relative step frequency (recall, f_r=f_s/f_DR).

where A_r is the relative amplitude and φ is the phase of the load amplitude relative to the forcing function.

The damping ratio resulting from the regression (ζ=0.172) was used to simulate optimal load interactions and to determine local CoT gradients for individual participants at and around their average preferred step frequencies (f_s,avg±0.20 Hz, 0.02 Hz intervals) for both 30% and 50% loading conditions. This was done to relate the local cost gradient to changes in average relative step frequency between rigid and compliant poles.

A Pearson correlation was performed to test the relationship between the change in relative step frequency and the local gradient of the CoT. We hypothesized that a negative correlation would be found, indicating that participants adjust step frequency to reduce CoT while carrying loads with a compliant pole. This hypothesis is consistent with three specific outcomes: (i) positive CoT slopes associated with negative changes in relative step frequency, (ii) negative CoT slopes associated with positive changes in relative step frequency, and (iii) zero CoT slopes associated with little or no change in relative step frequency. Least squares linear regressions were performed on the data to indicate what proportion of variance could be attributed to the linear relationship.

Linear mixed models were also used to more thoroughly evaluate changes in step frequency due to fixed effects such as loading level and pole type. The mixed model was chosen to control for repeated measures observed among different subjects, where subject was entered as a random effect. The statistical models were developed in JMP (SAS Institute Inc., Cary, NC, USA, version 14.1.0) using the restricted maximum likelihood method for parameter estimation and a compound symmetric covariance structure. First, changes in absolute step frequency were evaluated due to load effects with the rigid pole as well as the no-load condition (0%, 30% and 50% BW). As walking speed has a known effect on step frequency it was non-dimensionalized [ṽ=v/√(gL_max)] and included in the model as a covariate. The relationship between frequency and speed is generally non-linear (f_s∝v^0.58; Kuo, 2001; Bertram and Ruina, 2001). However, a simple linear regression was performed for step frequency versus non-dimensional velocity, and the relationship was deemed sufficiently linear after assessing residuals over the relatively narrow range of speeds chosen by participants.

Another model was used to evaluate changes in relative step frequency, where pole type (rigid versus compliant) and load (30% and 50% BW) were both included as fixed effects. Non-dimensional walking speed was again included in this model as a covariate after confirming linearity with regression and visual assessment of the residuals. Interaction terms relating pole type to load, pole type to ṽ and load to ṽ were also included in the model.

Our optimization model predicts differential subject responses in relative step frequency when carrying the compliant pole versus the rigid pole. Thus, differences in relative step frequency based on the least squares means of each pole type were evaluated for individual subjects with post hoc t-tests and 95% confidence intervals (CI). The significance of our mixed model effects was also evaluated with post hoc t-tests. We adjusted P-values (P_adj) based on Benjamini and Hochberg's (1995) method to control the false discovery rate during multiple significance testing and considered tests with P_adj<0.05 to be significant. Throughout the manuscript, unadjusted P-values are reported, and significance is indicated with asterisks. Additional details of the statistical models and general approach can be found in the Appendix.

During unloaded walking, there is a trade-off between two cost mechanisms that influence the model: work to swing the leg (costly at high frequencies) and work to extend the support legs (costly at low frequencies). Thus, there exists an optimal intermediate step frequency (1.90 Hz; Fig. 3A). This minimum cost solution agrees well with the least squares mean step frequency of participants walking with no load [mean (95% CI): 1.89 Hz (1.83–1.96 Hz)]. The optimization model predicts a slight increase in the optimal step frequency to carry non-zero loads, as leg extension work increases (further penalizing low frequencies) while the cost of leg swing remains unchanged. The model's optimal frequency increases to 1.98 Hz (a 4.21% increase) and to 2.02 Hz (a 6.32% increase) with 30% and 50% BW loads, respectively. These compare to mean participant step frequencies of 1.94 Hz (1.88–2.01 Hz) (a 2.65% increase, P<0.001) and 2.02 Hz (1.95–2.08 Hz) (a 6.88% increase, P<0.001*; Fig. 3B).

Optimized model outputs were compared with average trials of example subjects walking with varying relative step frequencies and a 50% BW load carried on compliant poles (Fig. 4). Vertical CoM and load position were normalized by participant leg length and plotted alongside outputs of the optimization model used to simulate circumstances appropriate to subject- and condition-specific parameters. The CoM's average height was adjusted to match that of the model curves, and the load's average height was determined such that average spring force matched the weight of the load. Spring force in Fig. 4 was normalized by body weight. In the examples shown, the model qualitatively matched the empirical trends reasonably well. Phase between the CoM and load positions gradually shifted from very low values (<10°) at f_r=0.76 (Fig. 4A) to ∼84° phase at f_r=1.12 (Fig. 4D). A more comprehensive comparison between model outputs and participant data can be found in Fig. S1.

Even with constant step frequency, cost still fluctuates according to the spring constant – and thus resonance – of the flexible loaded pole; the oscillation's magnitude, phase and resulting reaction forces transferred to the body are all greatly affected by relative step frequency. Fig. 5A shows how cost changes when the model takes the same step over a large range of pole spring constants (i.e. absolute step frequency stays constant, but relative step frequency changes as a result of resonance). The curves in Fig. 5A show cost for 50% BW loading of an average participant carrying a compliant pole, with low damping from Schroeder et al. (2018) and higher damping from the fit in Fig. 1. Similar trends were found for 30% BW loading, although the cost fluctuations were more subtle.

The global minimum cost occurs where relative step frequency is slightly above resonance (f_r=1.08 for low damping or f_r=1.30 for high damping) and the maximum is just below (f_r=0.90 for low damping or f_r=0.85 for high damping; Fig. 5A). In the low damping curve there is a local minimum that occurs at a relative step frequency of 0.60, just above a 2:1 harmonic, where the load exhibits low magnitude oscillations at twice step frequency (Fig. 5C). A 3:1 harmonic is observed where the load oscillates at three times step frequency (Fig. 5B). Although the harmonic frequency oscillations are too small to see in the position plot, their effect is visible in the spring force plot. Generally, out-of-phase load oscillations (90°≤φ≤180° at f_r>1) are associated with a lower cost of transport (CoT_cmp) relative to carrying a rigid load (CoT_rig). Conversely, in-phase load oscillations (0°≤φ≤90° at f_r<1) are associated with an increased CoT_cmp; however, this increase is slightly reduced at frequencies just above the 2:1 and 3:1 harmonic points (f_r=0.60, 0.40, respectively).

Relative step frequency is the ratio of absolute step frequency to resonant frequency of the pole–load system. In Fig. 3B, resonant frequency was constant (rigid pole) and the effect of step frequency on cost was explored. In Fig. 5A, step frequency was constant and the effect of resonant frequency on cost was explored. However, in practice, the two parameters have a simultaneous effect, as participants can choose both their step frequency and the pole they carry. In particular, spring constant varied from pole to pole, and this probably had a pronounced influence on the subjects' cost landscapes. Fig. 6 shows how poles with three spring constants (lowest, average and highest spring constant of poles in our sample) can have a dramatic effect on the shape of the total cost landscape.

Each panel in Fig. 6 has three cost curves. The blue curve indicates cost for carrying a compliant pole over a large range of relative step frequencies. This cost is identical in all three panels as frequency is normalized to pole resonance. The red curve indicates cost for carrying a rigid pole over various step frequencies (1.60 Hz≤f_s≤2.88 Hz). This cost is also identical in the three panels, although its relative step frequency domain is stretched depending on the resonance of the chosen spring constant. Finally, the black curve indicates how the first two costs interact to influence the total cost curve for an individual with a given pole spring constant and range of step frequencies. Even though the red and blue curves do not change cost in each of the three panels, the shape of the black curve is completely different as the red curve occurs over different segments of the blue curve and this influences the final cost. It is inappropriate to simply sum the two cost curves (i.e. red plus blue equals black), because the cost mechanisms are interdependent and may be optimized simultaneously. However, they can still be conceptualized as two influences that help to shape the total cost surface that a pole carrier navigates.

Individuals navigate different cost landscapes, not only as a result of differences in body mass and morphology, but also because of the properties of the pole they carry. As such, participants may adapt differently depending on the pole, even if energy minimization is the overarching optimization goal.

The effect of pole type and load on relative step frequency was statistically assessed with linear mixed models. Participants as a whole exhibited a slight but significant increase in relative step frequency when switching from the rigid pole to the compliant pole [coefficient (95% CI): β=0.014 (0.004–0.023), P=0.014*]. However, when non-dimensional walking speed (ṽ) was controlled for, this effect was diminished [β=0.009 (−0.001–0.019), P=0.083]. Additionally, one interaction term had a significant effect on frequency (pole type×ṽ) while two others did not (load×ṽ and pole type×load; see Appendix Table A1 for full model results). Load was also found to have a strong effect on relative frequency. However, this is unsurprising as load decreases resonant frequency which increases relative step frequency. Given our original hypotheses predict participants should respond differently (e.g. increase/decrease/no change to relative step frequency) depending on the local slope of their cost curve, we also tested individual responses to pole type (Appendix Table A2). Shifts in average relative step frequency (Δf_r) were used to quantify gait adaptations of individual participants carrying a compliant pole versus a rigid pole. The direction and magnitude of these shifts are compared with the local gradient, or slope, of each participant's CoT curve in Fig. 7. A significant negative correlation was found for the 50% loading level (R=−0.67, P=0.009*; Fig. 7B), with five subjects exhibiting significant changes in relative step frequency (indicated in Fig. 7B with magenta diamonds). Specifically, three subjects had a significant positive shift (as in Q2 from Fig. 7), two had a significant negative shift (Q3 and Q4 in Fig. 7) and the rest were not significantly different from zero (red diamonds). This result is consistent with our hypothesis that participants adjust relative step frequency to reduce the CoT predicted by our optimization model when carrying a compliant pole. Slightly under half of the variation could be attributed to the simple linear regression shown in Fig. 7B (y=−0.012x−0.002, R²=0.45).

For the 30% loading level, a weak, non-significant correlation was found (R=−0.23, P=0.430; Fig. 7C). Just over 5% of the variation in the data could be attributed to the simple linear regression (y=−0.004x−0.010, R²=0.053). Notably, changes in relative step frequency were less pronounced than in the 50% loading condition, even though they were associated with a similar range of CoT gradients, and even though two subjects exhibited significant positive shifts in their relative step frequency (magenta diamonds in Fig. 7C; Q1 and Q2 from Fig. 7A). The average magnitudes of the frequency changes were 0.020 and 0.029 for the 30% and 50% loading conditions; however, the 30% condition was largely influenced by an outlier with a particularly large shift (0.119). Without this outlier, the average shift in relative step frequency was only 0.012 for the 30% loading level. For both loading levels, the y-intercepts were not significantly different from zero (95% CI: −0.020–0.017 for 50% BW loading, −0.012–0.032 for 30% BW loading). These results are consistent with our hypothesis that a zero CoT slope should not motivate a shift in relative step frequency.

The optimization model predicts many attributes of human gait, e.g. the mass-specific CoT of unloaded walking varies over absolute step frequency as a convex function (i.e. ‘bowl shape’; Bertram, 2005). However, the minimum CoT (∼3 J kg⁻¹ m⁻¹) is somewhat higher than that of human metabolic data during walking at a similar velocity and step frequency (closer to 2 J kg⁻¹ m⁻¹; Bertram, 2005; Bastien et al., 2005). Although reasons behind the cost shift are unclear, the general trends are qualitatively consistent.

In our model, high step frequencies were dominated by a cost for swinging the leg and low frequencies by a cost to extend the leg (Fig. 3A). Previous studies have shown that average mechanical leg power and work increase proportionately to step length raised to the fourth power (Kuo, 2002; Donelan et al., 2002) for constant step frequency. Such models also predict leg work increases non-linearly at lower step frequencies, similar to our model.

Models that consider the cost of leg work and/or frequency-based cost functions for leg swing have been used to explore issues such as the velocity−step length relationship in walking (Kuo, 2001), the cost of swinging the leg isolated from gait (Doke et al., 2005), and running in reduced gravity simulations (Polet et al., 2018). In the former two, a force–rate cost has been proposed as a major determinant of energetic cost at higher frequencies, while the latter study indicated that a work-based cost gave more accurate predictions of kinematic outputs. Although more understanding of muscle contraction cost is needed, we opted for a work-based cost of swinging the leg simply so results could be compared more readily with leg extension work. Furthermore, force rate and work costs give similar predictions at moderate non-dimensional velocities: ∼0.38 (Kuo, 2001) and ∼0.43 in the current study.

The model predicts CoT for carrying rigid loads. Previous studies have shown that cost increases linearly with load carried (Griffin et al., 2003; Bastien et al., 2005), comparable to the current model (Fig. 3B). The optimal frequency is also predicted to increase, as leg extension work increases (i.e. higher penalty on low step frequencies) while leg swinging work is assumed constant. Changes in step frequency were statistically validated by our linear mixed model (P<0.001*), and subjects have been shown to increase frequency while carrying large loads in other studies (LaFiandra et al., 2003). However, such changes are not always statistically significant (Castillo et al., 2014). A potential for type II errors may exist because of the subtlety of these effects.

When step parameters (velocity, step frequency, etc.) are held constant, the model predicts minimal CoT at 1.08<f_r<1.30 (depending on damping). This agrees with work by Castillo et al. (2014) showing that large out-of-phase load oscillations (as in Fig. 5F) contribute to a minimum CoT at f_r=1.21. At this frequency, metabolic cost was reduced by ∼5% when carrying a compliant bamboo pole versus a steel pole. The idea that energetic cost should be reduced at relative step frequencies just above one directly contradicts previous studies suggesting a low spring constant is optimal (Kram, 1991; Ackerman and Seipel, 2014), and that step frequencies near resonance are always most costly (Ackerman and Seipel, 2014; Li et al., 2016a). Although our model predicts maximal cost to occur just below resonance, minimal cost is predicted just above resonance, and this implies the tuning of pole properties such as spring constant can play a critical role in facilitating optimal interactions.

Tuning spring constant was also found to be important for a backpack with elastic load suspension in the fore–aft direction (Li et al., 2016b), where out-of-phase oscillations were shown to modestly reduce mechanical power calculated from ground reaction forces. Kram (1991) found no difference in the cost of carrying load on a compliant structure made of polyvinyl chloride (f_r≈3) versus the expected cost of using a backpack with non-compliant straps. This confirms our model's prediction that the cost of carrying a compliant load approaches that of a rigid load at high relative step frequencies (Fig. 5A), although this may be a coincidence given Kram (1991) studied running, not walking.

Castillo et al. (2014) attributed their reduced cost findings to a relatively flat system mass trajectory, where load oscillations cancel out body oscillations. Although our results do not contradict this, we find it more insightful to consider the cost mechanisms in our model. Specifically, leg extension work is largely affected by the loading cycle felt as a force acting on the body. Fig. 5B–I shows the spring force of the pole transferred to the shoulder (solid green curves); damping forces are neglected in the underdamped system. These loading cycles indicate periods with relative on-loading and off-loading, where on-loading refers to spring forces more negative than the load weight, effectively heavier (i.e. solid green line below the dashed green line; Fig. 5B–I), and off-loading refers to spring forces less negative than the load weight, effectively lighter (i.e. solid green line above the dashed green line; Fig. 5B–I).

Generally, when load oscillations are out of phase, on-loading occurs during single stance and off-loading during double stance (f_r>1; Fig. 5F–H), resulting in a reduced CoT (relative to carrying a rigid load). Conversely, when load oscillations are in phase, off-loading occurs during single stance and on-loading during double stance (f_r<1; Fig. 5B–D), resulting in an increased CoT. A simple explanation concerns how leg extension work accrues over a step (leg swing cost is unchanged throughout Fig. 5A as step frequency is constant). The absolute value of leg extension power is concentrated near double stance, regardless of relative step frequency. This is consistent with the inverted pendulum description of walking – little external mechanical work occurs during single stance and the body moves over the foot as an inverted pendulum with the leg acting as a strut (Cavagna and Margaria, 1966).

Energy loss occurs at the transition from one inverted pendulum to the next and positive leg work is done to make up the loss and maintain steady gait (Donelan et al., 2002; Kuo, 2001; Lee et al., 2013). As most external work is done during double stance, leg extension cost is sensitive to loading during this transition. This helps explain why the in-phase/out-of-phase relationship is so well correlated to cost. When on-loading occurs during double stance, the leg bears the extra load while doing work on the system mass, and this is energetically expensive. When on-loading occurs during single stance, the leg simply bears the extra load isometrically at zero work cost (there is still a force rate accommodation, but this cost is small).

For higher damping, cost is minimized even though load oscillations are only partially out of phase (ɸ=135°, Fig. 5I). This is because fully out-of-phase oscillations at higher damping require a much higher relative frequency (Fig. 1B), far removed from resonance where the amplitude is dwindling. As a result, the minimum cost solution compromises on a relative frequency that is high enough to move the oscillations somewhat out of phase but low enough (i.e. closer to resonance) to garner meaningful amplitude such that off-loading during double stance can still affect cost. Although participants in our study carried loads with oscillations mostly in-phase (Fig. 4; Fig. S1), positive shifts in relative step frequency probably increased phase away from 0°, thus reducing the amount of off-loading during double stance and contributing to a lower cost (negative cost slope: magenta diamonds in Q2 of Fig. 7B,C). It is unclear why subjects did not increase frequency more in order to fully converge on optimal phase relationships. Other studies have suggested physiological noise as a limitation on individuals learning gait parameters minimizing energetic cost (Simha et al., 2019).

There are a few surprising nuances that arise from the low damping solutions in Fig. 5. Although an in-phase relationship is observed for most relative step frequencies below one, a local minimum occurs just above the 2:1 harmonic frequency (f_r=0.60, Fig. 5C). This minimum occurs because the harmonic frequency cancels out some of the on-loading felt by the in-phase, fundamental oscillation. In Fig. 5B, a relative step frequency slightly above the 3:1 harmonic occurs (f_r=0.40). In this case, both the step frequency and the harmonic frequency are in phase with the body. The harmonics have less influence on overall cost as they occur further away from the fundamental resonant frequency, and magnitude diminishes. However, they provide an interesting subtlety to the interpretations of load oscillations and consequences for associated cost mechanisms.

Another interesting result can be found at the maximum cost point (f_r=0.90; Fig. 5E). The oscillations of the CoM and load appear quite flat, despite their vicinity to resonance. In fact, the load amplitude is quite large (∼3.3 times that of the body; see Fig. 1A). At the same time, leg extension power continues well into the single stance portion of the step (unlike other solutions shown). This is consistent with a gait sometimes called ‘Groucho walking’ (Bertram et al., 2002), characterized by a flat body trajectory and compliant legs. Groucho walking is energetically costly as a result of work done to excessively flex/extend the leg and reduce vertical body oscillations (Ortega and Farley, 2005; Gordon et al., 2009; Kim and Bertram, 2018). Although this solution indicates maximal cost on the curve, it is energetically minimal for the particular relative step frequency, so this is the strategy selected by the optimization model. Our interpretation is that in-phase load oscillations near resonance are so costly for leg extension work that it is cheaper to employ a flat, Groucho-like gait to avoid the excessive oscillations. A more modest example of this strategy is observed at f_r=0.76 (Fig. 5D), providing evidence of a trade-off between costly load oscillations just below resonance and costly leg extension work associated with a flat gait. Foissac et al. (2009) observed similar gait in subjects carrying a flexible backpack load, where vertical trunk excursion was reduced at walking speeds stimulating resonance and energetic cost was increased.

Although Groucho walking is one way to avoid costly large-amplitude in-phase oscillations near resonance, another strategy is to simply reduce relative step frequency. Even though such adjustments bring phase even closer to zero, they reduce oscillation amplitudes mostly in phase anyway and this reduces cost (positive cost slope: f_r<0.9, Fig. 5A). It appears that most subjects did not use this strategy to reduce cost; however, one subject significantly reduced relative step frequency on a positive cost slope (Fig. 7B, magenta diamond in Q4).

Relative step frequency shifts at local cost gradients are summarized in Fig. 7. A significant correlation was found for 50% BW loading, but not for 30%. Given that CoT is generally higher for increased loading levels, it is possible that this plays a role in the sensitivity of individuals to adapt gait. However, our hypothesis was that participants follow the gradient downhill, regardless of cost. In other words, if there are cost savings available, why do experienced pole carriers not take advantage of them at 30% BW load?

It is likely that individuals – even experienced pole carriers – are only sensitive to cost savings at a certain threshold. Indeed, many participants walked at relatively flat cost gradients and did not adjust frequency much. It is unclear whether participants recognized there were no lower cost solutions nearby, or whether they were simply insensitive to shallow gradients. There are also examples of large CoT gradients where participants exhibited little response. The simple regression slope can be conceptualized as the sensitivity of the participant sample to their CoT gradient. However, it is likely that individuals have varying sensitivities to such gradients, characterized by individual slope values. A study with larger load levels or longer carrying durations might have clearer results. We were reluctant to overload participants even though many were used to carrying much greater loads (sometimes their own body weight or more) for substantial distances (kilometres). Indeed, for an average participant carrying a 100% BW load, our model predicts energy savings of approximately 18% (3.92 J kg⁻¹ m⁻¹ with a compliant pole versus 4.79 J kg⁻¹ m⁻¹ with a rigid pole).

Alternatively, the model may not account for all relevant costs associated with carrying a compliant bamboo pole during walking, such as costs to steady swaying load baskets or for balancing the pole on the shoulder (Li et al., 2019). We assumed experienced pole carriers are adept and that such costs are secondary to the work-based costs implemented in the model. Using a rigid pole for comparison to the compliant pole, as opposed to a load rigidly fixed to the body (backpack), helped account for such potential extra costs.

Experienced pole carriers are probably not only trying to optimize energetic costs. Individuals might choose to carry loads with a compliant suspension system to reduce peak reaction forces felt at the shoulder (Rome et al., 2005, 2006; Kram, 1991). Such reaction forces were approximated with kinematic data from our participants. On average, participants did see a reduction in peak shoulder forces when carrying the compliant pole versus the rigid pole (for 50% loading: 18.0%, or −75 N; for 30% loading: 12.7%, or −31 N). It is unclear how much these changes were influenced by shifts in relative step frequency. Our model predicts changes in the peak shoulder reaction force due to shifts in relative step frequency alone. For the case of 50% BW loading, an average reduction of −4.5 N (1.2%) was found; however, this average was largely dominated by two participants, without whom a slight increase in the reaction force was found of 0.4 N. For the case of 30% BW loading, only a −0.6 N (0.3%) reduction was found. Although shifts in relative step frequency had little influence on this effect for either load condition, participants probably experienced a sizeable reduction in peak shoulder forces when switching from a rigid to a compliant pole.

We hypothesized that changes in step frequency correlate with theoretical cost gradients predicted by a work-minimizing model. However, there are many reasons why step frequency changes occur. Stability is one alternative, as step frequency changes are often associated with destabilizing perturbations (McAndrew et al., 2010; Hak et al., 2012). Recent studies exploring coupled-oscillator models have shown that walking stability is decreased with lower spring constants (Ackerman and Seipel, 2011). At the same time, it is unclear what strategies, if any, experienced pole carriers use to stabilize load oscillations. The natural carrying style of resting a hand on top of the pole may be an effective stabilizing strategy, thus eliminating the need to adjust step frequency to avoid instability. Although the current study did not formally evaluate stability, there were no obvious examples of instability during trials. Further empirical studies are needed to investigate this complicating issue.

Another gait feature highly associated with changes in step frequency is walking speed. One goal of the study was to measure experienced bamboo pole carriers in a more ecological setting (i.e. not artificially in a laboratory); as such, subjects were not constrained to walk at fixed speeds, etc. However, we did attempt to control for speed variation in our statistical models. Non-dimensional walking speed (ṽ) had a strong significant effect on relative step frequency [β=0.447 (0.253–0.640), P<0.001*], and the inclusion of this covariate diminished the effect of pole type on frequency. However, an interaction between pole type and ṽ was also found to have a strong effect [β=0.285 (0.147–0.424), P<0.001*]. Our interpretation is that individuals vary walking speed as a means to economically change relative step frequency when switching pole type from rigid to compliant. In other words, it is less costly to increase step frequency at higher speeds than it is to do so at a constant speed (and vice versa). This is supported by studies showing that the metabolic cost of walking is minimized when speed co-varies with frequency for a given task constraint (Bertram, 2005).

There are numerous simplifying assumptions that may affect the results in this study. Morphology complexities were neglected in favour of a reductionist model. A telescopic leg replaced complicated flexion/extension of the knee and ankle during walking (e.g. soleus and gastrocnemius contraction toward the end of stance; Winter, 1991). Hip actuators could also have applied torque during stance and swing phases. This cost was somewhat accounted for in the model proposed by Doke et al. (2005) where torque work was derived from pendular motions (see Materials and Methods for details). Other simplifying factors include a point mass body (zero mass moment of inertia), no distribution of the load along the length of the pole, non-slip foot contacts and planar sagittal motion only. Although each of these issues has the potential to influence system dynamics and its energetic cost, there are no obvious indications that such considerations are vital to the topic. More complicated models can be used to probe such issues in future studies.

A practical limitation was access to equipment in the field. Metabolic measurements of energy consumption would provide more robust evidence of cost reductions supporting our hypotheses. Without direct measurements, all inferences regarding energy cost are theoretical. Although a primary goal of the study was to allow experienced pole carriers to walk ‘naturally’ (i.e. minimal experimental constraints), this limited the capacity to control for confounding variables and other complications. Future experiments could test learning strategies with more clear cost incentives in a controlled protocol.

We developed a trajectory optimization model to determine the theoretical CoT for individuals carrying compliant and rigid bamboo poles. We used the model to explore energetic consequences of carrying a rigid pole with various loads over a range of step frequencies. We also explored the energetic consequences of carrying poles with varying spring constants and relative step frequencies. The model considers costs due to leg extension work, leg swing work, and a FRS cost. This mechanistic perspective was used to interpret reduced costs associated with the alignment of force on-loading and off-loading the body at different portions of the gait cycle. As the majority of leg extension power is performed during double stance, it is beneficial to be off-loaded during this time, and this occurs most notably in a range of relative step frequencies slightly above resonance where the load and CoM oscillate out-of-phase. Higher order harmonics also affect the CoT, although these effects are modest. Finally, the model predicted changes in relative step frequency as a means of reducing energetic costs associated with the task. Pearson correlations revealed a significant negative correlation between the change in relative step frequency and the local CoT gradient with 50% BW loads. A weak, non-significant negative correlation was also found for 30% BW loads.

Ultimately, direct evidence of gait adaptation associated with empirical reductions in metabolic cost are required. Such experiments will provide additional validation to the cost mechanisms identified here. Regardless, a theoretical framework for understanding optimal body–pole–load interactions has been proposed that can explain various aspects of gait adaptation during load carriage with a flexible pole.

This section details results from the statistical models referred to in the Materials and Methods section.

Table A1 summarizes the repeated measures mixed linear models used during analysis in the manuscript. Model 1 includes data from individuals carrying the rigid pole as well as data from the no-load condition. Absolute step frequency was fitted to load level (0%, 30% and 50% subject BW) while controlling for non-dimensional walking speed (ṽ). Interaction terms between the covariate and the fixed effect were not included as they were not found to be significant. In model 2, the main fixed effects – load level (30% and 50% BW) and pole type (rigid and compliant) – were used to fit relative step frequency data for all conditions except the no-load condition. Non-dimensional walking speed was not included as interactions were found to be significant. However, an expanded version of the model is labelled (model 3) in Table A1. It includes non-dimensional walking speed and related interactions. All three models included subject as a random variable with compound symmetric covariance structures.

Table A2 summarizes differences in the least squares means of relative step frequency between the rigid and compliant pole conditions for individual subjects at both loading levels (30% and 50% BW) using model 2 from Table A1.

We appreciate the cooperation of the Mei Yen Commune, the officials administering Thái Nguyên Province, Vietnam, and the faculty and administration of Thái Nguyên University of Medicine and Pharmacy. We are especially indebted to the farmworkers who volunteered their time to assist with this study. We would also like to thank Andrew Biewener and Robert Dudley for their thoughtful comments and advice on an earlier version of this manuscript.