Quantifying workload using nonlinear dynamical measures of biomechanical parameters during cycling on a roller trainer

Ann-Kathrin Harsch; Alexander Kunert; Daniel Koska; Christian Maiwald

doi:10.1371/journal.pone.0285408

Abstract

The aim of the present study was to determine the effectiveness of nonlinear parameters in distinguishing individual workload in cycling by using bike-integrated sensor data. The investigation focused on two nonlinear parameters: The ML1, which analyzes the geometric median in phase space, and the maximum Lyapunov exponent as nonlinear measure of local system stability. We investigated two hypothesis: 1. ML1_α, derived from kinematic crank data, is as good as ML1_F, derived from force crank data, at distinguishing between individual load levels. 2. Increasing load during cycling leads to decreasing local system stability evidenced by linearly increasing maximal Lyapunov exponents generated from kinematic data.

A maximal incremental cycling step test was conducted on an ergometer, generating complete datasets from 10 participants in a laboratory setting. Pedaling torque and kinematic data of the crank were recorded. ML1_F, ML1_α, and Lyapunov parameters (λ_st, λ_lt, ι_st, ι_lt) were calculated for each participant at comparable load levels. The results showed a significant linear increase in ML1_α across three individual load levels, with a lower but still large effect compared to ML1_F. The contrast analysis also confirmed a linearly increasing trend for λ_st across three load levels, but this was not confirmed for λ_lt. However, the intercepts ι_st and ι_lt of the short- and longterm divergence showed a statistically significant linear increase across the load levels.

In summary, nonlinear parameters seem fundamentally suitable to distinguish individual load levels in cycling. It is concluded that higher load during cycling is associated with decreasing local system stability. These findings may aid in developing improved e-bike propulsion algorithms. Further research is needed to determine the impact of factors occurring in field application.

Citation: Harsch A-K, Kunert A, Koska D, Maiwald C (2023) Quantifying workload using nonlinear dynamical measures of biomechanical parameters during cycling on a roller trainer. PLoS ONE 18(5): e0285408. https://doi.org/10.1371/journal.pone.0285408

Editor: Angelo Marcelo Tusset, Federal University of Technology - Parana, BRAZIL

Received: August 26, 2022; Accepted: April 21, 2023; Published: May 9, 2023

Copyright: © 2023 Harsch et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant datasets are now available at the Zenodo database (DOI: 10.5281/zenodo.7338904).

Funding: The study is part of a research project which was grant-aided by the German Federal Ministry of Economic Affairs and Climate Action (ID: VF150010). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

When using e-bikes, the control and regulation of propulsion is important for ensuring a natural riding experience and for optimizing range. However, users can often only specify a level of support relative to their own pedaling power by selecting a predetermined level. Thus, the user’s individual effort is not factored into the algorithms generating output levels. To develop more sophisticated algorithms based on individual effort, parameters are needed that allow the adequate characterization of individual load properties. Preferably, these should be based on information from sensors that are mechanically integrated into the bike, rather than additional physiological devices like heart rate monitors.

Biomechanical studies in cycling often analyze pedaling movement behavior to optimize force transmission or training techniques using linear and nonlinear methods. Nonlinear methods are applied to examine the relationship between external load and movement behavior by analyzing changes in cadence, pedal torque on incline, or higher power output [1–5]. Nonlinear dynamic methods can offer an added value in this context, since they consider the spatio-temporal evolution of a system. Compared to linear methods, fluctuations do not necessarily represent noise, but may indicate underlying properties of the investigated system [6].

A basis for the application of nonlinear methods in pedaling was provided by Quintana-Duque and Saupe [7], who showed that knee motion data in indoor pedaling underlies deterministic chaos and that a chaotic dynamical system can be assumed. More advanced nonlinear methods, for example, can distinguish different body positions via time variability (changing behavior in cycle duration) [8]. Furthermore, non-physiological parameters, such as crank torque variability [9], seem to be effected by pedaling intensity.

In previous studies [10, 11], we analyzed crank forces and developed a nonlinear parameter that correlates with the individual load when, for example, changes in terrain inclination are prevalent. This parameter, called ML1, is defined by the vector length of the geometric median of the force points in two dimensional phase space. It is also associated with a physiological parameter (heart rate), but shows less delay in response to terrain-induced load changes. This was already demonstrated in laboratory and field tests. Therefore, ML1 appears to be suitable for being integrated into control systems of e-bike propulsion. However, attaching a force sensor to the crank is expensive. A more affordable alternative could be to integrate inertial sensors, such as gyroscopes, to the crank arm. Therefore, the first question of the present work is to investigate the extent to which the geometric median calculated via kinematic signals (ML1_α) agrees with the geometric median from force data (ML1_F) and can also distinguish individual load levels.

Variability, stability, and/or complexity of a movement are examined within the context of nonlinear dynamic system analysis in biomechanics [12]. These terms are related, but should not be interchanged. In the present paper, we follow van Emmerik et al. [6], from the field of gait analysis, and define the local stability of a movement system as the resistance to small perturbations. In gait analyses, local system stability is often quantified using the maximum Lyapunov exponent λ_max [13–15]. λ_max describes the rate of divergence of originally proximate points on trajectories in multidimensional phase space. Thus, λ_max characterizes the system’s sensitivity to local perturbations [16]. The question remains whether the Lyapunov approach can be used to quantify load-dependent changes in pedaling behavior when kinematic data from accelerometers on the crank are considered (research question II).

The overall aim of this paper is to determine the extent to which nonlinear methods in combination with kinematic data, are suited to describe individual load in cycling. In addition to the first research question, regarding the extent of agreement between ML1_α and ML1_F (1), we investigate how the local dynamic stability as quantified by the maximum Lyapunov exponent is able to distinguish between individual load levels (2).

