Skip to main content
Advertisement
  • Loading metrics

Significance of trends in gait dynamics

  • Klaudia Kozlowska ,

    Contributed equally to this work with: Klaudia Kozlowska, Miroslaw Latka, Bruce J. West

    Roles Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft

    Affiliation Department of Biomedical Engineering, Wroclaw University of Science and Technology, Wroclaw, Poland

  • Miroslaw Latka ,

    Contributed equally to this work with: Klaudia Kozlowska, Miroslaw Latka, Bruce J. West

    Roles Conceptualization, Formal analysis, Investigation, Supervision, Validation, Writing – original draft

    Miroslaw.Latka@pwr.edu.pl

    Affiliation Department of Biomedical Engineering, Wroclaw University of Science and Technology, Wroclaw, Poland

  • Bruce J. West

    Contributed equally to this work with: Klaudia Kozlowska, Miroslaw Latka, Bruce J. West

    Roles Formal analysis, Methodology, Writing – review & editing

    Affiliation Office of the Director, Army Research Office, Research Triangle Park, USA

Abstract

Trends in time series generated by physiological control systems are ubiquitous. Determining whether trends arise from intrinsic system dynamics or originate outside of the system is a fundamental problem of fractal series analysis. In the latter case, it is necessary to filter out the trends before attempting to quantify correlations in the noise (residuals). For over two decades, detrended fluctuation analysis (DFA) has been used to calculate scaling exponents of stride time (ST), stride length (SL), and stride speed (SS) of human gait. Herein, rather than relying on the very specific form of detrending characteristic of DFA, we adopt Multivariate Adaptive Regression Splines (MARS) to explicitly determine trends in spatio-temporal gait parameters during treadmill walking. Then, we use the madogram estimator to calculate the scaling exponent of the corresponding MARS residuals. The durations of ST and SL trends are determined to be independent of treadmill speed and have distributions with exponential tails. At all speeds considered, the trends of ST and SL are strongly correlated and are statistically independent of their corresponding residuals. The averages of scaling exponents of ST and SL MARS residuals are slightly smaller than 0.5. Thus, contrary to the interpretation prevalent in the literature, the statistical properties of ST and SL time series originate from the superposition of large scale trends and small scale fluctuations. We show that trends serve as the control manifolds about which ST and SL fluctuate. Moreover, the trend speed, defined as the ratio of instantaneous values of SL and ST trends, is tightly controlled about the treadmill speed. The strong coupling between the ST and SL trends ensures that the concomitant changes of their values correspond to movement along the constant speed goal equivalent manifold as postulated by Dingwell et al. 10.1371/journal.pcbi.1000856.

Author summary

During treadmill walking, the subject’s stride time (ST) and stride length (SL) must yield a stride speed (SS) which can fluctuate over a narrow range centered on the treadmill belt’s speed. The fact that both ST and SL are persistent is an intriguing property of human gait. For persistent fluctuations any deviation from the mean value is likely to be followed by a deviation in the same direction. To trace the origin of such persistence, we used a novel approach to determine trends in spatio-temporal gait parameters. We find that the trends of ST and SL of a subject are strongly correlated and are statistically independent of their corresponding residuals. Moreover, the trend speed, defined as the ratio of instantaneous values of SL and ST trends, is tightly controlled about the treadmill speed. The persistence of gait parameters stems from superposition of large scale trends and small scale fluctuations.

Introduction

Over two decades ago Hausdorff et al. [1, 2] discovered long-range, persistent correlations in stride duration (time) of human gait. Their choice of fractional Brownian motion (FBM) [3] for modelling such correlations has significantly influenced the way in which fluctuations of spatio-temporal gait parameters (SL, ST, and SS) were subsequently quantified and interpreted. Determining the source of ubiquitous trends observed in physical, social, or biological systems is a recurring problem in fractal time series analysis. In other words, one needs to establish whether trends arise from the intrinsic dynamics of a system or have an external origin. In the latter case, a possible approach is first to recognize trends and then filter them out before attempting to quantify correlations in the noise (residuals). In their original study Hausdorff et al. applied detrended fluctuation analysis (DFA) [4] to estimate the scaling (Hurst) exponent which characterizes properties of fractal time series. Their choice of DFA implied that a stride duration time series was made up of persistent (scaling exponent greater than 0.5) fractal fluctuations superposed on trends which are irrelevant from the point of view of fractal analysis.

The papers of Hausdorff et al. spurred significant interest in the emerging field of fractal physiology [5, 6]. In particular, they focused research on the intrinsic variability of physiological time series and their persistence [7]. Some argued that persistent fluctuations are a manifestation of the adaptability of underlying control systems and the absence of long-range correlations, or anti-persistence, indicates disease or pathology [810]. Delignières and Torre demonstrated that when healthy humans walk over ground in sync with a metronome, their stride times become anti-persistent [11]. During treadmill walking, a subject’s stride time (ST) and stride length (SL) must yield a stride speed (SS) that fluctuates over only a narrow range centered on the treadmill belt’s speed. It turns out that subjects’ SS is anti-persistent, while the other two gait parameters, ST and SL, are persistent. On the other hand, all three parameters are anti-persistent during treadmill walking with either auditory [12] or visual [13] cueing (alignment of step lengths with markings on the belt). In light of these latter experimental findings, anti-persistence, rather than being a manifestation of pathology, seems to be indicative of tight control.

While persistence of ST and SL during treadmill walking is intriguing, an even more puzzling property of gait dynamics is the weak coupling between stride time and stride length measured by cross-correlation. The coefficient of correlation between these dynamical variables for uncued treadmill walking is 0.28 and doubles to approximately 0.55 under the influence of a persistent fractal metronome [14]. Herein, we employ an explicit algorithm to determine trends in spatio-temporal gait parameters and examine their statistical properties as well as those of the corresponding residuals. In the process, light is shed on the origin of ST and SL persistence.

Results

ST and SL trends

The examples of stride time, length and speed time series are given in Fig 1. The solid, thick lines in this figure represent trends approximated using the piecewise linear variant of the MARS model.

thumbnail
Fig 1. Stride time, length, and speed of two healthy young subjects during treadmill walking at a preferred walking speed: (A) subject 6, trial 2, (B) subject 1, trial 2.

Thick solid lines show trends calculated using the piecewise linear variant of the MARS model. The thin solid lines correspond to trend speed (the ratio of values of stride length and stride time trends). Experimental data come from the study of Dingwell et al. [15].

https://doi.org/10.1371/journal.pcbi.1007180.g001

The Kruskal-Wallis test showed that the durations of normalized trends (durations divided by the subject’s mean value of stride time or stride length) were independent of treadmill speed (pST = 0.30, pSL = 0.81). Therefore, we aggregated data from all the trials. In Fig 2 we present the probability density functions (PDFs) of 1607 ST and 1396 SL trend durations. The Kolmogorov-Smirnov test showed that the distributions of both gait parameters were not different (p = 0.81). The average values and standard deviations of trend durations were close to each other (23.2 and 26.3 for ST, and 25.2 and 35.5 for SL).

thumbnail
Fig 2. Probability density function (PDF) of normalized duration of trends for: (A) stride time and (B) stride length.

Trend durations turned out to be independent of treadmill speed. Consequently, we created two aggregated data sets from trials at different speeds. The solid curves in both subplots show the best-fit exponential PDFs. The probability plots, presented as insets, show the differences between the experimental and exponential cumulative distribution functions (CDFs).

https://doi.org/10.1371/journal.pcbi.1007180.g002

The fundamental property of the exponential distribution, whose probability density function (PDF) is equal to exp(−λx)λ, is that its average value and standard deviation are both equal to 1/λ. Therefore, in Fig 2 we also present the exponential PDFs fitted to the experimental data. For both ST and SL, λ = 0.043. It is apparent from the probability plots, presented as insets in Fig 2, that the exponential distribution fairly approximates the tail of the experimental distributions.

We defined a normalized slope of a trend as the change of gait parameter normalized by the product of average value of this parameter during the trial and the normalized trend duration. The normalized slopes were independent of the treadmill’s speed. Consequently, for both ST and SL, we merged the data from all trials. In Fig 3 we present the PDFs of normalized ST and SL slopes. While the the Kolmogorov-Smirnov test showed that ST and SL slopes were drawn from the same distribution (p = 0.16), the Anderson-Darling test indicated just the opposite (p = 0.04). The thick lines in this figure are the Cauchy’s PDFs fitted to experimental data. The probability plots indicate that the fits are of good quality. The Kolmogorov-Smirnov test confirmed that the ST (p = 0.11) and SL (p = 0.06) slopes were drawn from a Cauchy distribution. The Anderson-Darling test gave the opposite results (p = 6 × 10−3 in both cases). The scale parameter γ of the Cauchy distribution was equal to 1.2 × 10−3 and 1.5 × 10−3 for ST and SL, respectively. The location parameter (location of the maximum of PDF) was equal to 8 × 10−4 for ST and 6 × 10−4 for SL.

thumbnail
Fig 3. Probability density function (PDF) of normalized slopes of trends in: (A) stride time and (B) stride length.

Slopes turned out to be independent of treadmill speed. Consequently, we created two aggregated data sets from trials at different speeds. The solid curves in both subplots show the best-fit Cauchy PDFs. The probability plots, presented as insets, show the differences between the experimental and Cauchy’s cumulative distribution functions (CDFs).

https://doi.org/10.1371/journal.pcbi.1007180.g003

We also analyzed the properties of MARS trends in phase-randomized SL and ST surrogate data. Both for SL and ST, we aggregated data from all trials and for each time series we generated a cross-correlated SL-ST surrogate (CCS) and an independently randomized surrogate (IS). For such surrogates, we built the MARS models. Then, using the Kolmogorov-Smirnov tests we compared the distribution of the normalized trend durations of the experimental data with those of CCS and IS. We repeated such comparisons using n = 50 different realizations of surrogate data sets. The same procedure was performed for the slope of trends.

For SL, the trend duration distribution was not statistically different from those of CCS and IS.

For ST, the CCS and IS distributions were the same as the experimental one in 92% and 46% cases, respectively. The observed differences were not large. For CCS, the mean trend duration was equal to 21.4 ± 0.5, range 20.5 − 22.8. For IS, the mean was equal to 20.9 ± 0.3, range 20.3 − 21.8. Please recall that the experimental data mean was equal to 23.2.

The trend slope distribution for SL and ST surrogate data exhibit similar properties.

Scaling exponents

The mean values of ST, SL, and SS scaling exponents at different treadmill speeds are presented in Table 1 along with the corresponding standard deviations. All indices that are statistically smaller than 0.5 are marked with the dagger symbol. The parameters α(1), α(2), and α(3) are exponents calculated using detrended fluctuation analysis of order 1 (DFA1), 2 (DFA2), and 3 (DFA3), respectively. The parameter α(MD) was computed using the madogram estimator (MD). Scaling analysis was performed for the raw experimental time series and time series detrended using the piecewise linear variant of the MARS model. In the latter case, the exponents have an L subscript. Fig 4 displays the results of the scaling analysis for ST and SL time series.

thumbnail
Fig 4. Scaling exponents of stride time (A) and stride length (B) as a function of treadmill speed.

α(1), α(2), α(3), and α(MD) are exponents calculated using detrended fluctuation analysis of order 1 (DFA1), 2 (DFA2), 3 (DFA3), and madogram estimator, respectively. Calculations were performed for raw experimental time series and the time series detrended using the piecewise linear variant of the MARS model. In the latter case the exponents have an L subscript. The horizontal red line delineates persistent (> 0.5) and anti-persistent (< 0.5) scaling.

https://doi.org/10.1371/journal.pcbi.1007180.g004

thumbnail
Table 1. Scaling exponents of spatio-temporal gait parameters for five treadmill speeds expressed as the percentage of preferred walking speed (PWS).

PWS is the speed at which a subject choose to walk on treadmill. α(1), α(2), and α(3) are exponents calculated using detrended fluctuation analysis of order 1 (DFA1), 2 (DFA2), 3 (DFA3), and madogram (MD), respectively. Data are presented as mean ± standard deviation. Scaling analysis was performed for the raw experimental time series and time series detrended using the piecewise linear variant of the MARS model. In the latter case the exponents have an L subscript. All indices that are statistically smaller than 0.5 are marked with the dagger symbol.