Materials and methods

Participants

15 healthy participants were informed verbally and in writing about the content of the study and then gave written approval to participate. It was approved by the Ethics Committee at Chemnitz University of Technology (reference number V-337–17-CM-Radsport-17052019).

Previous studies have investigated heterogeneous groups [10, 11]. Therefore, we did not focus on typical e-bike users to generate a more homogenous group regarding age and ability. To ensure good endurance capacity of the participants, only active male cyclists and triathletes with at least 5 h/week training volume were included in the study. Training volume and anthropometric data for the cyclists in this study are shown in Table 1.

Download:

Table 1. Training volume and anthropometric data for the cyclists (all male).

https://doi.org/10.1371/journal.pone.0285408.t001

Experimental setup

The participants performed a maximal incremental cycling test. Power output started at 100 W and was increased incrementally by 20 W each 3 minutes until the cyclists stopped the test due to exhaustion or permanently fell below a cadence of 60 cycles per minute. Therefore, each participant finished an individual number of intervals. No further restrictions regarding cadence (for example by using a metronome) were implemented, since Warlop et al. (2013) indicate that this can modify nonlinear properties of the pedaling behavior [12]. All tests were performed on the same bicycle (Stein Bikes Mauna Loa 29AL, Stein Bikes Zweiradhandel GmbH, Chemnitz) mounted on a roller trainer (Cyclus 2, RBM elektronik-automation GmbH). The power output was controlled using the Cyclus 2 software. All participants used the same non-cleat platform pedals to ensure comparable force transmission. Before starting the test, participants determined their individual saddle height. This could be adjusted after a 10 min warm up period, which was standardized to 100 W.

The raw pedaling torque was measured using a Stages Shimano XT Powermeter (Stages Inc., USA) with a crank length of 175 mm and a sampling rate of 64 Hz. Heart rate was recorded using a standard chest strap (Berry King, bestbeans UG, Germany). Raw data were transferred to a Windows PC using bluetooth protocol and were recorded in LabView (National Instruments, USA). Crank kinematics were captured using a gyroscope with an internal data logger (Envisible Dialogg Steinbeis-Forschungszentrum Human Centered Engineering, Chemnitz, DEU), which was attached to the left crank arm and recorded angular velocities at a sampling rate of 500 Hz. To control for fatigue, respiratory gases were analyzed with a spirometer (Oxycon Pro, Erich Jaeger, Germany). Exclusion criterion was a respiratory quotient greater than 1.

Data analysis

All data analyses were performed in R (version 4.2), using R Studio [17]. The analyses included crank angular velocity and force data.

Ten complete datasets were generated. Three participants had to be excluded due to technical difficulties and missing recorded data, caused by battery failure on the crank or chest strap. Two other trials were excluded because of implausible gyroscope data, and because of an implausible increase of one participant’s heart rate to 230 beats/min.

Raw velocity data were first filtered using a fourth-order low-pass Butterworth filter with a cut-off frequency of 120 Hz. To facilitate the calculation of Lyapunov parameters, time normalization was performed for the entire time series of each of the 3 minute power intervals. This involved counting the number of cycles (revolutions of the crank) in the current interval, which are defined as n_cyc, and normalizing the entire 3-minute time series to the number of cycles multiplied by 100. Thus, the number of points per power interval I_p was: (1)

Hence, information about the relative duration of one cycle was preserved.

The number of fulfilled power intervals was different for each participant, depending on the participant’s maximum achieved power. The relative power interval I_p,rel was determined to ensure comparability between individually experienced loads. I_p,rel includes the interval at 25%, 50% and 75% of the power difference of the start and the highest fully completed interval (I_p,max). This power difference is defined as ΔP = I_p,max − I_p,start. For example, the 75% relative power interval I_p,75 is then calculated as: (2)

For percentages between two intervals, the next higher level was set as I_p,rel.

Assumptions and testing for nonlinearity.

One assumption of the nonlinear methods applied in this study is the stationarity [16]. To ensure stationarity of the angular velocity signal, an acceleration measure was defined using pointwise differentiation of the gyroscope’s angular velocity signal Δω_i: (3)

Another prerequisite for the use of nonlinear analysis is that the data originates from a nonlinear process. Therefore, we first determined whether the data were subject to linear stochastic noise or nonlinearity by using the surrogate data method suggested in [16, 18]. This method was implemented in the R-package nonlinearTseries [19]. The function surrogateTest() was applied to the whole time series of the filtered and normalized differentiated angular velocity data for each participant.

Within this function, the null hypothesis is tested, that the data originates from a Gaussian linear stochastic process. Surrogate data are generated by using a phase randomization procedure and values of a discriminating statistic for the original (input) data and the surrogate data are compared. The whole time series of the filtered and normalized acceleration data for each participant was used as original (input) data for the test. Graphical results of the discriminating statistic and the original series can be found in the (see: S2 Fig). The null hypothesis was rejected for all participants, therefore the prerequisites for applying nonlinear methods are fulfilled for the tested data.

Nonlinear parameter ML1.

The novel ML1 method was applied to the pedaling forces using a sliding window algorithm [10, 11]. The window width was set to 50 normalized cycles (interpolation to 100 samples per cycle), representing 5000 points in phase space (see Fig 1). This resulted in a series of ML1 values in which each value represents the nonlinear characteristic of the last 50 force-time cycles. For each I_p,rel, the mean of all ML1 values was calculated and used for further statistical analyses.

Download:

Fig 1. Method for calculating ML1 using the example ML1_F.

A sliding window considering 50 full z-standardized and time-normalized pedaling cycles to calculate ML1_F (see bottom of Fig 1; 30 instead of 50 pedaling-cycles are considered). The geometric median in 2-dimensional phase-space with a 25% delay was determined (see right side of Fig 1) for these cycles.

https://doi.org/10.1371/journal.pone.0285408.g001

As input data for the ML1, the differentiated acceleration data (see Eq 3) of each I_p,rel was used to obtain the ML1_α. To determine ML1_F, force values were determined for the pedaling torque and time-normalized to 100 points per cycle. Furthermore, the force values were z-normalized across each I_p,rel of the incremental cycle test and the ML1 method was applied.

Nonlinear parameter Lyapunov-exponents.

As mentioned earlier, λ_max characterizes the system’s sensitivity to local perturbations [16]. By calculating the divergence rate of points in phase space, the stability of the system is quantified. According to [20], the average divergence of neighboring points in phase space for discrete time steps can be described as: (4) where d_j(i) is the distance between the jth pair of nearest neighbors after i discrete time-steps, i.e., i ⋅ Δt seconds, with Δt denoting the sampling period of the time series. C_j describes the initial separation and λ_max depicts the growth rate.

In biomechanical analysis, experimental time series are investigated, thus the true maximum Lyapunov exponent cannot be determined—unlike fully describable systems of differential equations. λ_max is therefore defined as the slope of the linear range of the divergence curve. When taking the logarithm on both sides of Eq 4 and using a least-square fit to the average line over all values of j, denoted by 〈 〉, the largest Lyapunov exponent is defined by: (5)λ_max ≈0 corresponds to a stable system. If λ_max is positive, the system is divergent and thus unstable [16].

For experimental data, a distinction is made between the short-term and long-term divergence exponents. When observing the divergence curve in gait analysis, the short-term Lyapunov exponent λ_st represents the slope of the nearly linear range from the beginning up to the first step or stride. The long-term Lyapunov exponent λ_lt represents the gradient from the middle to final range (approximately 4th to 10th stride in gait). For the current analysis in cycling, we consider a stride analogous to a revolution of the crank as one cycle. The two Lyapunov measures should not be interpreted in the same way. According to Terrier and Reynard [21], in gait analysis, λ_st can be considered as a stability measure to prevent falls, while λ_lt is a measure of (gait) complexity.

During the present data analysis, λ_st and λ_lt were evaluated using the R-package tseriesChaos [22]. In a first step, state spaces with corresponding time delays τ were calculated using the Average Mutual Information (AMI) method. Compared to the autocorrelation function, the AMI also considers nonlinear correlations [23]. Different AMI lags running from 0 to 150 are shown in Fig 2A for all included time series. The first minimum was considered the appropriate time delay [16]. The distribution of all first minima is shown in the histogram in Fig 2B. For further calculation, τ = 29 was assumed for all participants and relative power intervals. To generate a suitable state space, an appropriate embedding dimension was needed in addition to τ. To find an optimal embedding dimension, we applied the False Nearest Neighbors (FNN) method. Figuratively, close points on an unfolded attractor can be far away from each other in the next higher dimension and are defined as False Nearest Neighbors [24]. The proportion of these FNN in one dimension is determined within the applied function. Close time separations or distances between two points are highly correlated and must be excluded when determining FNN. In simple terms, this is defined by the Theiler window, which was set to 3*τ [25]. The diameter of the neighborhood in which (False) Neighbors should be found was set to half the standard deviation of the time series [23]. Based on these FNN calculations (see Fig 2C), an embedding dimension of 4 was specified for all participants and intervals.

Download:

Fig 2. Method for calculating Lyapunov-exponent.

A: Average Mutual Information (AMI) for all time series to determine overall delay. B: Distribution of considered delays as first AMI minimum. Overall delay was set to τ = 29. C: Calculation of False Nearest Neighbors to define overall embedding dimension, which was stated as ed = 4. D: Logarithmic divergence for calculating the short-term Lyapunov exponent as the slope for the first 100 points (λ_st) and long-term Lyapunov exponent as slope for 400 to 1000 points (λ_lt). Short-term and long-term intercepts were also considered.

https://doi.org/10.1371/journal.pone.0285408.g002

Finally, the divergence and the regression coefficients for λ_st and λ_lt were calculated with the functions lyap_k() and lyap(). Similar to gait analysis [13–15, 26–28], λ_st was determined as the slope of the exponential divergence for half a cycle (approximately 50 points). λ_lt was calculated as the slope of the divergence for cycles 4–10 (approximately points 400–1000). This method is shown in Fig 2D. The regression coefficients of the short-term and long-term divergence also include the intercepts ι_st and ι_lt. Lyapunov parameters were calculated considering the whole 3-min time series of the I_p,rel.

To obtain a common measure of the time series’ dispersion, the standard deviation SD was determined for each I_p,rel.

Statistical analysis

The calculated parameters (ML1_F, ML1_α, λ_st, λ_lt,ι_st, ι_lt, SD) were analyzed descriptively and tested for linear increase across the three I_p,rel using a contrast analysis. The contrast analysis is a special case of the analysis of variance. The alternative hypothesis is formulated via contrasts (i.e. lambda weights). A linear increase over the three individual relative power intervals I_p,rel was assumed for all determined parameters. Consequently, lambda weights of λ = {−1, 0, 1} were chosen. Contrast analyses of the parameters were performed using the calc_contrast()-function of the package cofad [29] for within-subject designs. Assumptions of normal distribution and homoscedasticity were assessed using boxplots. Symmetric boxes were interpreted as normally distributed data and interquartile ranges of equal length were interpreted as confirmation of variance homogeneity between groups.

Results

Incremental cycling test

Results of the incremental cycling test are shown in Table 2, reporting power outputs (in W per kg bodyweight), including means and standard deviations across the I_p,rel. The absolute mean values of the power outputs increased consistently over the relative power intervals and the standard deviations remained approximately constant. The corresponding table including raw power and the figure with individual depicted development across the I_p,rel can be found in the supporting information (see: S1 Fig and S1 Table).

Download:

Table 2. Power outputs across relative power intervals.

https://doi.org/10.1371/journal.pone.0285408.t002

Arithmetic mean and standard deviations across three relative power intervals (I_p,25, I_p,50, I_p,75) of respective nonlinear parameters are shown in Figs 3–5. Single participants are displayed as dashed lines. All necessary requirements of the statistical tests were met. Detailed results of the contrast analysis are also presented in Table 3.

Download:

Fig 3. Results of the ML1 calculation.

Individual participants are displayed as dashed lines. A: ML1_F across relative power intervals using crank force data. B: ML1_α across relative power intervals using crank acceleration data.

https://doi.org/10.1371/journal.pone.0285408.g003

Download:

Fig 4. Results of the Lyapunov calculation.

A: Short-term Lyapunov exponents as slope λ_st across relative power intervals show a slight increasing trend. B: Long-term Lyapunov exponents as slope λ_lt across relative power intervals show no trend. C: Short-term Lyapunov intercepts ι_st show an increasing trend across relative power intervals. D: Long-term Lyapunov intercepts ι_lt also show an increasing trend across relative power intervals.

https://doi.org/10.1371/journal.pone.0285408.g004

Download:

Fig 5. Results of the calculation of the standard deviation.

An increase of the standard deviation across relative power intervals was present for 5 out of 10 participants.

https://doi.org/10.1371/journal.pone.0285408.g005

Download:

Table 3. Results of the contrast analysis.

https://doi.org/10.1371/journal.pone.0285408.t003

ML1

The results of ML1_F are shown in Fig 3A. A linear increase was found for mean ML1_F across the relative power intervals for all but one participant. Examining all participants, an interindividual variability across the relative power intervals was present, but remained fairly stable with a standard deviation of 0.042, 0.048, and 0.060. Further, the contrast analysis confirmed the presumed linear increase of ML1_F across the steps of the three I_p,rel (p-value = 8 × 10⁻⁴, t-value = 4.489) with a large effect (g-effect = 1.420). The same tendency was observed for ML1_α (Fig 3B). However, the standard deviation increased with increasing load level. The contrast analysis also confirmed a statistically significant linear increase between ML1_α and relative power intervals with p = 0.0151 (t-value = 2.571 and a large effect g = 0.813).

Lyapunov parameter

The four Lyapunov parameters λ_st, λ_lt, ι_st and ι_lt across the individual power interval are shown in Fig 4. No linear increase was found for the divergence parameter λ_lt (Fig 4B). In addition, the contrast analysis did not show significant linear differences between the relative power intervals (p_λ,lt = 0.421). The effect and test sizes were small (g_λ,lt = 0.065 and t_λ,lt = 0.206).

For λ_st (see Fig 4A), however, a linear trend is present ( = 0.028, t_λ,st = 2.196 and g_λ,st = 0.695). A distinct linear trend is visible for both intercepts ι_st and ι_lt in Fig 4C and 4D. This trend is also statistically significant with a higher effect g for ι_st (p_ι,st = 2 × 10⁻⁵ with g_ι,st = 2.327 and p_ι,lt = 1 × 10⁻⁴ with g_ι,lt = 1.845).

Standard deviation

As shown in Fig 5, there was no clear increase of the standard deviation across relative power intervals. Even though 5 out of 10 participants showed an increasing trend across all three intervals, this was not confirmed for the other 5 participants. Consequently, the contrast analysis for a linear relationship between the relative power intervals and the standard deviation shows no significant result (p = 0.0563).

Discussion

The present paper aimed to answer two main research questions: 1. To what extent can the underlying force data of the previously established ML1_F be substituted by kinematic data to calculate ML1_α and still distinguish individual relative power intervals in cycling? 2. Can other variables, such as Lyapunov parameters, determined via kinematic data also distinguish individual relative power intervals? To answer these questions, an acceleration time series was calculated from the output of a crank-based gyroscope in a cycling laboratory step test. Previous results from laboratory and field studies were confirmed: a statistically significant, linear relationship between ML1_F and the relative power intervals was found although the effects in the current step test were slightly smaller. This indicates a very good reproducibility of ML1_F even under different experimental conditions. Additionally, a linear increase for ML1_α using acceleration data generated from the gyroscope was shown across the three increasing individual load levels. This shows the robustness of ML1 towards various signals, as ML1_α seems able to distinguish load levels. The use of a simpler gyroscope or accelerometer on the crank appears justified.

Furthermore, statistically significant linear increase for the short-term Lyapunov exponent λ_st, as well as for short and long-term intercepts of divergences ι_st and ι_lt were found for increasing individual load levels I_p,25, I_p,50 and I_p,75. The same results were not found for the long-term Lyapunov exponent λ_lt. This supports the idea that increasing load during cycling can be described using Lyapunov exponents calculated from kinematic data.