https://doi.org/10.1371/journal.pcbi.1007180.t001

In the majority of cases, the gait parameters have the following properties:

  • The ST, SL, and SS scaling indices do not depend on the treadmill’s speed.
  • For all three parameters, there are no statistically significant differences between the scaling exponents α(n) and α(MD).
  • The scaling exponents of ST and SL MARS residuals are significantly smaller than those of the corresponding experimental time series.

The exceptions are listed below.

For stride speed, the differences in scaling exponents for different speeds were observed for:

  • DFA1: 100% PWS vs 120% PWS (p = 0.04)
  • DFA2: 80% PWS vs 100% PWS (p = 0.01), 100% PWS vs 120% PWS (p = 0.02)
  • DFA3: 80% PWS vs 100% PWS (p = 0.02)
  • MD: 110% PWS vs 120% PWS (p = 0.01), 100% PWS vs 120% PWS (p = 0.01).

The scaling exponents α(n) and α(MD) were statistically different for:

  • 90% PWS: α(1)α(3) (p = 0.04)
  • 90% PWS: α(1)α(MD) (p = 0.02)
  • 120% PWS: α(1)α(3) (p = 0.002).

For stride speed for 90% PWS: α(3)α(MD) (p = 0.02).

To identify the origin of ST and SL persistence, we determined the piecewise linear MARS trends and the corresponding residuals for each trial. We then created an ensemble of 100 composite signals constructed from the original trends and randomly shuffled residuals. We performed DFAn and madogram analyses on such ensembles. For both gait parameters and all treadmill speeds, α(n) and α(MD) were greater than 0.5. Thus, the outcome of scaling analysis was not determined by the properties of the surrogate noise, which was devoid of any correlations. The results of this numerical experiment for the time series from Fig 1 are displayed in Fig 5.

thumbnail
Fig 5. Boxplots of scaling exponents α(n) and α(MD) for ensembles of signals comprised of MARS trends extracted from experimental time series and randomly shuffled corresponding MARS residuals.

The red disks to the left of each boxplot represent the values of the scaling exponents of the original time series used to generate the ensembles. The calculations were carried out for both stride time (ST) and stride length (SL) time series shown in Fig 1A and 1B. Even though there were no correlations in shuffled MARS residuals, the median values of the scaling exponents of composite signals were significantly greater than 0.5 and similar to the scaling exponents of the original time series. Thus, this numerical experiment shows that DFA detrending is incapable of removing trends from ST and SL time series. The values of the scaling exponents of the experimental time series shown in Fig 1A were: α(1) = 1.02, α(2) = 0.84, α(3) = 0.78, and α(MD) = 0.89 for SL, and α(1) = 1.02, α(2) = 0.91, α(3) = 0.87, and α(MD) = 0.88 for ST. For time series in Fig 1B: α(1) = 0.63, α(2) = 0.67, α(3) = 0.69, and α(MD) = 0.61 for SL, and α(1) = 0.64, α(2) = 0.57, α(3) = 0.56, and α(MD) = 0.51 for ST.

https://doi.org/10.1371/journal.pcbi.1007180.g005

We also calculated α(MD) for phase-randomized surrogates. In particular, for each time series at a given treadmill speed we generated a CCS and an IS. Then, we used either a paired t-test or Wilcoxon signed rank test to compare the scaling index of experimental and surrogate data. We repeated this procedure using n = 50 different realizations of surrogates. The statistically significant differences were sporadic. For SL at the PWS, they were observed in 2% and 10% cases for CCS and IS, respectively. For ST, the corresponding values were equal to 2% and 14%.

Scaling exponent of short fractional Brownian Motion time series

In Fig 6 we show the dependence of the scaling (Hurst) exponent of fractional Brownian motion on the length of the data window. Two ensembles of 500 random walks of length 260 with the Hurst exponent equal to H = 0.40 and H = 0.75 were generated using the Matlab function wfbm. Each trajectory was divided into non-overlapping windows of length k from k = 40 to k = 260 with step 20. For each window, the madogram estimator and detrended fluctuation analysis of order n = 1 to n = 3 were used to compute α(MD) and α(n), respectively. The boxplots of scaling exponents for all four methods are plotted as a function of window length for H = 0.40 (Fig 6A) and H = 0.75 (Fig 6B). These particular values of H were chosen because they do not differ significantly from the average value of scaling exponents of MARS residuals and those of ST and SL time series.

thumbnail
Fig 6. Dependence of scaling exponent of fractional Brownian motion on data window length.

Two ensembles of 500 random walks of length 260 with the Hurst exponent equal to H = 0.40 and H = 0.75 were generated. Each trajectory was divided into non-overlapping windows of length k from k = 40 to k = 260 with step size of 20. For each window, the madogram estimator and detrended fluctuation analysis of order n = 1 to n = 3 were used to compute α(MD) and α(n), respectively. The boxplots of scaling exponents for all four methods are plotted as a function of window length. Subplot (A) shows the results for H = 0.40 and subplot (B) for H = 0.75. For the smallest window with k = 40, the boxplots for α(2) and α(3) were not shown due to the large number of outliers.

https://doi.org/10.1371/journal.pcbi.1007180.g006

We can see in Fig 6 that for short time series DFAn overestimates the Hurst exponent. For example, for DFA1 and k = 40, the median values were equal to α(1) = 0.53 (33% difference) and α(1) = 0.82 (9% difference) for H = 0.40 and H = 0.75, respectively. As window length increases, the estimates of α(n) converge to the true value. For k = 260, the differences dropped to 8% (α(1) = 0.43) and 1% (α(1) = 0.76), for H = 0.40 and H = 0.75, respectively. DFA2 and DFA3 were less accurate. For the smallest window, we did not present the boxplots of scaling exponents for these two methods because of the large number of outliers that would obscure Fig 6. The madogram estimator was the most accurate algorithm. Even for k = 40, the estimates were good: α(MD) = 0.41 and α(MD) = 0.75 for H = 0.40 and H = 0.75, respectively. For H = 0.40 and the largest window (k = 260), the medians of scaling exponents calculated using four methods (madogram and DFAn) were different from each other (p = 1 × 10−42). For H = 0.75 and k = 260, α(MD) was significantly smaller than α(n) for all three orders of detrending. α(1) was smaller than α(3) (p = 7 × 10−12).