Nevertheless, there are some limiting factors that should be mentioned. λ_st is calculated as a possible surrogate for local system stability, while λ_lt is more related to the complexity of the attractor of the corresponding movement [21]. Transferred to the present study, decreasing local motion stability would have to be assumed with increasing load during cycling. However, the complexity of the attractor in phase space does not increase. At the same time, the intercepts ι_st and ι_lt show a significant linear increase with increasing load. While the slope (λ) is a measure of change of divergence in phase space, the intercept ι rather describes the ‘initial distance’ or offset of diverging points in this phase space. This could be important within the present calculation of the Lyapunov parameters using the implementation of [22], since the radius for where to find the nearest neighbors (eps) has to be defined. In the presented algorithm, this radius was set to half the standard deviation of the corresponding time series. As shown in Fig 5, half the participants show an increase over the relative power intervals. Therefore, the following applies for these participants: The larger the standard deviation, the larger the radius in which the nearest neighbors can be found. However, with a fixed assumed eps (half mean of overall SD) for all conditions, the linear tendency of the intercepts still occurs. Thus, this does not explain the linear increase of ι.

Another point to mention is that methodological differences in calculating Lyapunov exponents lead to a lack of comparability [21, 30]. The methodology used in this analysis was transferred from gait analysis and applied to crank acceleration data during cycling. Even though our initial hypothesis was confirmed, the methodology should be adapted and standardized. In this way, different conditions, such as field vs. laboratory settings, may be more comparable.

Moreover, the pedaling motion is highly constrained, in particular when using an ergometer, as the crank is (nearly) exclusively accelerated in the sagittal plane. Without any axial inclination and a constant crank length in the static position on the ergometer, the motion is restricted to a static circle. Nevertheless, oscillations and vibrations may even occur for a static ergometer, and may have an impact on the results. To exclude such artefacts, a frequency analysis with a Fast Fourier Transformation was performed for each time series of the relative power intervals. No indications for unexplainable peaks could be found. Since the participants executed a cadence between 60 and 90 cycles per minute, a frequency peak can be distinguished at approximately 1 to 1.3 Hz (see Fig 6). Using a 42 teeth chainring, another frequency peak at about 50 to 55 Hz can be explained by multiplying the number of teeth by the cadence. No other unexplainable characteristic frequencies were found. The mentioned frequencies were eliminated using the low-pass filter.

Download:

Fig 6. Spectral analysis.

Exemplary time series with labeled characteristic frequencies.

https://doi.org/10.1371/journal.pone.0285408.g006

The high number of restrictions in the current laboratory study limits the transfer of our results to field applications. In the field, various confounding environmental variables, such as wind, uneven surfaces, changes in direction, and acceleration, etc., can affect the kinematic signal, in particular through different vibrations. When calculating Lyapunov exponents, these influences add to the complexity of the system. It is impossible to say whether and in what manner these conditions affect the results of the Lyapunov exponents. Acceleration signals and ML1_α may be more exposed to external factors than force data and ML1_F. This may also be the reason for the clearer linear increase of ML1_F compared to ML1_α in the results even in the current laboratory test. This is also suggested in [11], where a replication of ML1_F within indoor and outdoor tests was shown—despite external disturbance factors and additional vibrations. However, the disadvantage of ML1_α in field applications remains hypothetical at this point and should be further investigated in field studies.

With respect to future applications of the present findings, further consideration should be given to the extent to which they can be used to optimize the algorithms of e-bike propulsion. Currently, every available time series was considered in its entirety to maintain the system’s stability. These 3-minutes intervals contained approximately 255 cycles, i.e. 255 revolutions of the crank, and were thus too long for a real-time application. A sensitivity analysis should be performed to examine the required length of the time series that is needed to still distinguish relative power intervals in the same way as when considering the entire time series. The numbers of cycles should be determined from the middle of the corresponding relative power interval, in order to avoid possible short-term adjustments to load changes. According to a first assessment, a length of about 75 to 100 cycles seems to be necessary. This corresponds to a slightly longer time series than the length proposed in [11] for ML1_F. Nevertheless, further research is needed. This also applies to the generalizability of our findings, since our sample was fairly small and represented a rather homogeneous group of participants (males only).

Conclusions

To efficiently optimize e-bike algorithms, it is beneficial to determine the individual load of the cyclist using direct mechanical bicycle data, so there is no need for additional devices. Our study found that nonlinear dynamical methods are suitable for detecting individual load in cycling using mechanical kinematic data. We observed a linear increase in ML1, the short term Lyapunov Exponent λ_st and slopes (ι_st and ι_lt) across three relative power intervals based on data from a laboratory ergometer step test.

For both parameters, intraindividual differences were shown across three relative power intervals, but with notable interindividual differences of the participants within the same relative power interval. Thus, there is a noticeable trend for individual participants across increasing load levels, but with different respective base levels. However, the interindividual difference is not related to the individually achieved absolute power output of each participant in the ergometer step test. The nonlinear methods used in this context may accurately describe the individual load.

However, there are limitations to the transferability of these findings, as the data from the gyroscope used in this study were collected in a laboratory setting. In field conditions, or simply when not cycling on an ergometer, there may be additional influences (such as stronger vibrations), that affect the system stability. The high variance of ML1_α at the highest relative power interval (I_p,75) may be an indication of the presence of increasing interactions of hidden influences, even since no unexplained frequencies were found in the FFT. Further field tests are needed to determine the extent to which such influences may affect the system stability described by λ_st.

In summary, results of this paper show that λ_st in addition to the ML1_F may be useful to optimize of e-bike algorithms with respect to individual load based on kinematic data. However, further field tests and analysrs of real-time capability are needed to validate these findings.

Supporting information

S1 Fig. Absolute power output over relative power intervals.

The absolute mean values of the power outputs in W/kg increase equally over the relative power intervals and the variances remain constant.

https://doi.org/10.1371/journal.pone.0285408.s001

(TIFF)

S2 Fig. Discriminating statistics of the surrogate and original data.

For each participant, the statistic value of the original data is outside the values of the surrogate data.

https://doi.org/10.1371/journal.pone.0285408.s002

(TIFF)

S1 Table. Absolute power outputs across relative power intervals in W and W/kg.

https://doi.org/10.1371/journal.pone.0285408.s003

(PDF)

Acknowledgments

The authors would like to thank Lisa Peterson for linguistic correction as a native speaker.

References

1. Kautz AS, Feltner ME, Coyle EF, Baylor AM. The pedaling technique of elite endurance cyclists: changes with increasing workload at constant cadence. Journal of applied biomechanics. 1991;7.1:29–53.
- View Article
- Google Scholar
2. Hansen EA, Jørgensen LV, Jensen K, Fregly BJ, Sjøgaard G. Crank inertial load affects freely chosen pedal rate during cycling. Journal of biomechanics. 2002;35(2):277–285. pmid:11784546
- View Article
- PubMed/NCBI
- Google Scholar
3. Bertucci W, Grappe F, Girard A, Betik A, Rouillon JD. Effects on the crank torque profile when changing pedalling cadence in level ground and uphill road cycling. Journal of biomechanics. 2005;38(5):1003–1010. pmid:15797582
- View Article
- PubMed/NCBI
- Google Scholar
4. Emanuele U, Denoth J. Influence of road incline and body position on power-cadence relationship in endurance cycling. European journal of applied physiology. 2012;112(7):2433–2441. pmid:22045414
- View Article
- PubMed/NCBI
- Google Scholar
5. Rossato M, Bini RR, Carpes FP, Diefenthaeler F, Moro ARP. Cadence and workload effects on pedaling technique of well-trained cyclists. International journal of sports medicine. 2008;29(9):746–752. pmid:18302076
- View Article
- PubMed/NCBI
- Google Scholar
6. van Emmerik REA, Ducharme SW, Amado AC, Hamill J. Comparing dynamical systems concepts and techniques for biomechanical analysis. Journal of sport and health science. 2016;5(1):3–13. pmid:30356938
- View Article
- PubMed/NCBI
- Google Scholar
7. Quintana-Duque JC, Saupe D. Evidence of Chaos in Indoor Pedaling Motion using Non-linear Methods. In: Peters DM, O’Donoghue P, editors. Performance Analysis of Sport IX. Hoboken: Taylor and Francis; 2013. p. 203–214.
8. Padulo J, Di Capua R, Viggiano D. Pedaling time variability is increased in dropped riding position. European journal of applied physiology. 2012;112(8):3161–3165. pmid:22183087
- View Article
- PubMed/NCBI
- Google Scholar
9. Herrero-Molleda A, Flora P, Garcia-Lopez J. The influence of cyclists’ competitive level and exercise intensity on crank torque variability during pedalling. ISBS Proceedings Archive. 2021;39(1):364.
- View Article
- Google Scholar
10. Kunert A, Ott M, Reuter T, Koska D, Maiwald C. Phase space methods for non-linear analysis of pedalling forces in cycling. PloS one. 2019;14(4):e0198914. pmid:30998746
- View Article
- PubMed/NCBI
- Google Scholar
11. Kunert A, Koska D, Maiwald C. Investigation of the pedaling forces in a cycling field test using the non-linear measure ML1. Proceedings of the Institution of Mechanical Engineers, Part P: Journal of Sports Engineering and Technology. 2020;2:175433712097595.
- View Article
- Google Scholar
12. Warlop TB, Bollens B, Crevecoeur F, Detrembleur C, Lejeune TM. Dynamics of revolution time variability in cycling pattern: voluntary intent can alter the long-range autocorrelations. Annals of biomedical engineering. 2013;41(8):1604–1612. pmid:23712680
- View Article
- PubMed/NCBI
- Google Scholar
13. Dingwell JB, Cusumano JP, Cavanagh PR, Sternad D. Local dynamic stability versus kinematic variability of continuous overground and treadmill walking. Journal of biomechanical engineering. 2001;123(1):27–32. pmid:11277298
- View Article
- PubMed/NCBI
- Google Scholar
14. Dingwell JB, Marin LC. Kinematic variability and local dynamic stability of upper body motions when walking at different speeds. Journal of biomechanics. 2006;39(3):444–452. pmid:16389084
- View Article
- PubMed/NCBI
- Google Scholar
15. Bruijn SM, Meijer OG, Beek PJ, van Dieën JH. Assessing the stability of human locomotion: a review of current measures. Journal of the Royal Society, Interface. 2013;10(83). pmid:23516062
- View Article
- PubMed/NCBI
- Google Scholar
16. Kantz H, Schreiber T. Nonlinear time series analysis. Cambridge university press; 2004.
17. R Core Team. R: A Language and Environment for Statistical Computing; 2020. Available from: https://www.R-project.org/.
18. Schreiber T, Schmitz A. Surrogate time series. Physica D: Nonlinear Phenomena. 2000;142(3):346–382.
- View Article
- Google Scholar
19. Garcia CA. nonlinearTseries: Nonlinear Time Series Analysis; 2020. Available from: https://CRAN.R-project.org/package=nonlinearTseries.
20. Rosenstein MT, Collins JJ, de Luca CJ. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena. 1993;65(1-2):117–134.
- View Article
- Google Scholar
21. Terrier P, Reynard F. Maximum Lyapunov exponent revisited: Long-term attractor divergence of gait dynamics is highly sensitive to the noise structure of stride intervals. Gait & posture. 2018;66:236–241.
- View Article
- Google Scholar
22. Antonio, Fabio Di Narzo. tseriesChaos: Analysis of Nonlinear Time Series; 2019. Available from: https://CRAN.R-project.org/package=tseriesChaos.
23. Hegger R, Kantz H, Schreiber T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos (Woodbury, NY). 1999;9(2):413–435. pmid:12779839
- View Article
- PubMed/NCBI
- Google Scholar
24. Kennel MB, Brown R, Abarbanel HDI. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical review. 1992;45(6):3403. pmid:9907388
- View Article
- PubMed/NCBI
- Google Scholar
25. Theiler J. Estimating fractal dimension. Journal of the Optical Society of America A. 1990;7(6):1055.
- View Article
- Google Scholar
26. Dingwell JB, Cusumano JP. Nonlinear time series analysis of normal and pathological human walking. Chaos (Woodbury, NY). 2000;10(4):848–863. pmid:12779434
- View Article
- PubMed/NCBI
- Google Scholar
27. Kang HG, Dingwell JB. Intra-session reliability of local dynamic stability of walking. Gait & posture. 2006;24(3):386–390. pmid:16413784
- View Article
- PubMed/NCBI
- Google Scholar
28. Mehdizadeh S, Arshi AR, Davids K. Effect of speed on local dynamic stability of locomotion under different task constraints in running. European journal of sport science. 2014;14(8):791–798. pmid:24720520
- View Article
- PubMed/NCBI
- Google Scholar
29. Burkhardt M, Titz J. cofad: Contrast Analyses for Factorial Designs; 2020. Available from: https://CRAN.R-project.org/package=cofad.
30. Mehdizadeh S. The largest Lyapunov exponent of gait in young and elderly individuals: A systematic review. Gait & posture. 2018;60:241–250.
- View Article
- Google Scholar