Dependence of ST and SL scaling exponents on data window length

Fig 6 serves as a backdrop for the analysis of analogous calculations performed for stride time and stride length time series. We set the length of the largest window to be smaller than the minimal number of strides observed for the trial at a given treadmill’s speed (such lengths varied from 220 to 260). Consequently, the scaling exponents from all subjects could be included in the analysis. For brevity, we confine the presentation of results to the most accurate algorithm—the madogram estimator.

Fig 7 shows the dependence of ST and SL scaling exponents on data window length for the 80% PWS experiment. We can see that for k = 40 the median of ST scaling scaling exponent α(MD) = 0.56 is significantly smaller than α(MD) = 0.76 (p = 4 × 10−5) for the largest window with k = 260. The same effect is observed for stride length: α(MD) = 0.58 and α(MD) = 0.71 for k = 40 and k = 260, respectively (p = 1 × 10−4).

thumbnail
Fig 7. Dependence of scaling exponents of stride time (A) and stride length (B) time series on data window length.

Each time series for 80% PWS experiment was divided into non-overlapping windows of length k from k = 40 to k = 260 with step size of 20. For each window, the madogram estimator and detrended fluctuation analysis of order n = 1 to n = 3 were used to compute α(MD) and α(n), respectively. The boxplots of scaling exponents for all four methods are plotted as a function of window length. For the smallest window with k = 40, the boxplots for α(2) and α(3) were not shown due to the large number of outliers.

https://doi.org/10.1371/journal.pcbi.1007180.g007

In the other four trials, for the smallest window with k = 40 the values of α(MD) for ST were equal to: 0.57, 0.60, 0.62, 0.60 for 90-120% PWS, respectively. For SL, the corresponding values were equal to: 0.57. 0.59, 0.57, 0.59. In all cases, the differences between the smallest and largest window were statistically significant. For k = 40, the median of α(MD) calculated for all five trials was equal to 0.60 and 0.58 for ST and SL, respectively.

Scaling exponents of parts of gait time series with small slope trends

Fig 8 presents the boxplots of scaling exponent α(MD) of the parts of ST and SL time series with small slope trends. More specifically, for each trial we found all segments of MARS trends whose normalized length was greater than 40 and the absolute value of normalized slope was smaller than 0.001 (Fig 3 shows the distribution of ST and SL slopes). Then, we used the madogram estimator to calculate the scaling exponent of the parts of experimental time series with such trends. The boxplots show the exponents aggregated from all trials (80%-120% PWS). There were 267 values for ST (median 0.57) and 249 values for SL (median 0.56).

thumbnail
Fig 8. Boxplots of scaling exponent α(MD) of the parts of stride time (ST) and stride length (SL) series with small slope trends.

For each experimental time series, we found all segments of MARS trends whose normalized length was greater than 40 and the absolute value of normalized slope was smaller than 0.001. Then, we extracted the parts of time series with such trends. We used the madogram estimator to calculate the scaling exponent of selected parts. The boxplots show the exponents aggregated from all trials (80%-120% PWS).

https://doi.org/10.1371/journal.pcbi.1007180.g008

In Table 2 we collected the values of α(MD) for subject 2 for whom in 6 trials the MARS algorithm did not detect linear trends. The median of these values was equal to 0.56.

thumbnail
Table 2. The values of α(MD) for subject 2’s trials with no MARS trends.

The scaling exponents are shown for stride length (SL) and stride time (ST).

https://doi.org/10.1371/journal.pcbi.1007180.t002

Correlation coefficients

Fig 9 shows the boxplots of Pearson’s correlation coefficient ρ between stride time and stride length for all five treadmill speeds. Correlations were calculated for: experimental (raw) time series, trends determined using the MARS algorithm, and MARS residuals (noise). ρ was independent of speed: praw = 0.12, ptrend = 0.35, and pnoise = 0.17. For the PWS, ρraw = 0.50 ± 0.17, ρtrend = 0.79 ± 0.22, and ρnoise = 0.21±0.10. Thus, coupling between trends was determined to be very strong. The same holds true for the cross-correlated phase-randomized surrogates. At the PWS, was not statistically different from ρtrend. One can see in Fig 10A that the correlation coefficient between the piecewise linear MARS trends in the independently phase-randomized SL and ST surrogates is close to zero.

thumbnail
Fig 9. Pearson’s correlation coefficient between stride time and stride length at five different treadmill speeds.

Correlations were calculated for: experimental (raw) time series, trends determined using the MARS algorithm, and MARS residuals (noise).

https://doi.org/10.1371/journal.pcbi.1007180.g009

thumbnail
Fig 10.

(A) Correlation coefficient between the piecewise linear MARS trends in stride length and stride time. (B) Coefficient of variation of trend speed for the same data as in the top subplot. The boxplots are shown for the original time series, phase-randomized surrogates, and cross-correlated (SL-ST) phase-randomized surrogates.

https://doi.org/10.1371/journal.pcbi.1007180.g010

Trend speed

The examples of time evolution of trend speed v(trend) are shown in Fig 1. The boxplots of the trend speed control parameter TSC in Fig 11 show that v(trend) is tightly controlled about the treadmill speed. For example, at the preferred walking speed, TSC = 0.13. For the same speed, the coefficient of variation of v(trend) was equal to 0.6% (see Fig 12) and was about three times smaller than that of ST (1.6%), SL (1.8%), and SS (1.7%). For all treadmill speeds, the trend speed COV for the cross-correlated phase-randomized surrogates was not different from that of the experimental data and much smaller that of independently phase-randomized surrogates (Fig 10B).

thumbnail
Fig 11. Boxplots of trend speed control parameter (Eq (6)) for five treadmill speeds.

The results are shown for the original time series, phase-randomized surrogates, and cross-correlated (SL-ST) phase-randomized surrogates.

https://doi.org/10.1371/journal.pcbi.1007180.g011

thumbnail
Fig 12. Boxplots of coefficient of variation for trend speed, stride speed, stride length, and stride time.