[ref1] 1. Kautz AS, Feltner ME, Coyle EF, Baylor AM. The pedaling technique of elite endurance cyclists: changes with increasing workload at constant cadence. Journal of applied biomechanics. 1991;7.1:29–53.
View Article
Google Scholar

[2] View Article

[3] Google Scholar

[ref2] 2. Hansen EA, Jørgensen LV, Jensen K, Fregly BJ, Sjøgaard G. Crank inertial load affects freely chosen pedal rate during cycling. Journal of biomechanics. 2002;35(2):277–285. pmid:11784546
View Article
PubMed/NCBI
Google Scholar

[5] View Article

[6] PubMed/NCBI

[7] Google Scholar

[ref3] 3. Bertucci W, Grappe F, Girard A, Betik A, Rouillon JD. Effects on the crank torque profile when changing pedalling cadence in level ground and uphill road cycling. Journal of biomechanics. 2005;38(5):1003–1010. pmid:15797582
View Article
PubMed/NCBI
Google Scholar

[9] View Article

[10] PubMed/NCBI

[11] Google Scholar

[ref4] 4. Emanuele U, Denoth J. Influence of road incline and body position on power-cadence relationship in endurance cycling. European journal of applied physiology. 2012;112(7):2433–2441. pmid:22045414
View Article
PubMed/NCBI
Google Scholar

[13] View Article

[14] PubMed/NCBI

[15] Google Scholar

[ref5] 5. Rossato M, Bini RR, Carpes FP, Diefenthaeler F, Moro ARP. Cadence and workload effects on pedaling technique of well-trained cyclists. International journal of sports medicine. 2008;29(9):746–752. pmid:18302076
View Article
PubMed/NCBI
Google Scholar

[17] View Article

[18] PubMed/NCBI

[19] Google Scholar

[ref6] 6. van Emmerik REA, Ducharme SW, Amado AC, Hamill J. Comparing dynamical systems concepts and techniques for biomechanical analysis. Journal of sport and health science. 2016;5(1):3–13. pmid:30356938
View Article
PubMed/NCBI
Google Scholar

[21] View Article

[22] PubMed/NCBI

[23] Google Scholar

[ref7] 7. Quintana-Duque JC, Saupe D. Evidence of Chaos in Indoor Pedaling Motion using Non-linear Methods. In: Peters DM, O’Donoghue P, editors. Performance Analysis of Sport IX. Hoboken: Taylor and Francis; 2013. p. 203–214.

[ref8] 8. Padulo J, Di Capua R, Viggiano D. Pedaling time variability is increased in dropped riding position. European journal of applied physiology. 2012;112(8):3161–3165. pmid:22183087
View Article
PubMed/NCBI
Google Scholar

[26] View Article

[27] PubMed/NCBI

[28] Google Scholar

[ref9] 9. Herrero-Molleda A, Flora P, Garcia-Lopez J. The influence of cyclists’ competitive level and exercise intensity on crank torque variability during pedalling. ISBS Proceedings Archive. 2021;39(1):364.
View Article
Google Scholar

[30] View Article

[31] Google Scholar

[ref10] 10. Kunert A, Ott M, Reuter T, Koska D, Maiwald C. Phase space methods for non-linear analysis of pedalling forces in cycling. PloS one. 2019;14(4):e0198914. pmid:30998746
View Article
PubMed/NCBI
Google Scholar

[33] View Article

[34] PubMed/NCBI

[35] Google Scholar

[ref11] 11. Kunert A, Koska D, Maiwald C. Investigation of the pedaling forces in a cycling field test using the non-linear measure ML1. Proceedings of the Institution of Mechanical Engineers, Part P: Journal of Sports Engineering and Technology. 2020;2:175433712097595.
View Article
Google Scholar

[37] View Article

[38] Google Scholar