https://doi.org/10.1371/journal.pcbi.1007180.g012

Discussion

Detrended fluctuation analysis is one of the most frequently used algorithms for fractal analysis of experimental time series. At the time of writing, the seminal paper of Peng et al. [4] has been cited almost 1800 times. It is difficult to overestimate the importance of this method in general and in quantitative gait analysis in particular [1, 2, 14, 1627].

Early on, numerical investigations showed that simple non-stationarities, such as sinusoidal periodicity or monotonic global trends [28, 29], as well as more complex sources of non-stationarity [30], affected DFA estimates of the scaling exponent. Bryce and Sprague correctly pointed out that, paradoxically, rather than raising doubts about the purported universal applicability of DFA for trended processes, those papers were actually used to demonstrate the efficacy of the method [31].

In DFA, detrending is performed on the integrated signal. Consequently, global linear trends can be eliminated using second or higher order DFA [28, 32]. Fig 1A shows that prominent piecewise linear trends in stride time and stride length may appear even during treadmill walking. Thus, it is surprising that DFA1, which cannot mitigate even global linear trends, was used in most of the previous studies of scaling in spatio-temporal gait parameters. There are two possible explanations for the prevalence of DFA1 in gait analysis. First, while global linear trends distinctly manifest themselves as crossovers in DFA fluctuation function, no apparent crossovers are observed in gait data. Second, the scaling exponents determined using DFA1 and DFA2 are not statistically different, as is evident from the data collected in Table 1, as well as from the visualizations of their distributions in Fig 4. This agreement was presumably used in favor of DFA1.

We demonstrated that the order of DFA detrending did not have a statistically significant effect on the values of the scaling exponents of the experimental time series. Moreover, α(0), α(1), and α(2) were not significantly different from α(MD). Thus, we infer that DFA detrending is ineffective, and the observed strong persistence of ST and SL can not stem from the properties of residuals. This interpretation is supported by DFA of signals made up of the MARS trends of experimental time series and randomly shuffled corresponding MARS residuals (shuffled residuals are obviously devoid of any correlations). The fact that α(n) for both ST and SL were sharply greater than 0.5 (see Fig 5A) showed that, as in the case of experimental data, strong persistence was associated with the presence of trends.

Let us begin the discussion of significance of trends in gait dynamics with a summary of scaling properties of gait parameters. We focus on the results obtained using the most accurate algorithm—the madogram estimator:

  • ST and SL scaling exponents α(MD) increase with the data window length (Fig 7). For the shortest analyzed window (k = 40), the medians of all trials are respectively equal to 0.60 and 0.58. For k = 260, fluctuations of gait parameters are strongly persistent (0.75 and 0.71).
  • The median of scaling exponent of the parts of the experimental time series with small MARS slopes are equal to 0.57 and and 0.56 for ST and SL, respectively (Fig 8, see also Table 2).
  • For both the ST and SL MARS residuals, α(MD) averaged over the treadmill speed is equal to 0.48 (Table 1).

The first two properties strongly indicate that the values of scaling indices of stride time and length are determined by the superposition of large scale trends and small scale fluctuations (note that the fractal dimension of a piecewise linear curve is 1).

We use the second and third property as evidence in support of the hypothesis that trends serve as control manifolds about which ST and SL fluctuate.

Before we discuss the second hypothesis in detail, let us note that during treadmill walking a subject’s stride time and stride length must yield a stride speed that fluctuates over a narrow range centered on the treadmill belt’s speed. In principle, there are infinitely many such combinations which satisfy this type of constraint. These combinations form a straight line in the ST-SL phase space—a “Goal Equivalent Manifold” (GEM) [15]. Dingwell et al. projected the deviation vector (calculated in ST-SL space with respect to the mean values of ST and SL) onto the GEM and the axis perpendicular to it. It turned out that fluctuations of tangential and transverse components were persistent and anti-persistent, respectively. Moreover, the tangential variability was higher than the transverse one. Thus, these statistical properties provided the evidence that subjects did not regulate ST and SL independently but instead adjusted them in a coordinated manner to maintain walking speed.

Herein we not only provide direct proof of such interdependence but also elucidate how speed control during treadmill walking might actually be accomplished. The distribution of the trend speed control parameter (TSC) in Fig 11 and coefficient of variation of v(trend) in Fig 12 show that the trend speed v(trend), defined as the ratio of instantaneous values of piecewise linear SL and ST MARS trends, is tightly controlled about the treadmill speed. The strong persistence of stride time and stride length does not impair a subject’s ability to maintain speed because of the strong coupling between their trends (e.g. ρtrend = 0.85 at PWS). The concomitant changes of instantaneous values of ST and SL trends correspond to movement along the GEM as postulated by Dingwell et al. in their model of redundant stride speed control [33, 34].

We also calculated the correlation coefficient between ST and SL trends as well as the trend speed control parameter using the data from the recent research of Roerdink et al. [35]. While the strength of correlation was comparable to the values presented in this paper, the TSC parameter was smaller, indicating even tighter control of v(trend). The detailed analysis of Roerdink’s data will be presented elsewhere.

Humans have an innate ability to synchronize their movements with rhythmic sound stimuli. Sejdic et al. argue that sound associated with bipedal gait influenced the evolution of human auditory-motor rhythmic abilities [36]. Walking in sync with a metronome has previously been investigated as a potential rehabilitation tool in patients with Parkinson’s disease [37] and stroke [38]. Isochronous metronome induces anti-persistence in stride time which raises doubts about the efficacy of such rehabilitation. Application of persistent metronomes to gait rehabilitation is an active field of research [2023, 39].

In this work we demonstrated that the persistence of stride time and length originates from trends that may be approximated by the piecewise linear MARS model. We are not aware of any previous studies that explicitly analyzed the properties of trends in human gait. The normalized trend durations of ST and SL have a distribution with an exponential tail. The exponential distribution describes the time between events in a Poisson point process, i.e. a process in which events occur continuously and independently at a constant average rate. Further research is needed to verify whether this kind of process is involved in the generation of gait trends. The scaling indices and the distributions of trend properties (durations and slopes) of cross-correlated phase-randomized SL-ST surrogates are either the same or very similar to those of the experimental data. This is a plausible indication that the speed control during treadmill walking is predominately linear as suggested by [15]. We believe that better understanding of gait control mechanisms will pave the way for novel, more effective strategies of gait rehabilitation.

Conclusions

Hausdorff et al. [1, 2] discovered long-range, persistent correlations in stride duration of human gait nearly a quarter century ago. This persistence was attributed to correlated noise superposed on trends. The evidence presented herein necessitates a fundamental revision of such an interpretation: the statistical properties of ST and SL time series stem from the superposition of large scale trends and small scale fluctuations. We demonstrated that the ST and SL trends serve as control manifolds and that the trend speed is tightly controlled about the treadmill speed. The strong coupling between the ST and SL trends ensures that the concomitant changes of their instantaneous values correspond to movement along the constant speed GEM as postulated by Dingwell et al. [33, 34].

Materials and methods

Experimental data

In our analysis, we used data from the study of Dingwell et al. [15], which were available in the Dryad repository [40]. Seventeen young healthy adults were asked to walk on a motor-driven treadmill for five minutes, at five different speeds, determined as a percentage of a subject’s preferred walking speed (PWS). PWS is the speed at which a subject choose to walk on treadmill. There were two trials for each speed (80%, 90%, 100%, 110%, 120% PWS). SL, ST, and SS were determined using a motion capture system. A comprehensive description of the experimental protocol and participants’ characteristics can be found in the original paper.

Multivariate Adaptive Regression Splines (MARS)

To model trends in spatio-temporal gait parameters we adopted MARS (Multivariate Adaptive Regression Splines)—a nonparametric adaptive regression method proposed by J. Friedman [41, 42]. Let be a time series of experimental values observed at instances . The MARS model f is a linear combination of basis functions hm: (1) We denote by the set of basis functions hm which are constructed in a very specific way. Let us consider two piecewise linear functions max(0,x−t) and max(0,t-x) with a knot at t. These two functions (linear splines) form a so-called reflected pair. The example of such a pair is given in Fig 13.

thumbnail
Fig 13. A pair of reflected basis functions with a knot at t = 5.

https://doi.org/10.1371/journal.pcbi.1007180.g013

We begin the construction of by creating a set C of N reflected pairs with knots equal to values xi: (2) The MARS model is build iteratively. A new basis function hM+1 is the product of one of the already constructed basis functions hl and one of the reflected pairs from set C: (3)

From all possible (M + 1)N products of this form, we select the one that gives the maximum reduction in sum-of-squares residual error. Here and are coefficients estimated by least squares, along with all the other M + 1 coefficients in the model. The first basis function h0(x) is constant equal to 1. The process of adding new terms, which is often called a forward model-building procedure (forward pass), is continued until the change in residual error is smaller than the predefined stopping condition or until the maximum number of terms is reached.

The model built during the forward pass is usually overfitted and needs to be pruned to obtain better generalization ability. During each stage of a backward deletion procedure (backward pass), the term whose removal causes the smallest increase in residual squared error is deleted. In this way, we find the best model with λ terms. The goal of the pruning process is to minimize the value of generalized cross-validation: (4) in order to estimate the optimal number of terms. M(λ) is the effective number of parameters in the model: (5) that accounts for both the number of linearly independent basis functions r and the number of knots K in the forward pass. c is a penalty for adding a new knot.

In this study, we used ARESLab, an open source MATLAB implementation of MARS [43], to approximate trends in gait spatiotemporal parameters using a piecewise linear, additive MARS model. In other words, every new basis function was just a linear combination of functions which belonged to a reflected pair (maximum order of interaction was equal to 1). The stopping condition for the forward phase was set to 0.001, and the generalized cross-validation knot penalty c was set equal to 2, the value most commonly used for additive models [42].

Trend speed

We define a trend speed as the ratio of values of piecewise linear SL and ST MARS trends at ith stride. To quantify deviations of the trend speed from the average stride speed < SS > during a given trial we introduce a trend speed control (TSC) parameter: (6) Thus, TSC is equal to zero when the trend and mean stride speed are equal to each other during the whole trial.

Scaling analysis

Detrended fluctuation analysis (DFA).

DFA was designed to detect long-time, power-law scaling of the second moment in the presence of additive, polynomial non-stationarity. Given a bounded time series , one creates an unbounded process: (7) where < y > denotes the mean value of the time series. Next, Yt is divided into windows of length n and a local least-squares polynomial fit is performed within each window. Let Pt indicate the piece-wise sequence of such polynomial fits. Then, the root-mean-square deviation from the trend, the fluctuation function, is calculated: (8) A straight line on a log-log graph of F(n) as a function of n is considered a manifestation of self-affinity as F(n) ∝ nα. For fractional Gaussian noise, the exponent (index) α varies between 0 and 1. In this case, for α < 0.5 fluctuations are anti-persistent and for α > 0.5 they are persistent. The n-th order polynomial regression in the DFA family is typically denoted as DFAn. To facilitate reproducibility of the results, we used the MATLAB function dfa from the WFDB Toolbox [44, 45] for DFA calculations.

Madogram estimator (MD).

The variogram of order p of a stochastic process with stationary increments is defined as one half times the expectation value of an increment at lag t [46] (9) For p = 1 we call such a structure function the madogram. As t → 0 (10) where H ∈ (0, 1] and constants β and cP are positive. For one dimensional time series , we may define the power variation of order p: (11) and using Eq (10) derive the following estimate of fractal dimension: (12) The fractal dimension and scaling exponent α are related by the following equation: (13) Despite its simplicity, the madogram estimator (p = 1) turns out to be particularly robust, especially for non-Gaussian processes. Gneiting er al. [46] compare the properties of fractal dimension estimators that were implemented in the R package fractaldim.

Surrogate data

We used the program surrogates from TISEAN Nonlinear Time Series Analysis library [47] to create SL/ST/SS independently phase-randomized surrogate time series (IS) as well as cross-correlated (SL-ST) phase-randomized surrogates (CCS).

Statistical analysis

We used the Shapiro-Wilk test to determine whether the analyzed data were normally distributed. The significance threshold for all the statistical tests was set to 0.05.