[ref12] 12. Warlop TB, Bollens B, Crevecoeur F, Detrembleur C, Lejeune TM. Dynamics of revolution time variability in cycling pattern: voluntary intent can alter the long-range autocorrelations. Annals of biomedical engineering. 2013;41(8):1604–1612. pmid:23712680
View Article
PubMed/NCBI
Google Scholar

[40] View Article

[41] PubMed/NCBI

[42] Google Scholar

[ref13] 13. Dingwell JB, Cusumano JP, Cavanagh PR, Sternad D. Local dynamic stability versus kinematic variability of continuous overground and treadmill walking. Journal of biomechanical engineering. 2001;123(1):27–32. pmid:11277298
View Article
PubMed/NCBI
Google Scholar

[44] View Article

[45] PubMed/NCBI

[46] Google Scholar

[ref14] 14. Dingwell JB, Marin LC. Kinematic variability and local dynamic stability of upper body motions when walking at different speeds. Journal of biomechanics. 2006;39(3):444–452. pmid:16389084
View Article
PubMed/NCBI
Google Scholar

[48] View Article

[49] PubMed/NCBI

[50] Google Scholar

[ref15] 15. Bruijn SM, Meijer OG, Beek PJ, van Dieën JH. Assessing the stability of human locomotion: a review of current measures. Journal of the Royal Society, Interface. 2013;10(83). pmid:23516062
View Article
PubMed/NCBI
Google Scholar

[52] View Article

[53] PubMed/NCBI

[54] Google Scholar

[ref16] 16. Kantz H, Schreiber T. Nonlinear time series analysis. Cambridge university press; 2004.

[ref17] 17. R Core Team. R: A Language and Environment for Statistical Computing; 2020. Available from: https://www.R-project.org/.

[ref18] 18. Schreiber T, Schmitz A. Surrogate time series. Physica D: Nonlinear Phenomena. 2000;142(3):346–382.
View Article
Google Scholar

[58] View Article

[59] Google Scholar

[ref19] 19. Garcia CA. nonlinearTseries: Nonlinear Time Series Analysis; 2020. Available from: https://CRAN.R-project.org/package=nonlinearTseries.

[ref20] 20. Rosenstein MT, Collins JJ, de Luca CJ. A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena. 1993;65(1-2):117–134.
View Article
Google Scholar

[62] View Article

[63] Google Scholar

[ref21] 21. Terrier P, Reynard F. Maximum Lyapunov exponent revisited: Long-term attractor divergence of gait dynamics is highly sensitive to the noise structure of stride intervals. Gait & posture. 2018;66:236–241.
View Article
Google Scholar

[65] View Article

[66] Google Scholar

[ref22] 22. Antonio, Fabio Di Narzo. tseriesChaos: Analysis of Nonlinear Time Series; 2019. Available from: https://CRAN.R-project.org/package=tseriesChaos.

[ref23] 23. Hegger R, Kantz H, Schreiber T. Practical implementation of nonlinear time series methods: The TISEAN package. Chaos (Woodbury, NY). 1999;9(2):413–435. pmid:12779839
View Article
PubMed/NCBI
Google Scholar

[69] View Article

[70] PubMed/NCBI

[71] Google Scholar

[ref24] 24. Kennel MB, Brown R, Abarbanel HDI. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Physical review. 1992;45(6):3403. pmid:9907388
View Article
PubMed/NCBI
Google Scholar

[73] View Article

[74] PubMed/NCBI

[75] Google Scholar

[ref25] 25. Theiler J. Estimating fractal dimension. Journal of the Optical Society of America A. 1990;7(6):1055.
View Article
Google Scholar

[77] View Article

[78] Google Scholar

[ref26] 26. Dingwell JB, Cusumano JP. Nonlinear time series analysis of normal and pathological human walking. Chaos (Woodbury, NY). 2000;10(4):848–863. pmid:12779434
View Article
PubMed/NCBI
Google Scholar

[80] View Article

[81] PubMed/NCBI

[82] Google Scholar

[ref27] 27. Kang HG, Dingwell JB. Intra-session reliability of local dynamic stability of walking. Gait & posture. 2006;24(3):386–390. pmid:16413784
View Article
PubMed/NCBI
Google Scholar

[84] View Article

[85] PubMed/NCBI

[86] Google Scholar

[ref28] 28. Mehdizadeh S, Arshi AR, Davids K. Effect of speed on local dynamic stability of locomotion under different task constraints in running. European journal of sport science. 2014;14(8):791–798. pmid:24720520
View Article
PubMed/NCBI
Google Scholar

[88] View Article

[89] PubMed/NCBI

[90] Google Scholar

[ref29] 29. Burkhardt M, Titz J. cofad: Contrast Analyses for Factorial Designs; 2020. Available from: https://CRAN.R-project.org/package=cofad.

[ref30] 30. Mehdizadeh S. The largest Lyapunov exponent of gait in young and elderly individuals: A systematic review. Gait & posture. 2018;60:241–250.
View Article
Google Scholar

[93] View Article

[94] Google Scholar

Figures

Abstract

Introduction

Materials and methods

Participants

Experimental setup

Data analysis

Assumptions and testing for nonlinearity.

Nonlinear parameter ML1.

Nonlinear parameter Lyapunov-exponents.

Statistical analysis

Results

Incremental cycling test

ML1

Lyapunov parameter

Standard deviation

Discussion

Conclusions

Supporting information

S1 Fig. Absolute power output over relative power intervals.

S2 Fig. Discriminating statistics of the surrogate and original data.

S1 Table. Absolute power outputs across relative power intervals in W and W/kg.

Acknowledgments

References