The dependence of trend durations of stride time and stride length on treadmill speed were examined with the Kruskal-Wallis test with Tukey’s post hoc comparisons.

We compared the values of Pearson’s correlation coefficient between stride time and stride length (experimental values, piecewise linear trends, and MARS residuals) at different speeds using ANOVA or the Kruskal-Wallis test (depending on normality of data). In both cases we used Tukey’s post hoc comparisons.

For a given gait parameter and treadmill speed, the differences in scaling exponents determined using DFAn and madogram estimator were assessed with either ANOVA or the Kruskal-Wallis test (with Tukey’s post hoc comparisons in both cases). Such statistical analysis was performed for both original and MARS detrended time series.

To determine whether fluctuations of spatio-temporal gait parameters were anti-persistent, we checked whether the corresponding scaling exponents were smaller than 0.5 using left-sided tests (t-test or the Wilcoxon signed rank test).

The Kolmogorov-Smirnov test was used to compare distributions and to verify whether a sample was drawn from a given distribution. In a few instances, we also used the Anderson-Darling test to show that the drawn conclusions might depend on the choice of a statistical test.

In the presented boxplots, the black and grey dots correspond to the outliers and far outliers, respectively.

The analysis of statistical properties of stride time and length trends was carried out using Mathematica 11.3 (Wolfram Research).

Software

GaitTrends repository [48] contains the MATLAB scripts that we developed while working on this paper. The code which is distributed under the GNU GPL license was tested using MATLAB R2016b and R2018a (Mathworks) on Windows 10. The included user’s manual shows how to reproduce all the results presented in this paper.

References

  1. 1. Hausdorff JM, Peng CK, Ladin Z, Wei JY, Goldberger AL. Is walking a random walk? Evidence for long-range correlations in stride interval of human gait. Journal of Applied Physiology. 1995;78(1):349–358. pmid:7713836
  2. 2. Hausdorff JM, Purdon PL, Peng CK, Ladin Z, Wei JY, Goldberger AL. Fractal dynamics of human gait: stability of long-range correlations in stride interval fluctuations. Journal of Applied Physiology. 1996;80(5):1448–1457. pmid:8727526
  3. 3. Feder J. Fractals. Springer US; 1988. Available from: https://doi.org/10.1007%2F978-1-4899-2124-6.
  4. 4. Peng CK, Buldyrev SV, Havlin S, Simons M, Stanley HE, Goldberger AL. Mosaic organization of DNA nucleotides. Physical Review E. 1994;49(2):1685–1689. pmid:9961383
  5. 5. Goldberger AL, Rigney DR, West BJ. Science in Pictures: Chaos and Fractals in Human Physiology. Scientific American. 1990;262(2):42–49. pmid:2296715
  6. 6. Bassingthwaighte JB, Liebovitch LS, West BJ. Fractal Physiology. Springer New York; 1994. Available from: https://doi.org/10.1007%2F978-1-4614-7572-9.
  7. 7. West BJ. Fractal physiology and the fractional calculus: a perspective. Frontiers in Physiology. 2010;1. pmid:21423355
  8. 8. Peng CK, Buldyrev SV, Hausdorff JM, Havlin S, Mietus JE, Simons M, et al. Non-equilibrium dynamics as an indispensable characteristic of a healthy biological system. Integrative Physiological and Behavioral Science. 1994;29(3):283–293. pmid:7811648
  9. 9. Hausdorff J, Mitchell S, Firtion R, Peng C, Cudkowicz M, Wei J, et al. Altered fractal dynamics of gait: reduced stride-interval correlations with aging and Huntington’s disease. J Appl Physiol (1985). 1997;82:262–9. pmid:9029225
  10. 10. Goldberger A, Amaral L, Hausdorff J, Ivanov P, Peng C, Stanley H. Fractal dynamics in physiology: alterations with disease and aging. Proc Natl Acad Sci U S A. 2002;99 Suppl 1:2466–72. pmid:11875196
  11. 11. Delignières D, Torre K. Fractal dynamics of human gait: a reassessment of the 1996 data of Hausdorff et al. Journal of Applied Physiology. 2009;106(4):1272–1279. pmid:19228991
  12. 12. Terrier P, Dériaz O. Persistent and anti-persistent pattern in stride-to-stride variability of treadmill walking: Influence of rhythmic auditory cueing. Human Movement Science. 2012;31(6):1585–1597. pmid:23164626
  13. 13. Terrier P. Fractal Fluctuations in Human Walking: Comparison Between Auditory and Visually Guided Stepping. Annals of Biomedical Engineering. 2016;44(9):2785–2793. pmid:26903091
  14. 14. Roerdink M, Daffertshofer A, Marmelat V, Beek PJ. How to Sync to the Beat of a Persistent Fractal Metronome without Falling Off the Treadmill? PLOS ONE. 2015;10(7):e0134148. pmid:26230254
  15. 15. Dingwell JB, Cusumano JP. Re-interpreting detrended fluctuation analyses of stride-to-stride variability in human walking. Gait & Posture. 2010;32(3):348–353. pmid:20605097
  16. 16. Terrier P, Turner V, Schutz Y. GPS analysis of human locomotion: Further evidence for long-range correlations in stride-to-stride fluctuations of gait parameters. Human Movement Science. 2005;24(1):97–115. pmid:15896861
  17. 17. Jordan K, Challis JH, Newell KM. Walking speed influences on gait cycle variability. Gait & Posture. 2007;26(1):128–134. pmid:16982195
  18. 18. Hausdorff JM. Gait dynamics fractals and falls: Finding meaning in the stride-to-stride fluctuations of human walking. Human Movement Science. 2007;26(4):555–589. pmid:17618701
  19. 19. Lamoth CJ, van Deudekom FJ, van Campen JP, Appels BA, de Vries OJ, Pijnappels M. Gait stability and variability measures show effects of impaired cognition and dual tasking in frail people. Journal of NeuroEngineering and Rehabilitation. 2011;8(1). pmid:21241487
  20. 20. Kaipust JP, McGrath D, Mukherjee M, Stergiou N. Gait Variability is Altered in Older Adults When Listening to Auditory Stimuli with Differing Temporal Structures. Annals of Biomedical Engineering. 2013;41(8):1595–1603. pmid:22956164
  21. 21. Rhea C, Kiefer A, Wittstein M, Leonard K, MacPherson R, Wright W, et al. Fractal gait patterns are retained after entrainment to a fractal stimulus. PLoS One. 2014;9:e106755. pmid:25221981
  22. 22. Rhea CK, Kiefer AW, D’Andrea SE, Warren WH, Aaron RK. Entrainment to a real time fractal visual stimulus modulates fractal gait dynamics. Human Movement Science. 2014;36:20–34. pmid:24911782
  23. 23. Marmelat V, Torre K, Beek PJ, Daffertshofer A. Persistent Fluctuations in Stride Intervals under Fractal Auditory Stimulation. PLoS ONE. 2014;9(3):e91949. pmid:24651455
  24. 24. Moon Y, Sung J, An R, Hernandez ME, Sosnoff JJ. Gait variability in people with neurological disorders: A systematic review and meta-analysis. Human Movement Science. 2016;47:197–208. pmid:27023045
  25. 25. Kuznetsov N, Rhea C. Power considerations for the application of detrended fluctuation analysis in gait variability studies. PLoS One. 2017;12:e0174144. pmid:28323871
  26. 26. Choi J, Kang D, Seo J, Tack G. Fractal fluctuations in spatiotemporal variables when walking on a self-paced treadmill. J Biomech. 2017;65:154–160. pmid:29096982
  27. 27. Ducharme S, Liddy J, Haddad J, Busa M, Claxton L, van ER. Association between stride time fractality and gait adaptability during unperturbed and asymmetric walking. Hum Mov Sci. 2018;58:248–259. pmid:29505917
  28. 28. Kantelhardt JW, Koscielny-Bunde E, Rego HHA, Havlin S, Bunde A. Detecting long-range correlations with detrended fluctuation analysis. Physica A: Statistical Mechanics and its Applications. 2001;295(3-4):441–454.
  29. 29. Hu K, Ivanov PC, Chen Z, Carpena P, Stanley HE. Effect of trends on detrended fluctuation analysis. Physical Review E. 2001;64(1). pmid:11461232
  30. 30. Chen Z, Ivanov PC, Hu K, Stanley HE. Effect of nonstationarities on detrended fluctuation analysis. Physical Review E. 2002;65(4). pmid:12005806
  31. 31. Bryce RM, Sprague KB. Revisiting detrended fluctuation analysis. Scientific Reports. 2012;2(1). pmid:22419991
  32. 32. Bashan A, Bartsch R, Kantelhardt JW, Havlin S. Comparison of detrending methods for fluctuation analysis. Physica A: Statistical Mechanics and its Applications. 2008;387(21):5080–5090.
  33. 33. Dingwell JB, John J, Cusumano JP. Do Humans Optimally Exploit Redundancy to Control Step Variability in Walking? PLoS Computational Biology. 2010;6(7):e1000856. pmid:20657664
  34. 34. Dingwell JB, Cusumano JP. Identifying stride-to-stride control strategies in human treadmill walking. PloS one. 2015;10(4):e0124879. pmid:25910253
  35. 35. Roerdink M, De Jonge CP, Smid LM, Daffertshofer A. Tightening up the control of treadmill walking: Effects of maneuverability walking: Effects of maneuverability range and acoustic pacing on stride-to-stride fluctuations. Frontiers in Physiology. 2019;10(MAR):257. pmid:30967787
  36. 36. Sejdić E, Fu Y, Pak A, Fairley JA, Chau T. The effects of rhythmic sensory cues on the temporal dynamics of human gait. PLoS ONE. 2012;7(8). pmid:22927946
  37. 37. Suteerawattananon M, Morris G, Etnyre B, Jankovic J, Protas E. Effects of visual and auditory cues on gait in individuals with Parkinson’s disease. Journal of the neurological sciences. 2004;219(1-2):63–69. pmid:15050439
  38. 38. Roerdink M, Lamoth CJ, van Kordelaar J, Elich P, Konijnenbelt M, Kwakkel G, et al. Rhythm perturbations in acoustically paced treadmill walking after stroke. Neurorehabilitation and Neural Repair. 2009;23(7):668–678. pmid:19307435
  39. 39. Hunt N, McGrath D, Stergiou N. The influence of auditory-motor coupling on fractal dynamics in human gait. Scientific Reports. 2014;4(1). pmid:25080936
  40. 40. Dingwell JB, Cusumano JP. Data from: Do humans optimally exploit redundancy to control step variability in walking? Dryad Digital Repository.; 2015. Available from: http://dx.doi.org/10.5061/dryad.sk55m.
  41. 41. Friedman JH. Multivariate Adaptive Regression Splines. The Annals of Statistics. 1991;19(1):1–67.
  42. 42. Hastie T, Tibshirani R, Friedman J. The Elements of Statistical Learning. Springer New York; 2009. Available from: https://doi.org/10.1007%2F978-0-387-84858-7.
  43. 43. Jēkabsons G. ARESLab: Adaptive Regression Splines toolbox; 2016. Available from: http://www.cs.rtu.lv/jekabsons/regression.html.
  44. 44. Moody GB, Mark RG, Goldberger AL. PhysioNet: a research resource for studies of complex physiologic and biomedical signals. In: Computers in Cardiology 2000. Vol.27 (Cat. 00CH37163). IEEE;.Available from: https://doi.org/10.1109%2Fcic.2000.898485.
  45. 45. Silva I, Moody GB. An Open-source Toolbox for Analysing and Processing PhysioNet Databases in MATLAB and Octave. Journal of Open Research Software. 2014;2. pmid:26525081
  46. 46. Gneiting T, Ševčíková H, Percival DB. Estimators of fractal dimension: Assessing the roughness of time series and spatial data. Statistical Science. 2012; p. 247–277.
  47. 47. Hegger R, Kantz H, Schreiber T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos: An Interdisciplinary Journal of Nonlinear Science. 1999;9(2):413–435. pmid:12779839
  48. 48. Kozlowska K, Latka M. GaitTrends. Software from: Significance of trends in gait dynamics.; 2020. Available from: https://github.com/MiroslawLatka/GaitTrends.git.