Front. Bioeng. Biotechnol. Frontiers in Bioengineering and Biotechnology Front. Bioeng. Biotechnol. 2296-4185 Frontiers Media S.A. 687951 10.3389/fbioe.2021.687951 Bioengineering and Biotechnology Original Research A Multivariate Polynomial Regression to Reconstruct Ground Contact and Flight Times Based on a Sine Wave Model for Vertical Ground Reaction Force and Measured Effective Timings Patoz et al. Reconstructed Contact and Flight Times Patoz Aurélien 1 2 * Lussiana Thibault 2 3 4 Breine Bastiaan 2 5 Gindre Cyrille 2 3 Malatesta Davide 1 Institute of Sport Sciences University of Lausanne, Lausanne, Switzerland Research and Development Department Volodalen Swiss Sport Lab, Aigle, Switzerland Research and Development Department Volodalen, Chavéria, France Research Unit EA3920 Prognostic Markers and Regulatory Factors of Cardiovascular Diseases and Exercise Performance Health Innovation Platform University of Franche-Comté, Besançon, France Department of Movement and Sports Sciences Ghent University, Ghent, Belgium

Edited by: Yang Liu, Hong Kong Polytechnic University, Hong Kong, SAR China

Reviewed by: Weiwei Yan, China Jiliang University, China

Shuo Chen, Tongji University, China

 *Correspondence: Aurélien Patoz, aurelien.patoz@unil.ch

This article was submitted to Biomechanics, a section of the journal Frontiers in Bioengineering and Biotechnology

 04 11 2021 2021 9 687951 30 03 2021 29 09 2021 Copyright © 2021 Patoz, Lussiana, Breine, Gindre and Malatesta. 2021 Patoz, Lussiana, Breine, Gindre and Malatesta

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Effective contact ( t c e ) and flight ( t f e ) times, instead of ground contact ( t c ) and flight ( t f ) times, are usually collected outside the laboratory using inertial sensors. Unfortunately, t c e and t f e cannot be related to t c and t f because the exact shape of vertical ground reaction force is unknown. However, using a sine wave approximation for vertical force, t c e and t c as well as t f e and t f could be related. Indeed, under this approximation, a transcendental equation was obtained and solved numerically over a t c e  x   t f e grid. Then, a multivariate polynomial regression was applied to the numerical outcome. In order to reach a root-mean-square error of 0.5 ms, the final model was given by an eighth-order polynomial. As a direct application, this model was applied to experimentally measured t c e values. Then, reconstructed t c (using the model) was compared to corresponding experimental ground truth. A systematic bias of 35 ms was depicted, demonstrating that ground truth t c values were larger than reconstructed ones. Nonetheless, error in the reconstruction of t c from t c e was coming from the sine wave approximation, while the polynomial regression did not introduce further error. The presented model could be added to algorithms within sports watches to provide robust estimations of t c and t f in real time, which would allow coaches and practitioners to better evaluate running performance and to prevent running-related injuries.

 running biomechanics sensors inertial measurement unit machine learning Introduction

Ground contact ( t c ) and flight ( t f ) times are key temporal parameters of running biomechanics. Indeed, Novacheck (1998) postulated that the presence of t f allowed distinguishing walking from running gaits. In other words, the duty factor (the ratio of t c over stride duration) is under 50% for running (Minetti, 1998; Folland et al., 2017). Moreover, t c was shown to be self-optimized to minimize the metabolic cost of running (Moore et al., 2019). These two parameters are obtained from foot-strike (FS) and toe-off (TO) events. More specifically, t c represents the time from FS to TO of the same foot, while t f is the time from TO of one foot to FS of the contralateral foot. Therefore, t c and t f rely on the accuracy of FS and TO detections, for which the use of force plates is considered the gold standard method. However, force plates could not always be available and used (Abendroth-Smith, 1996; Maiwald et al., 2009). In such case, alternatives would be to use a motion capture system (Lussiana et al., 2019; Patoz et al., 2020) or a light-based optical technology (Debaere et al., 2013). Nevertheless, even though these three systems can be used outside the laboratory (Purcell et al., 2006; Hébert-Losier et al., 2015; Ammann et al., 2016; Lussiana and Gindre, 2016), they suffer a lack of portability and are restricted to a specific and small capture volume, that is, they do not allow continuous temporal gait data collection throughout the entire training or race. To overcome such limitations, techniques to identify FS and TO events were developed using portative tools such as inertial measurement units (IMUs), which are easy to use, low cost, and suitable for field measurements and very practical to use in a coaching environment (Camomilla et al., 2018).

Different techniques to identify gait events are available and depend on the placement of the IMU on the human body (Moe-Nilssen, 1998; Lee et al., 2010; Flaction et al., 2013; Giandolini et al., 2014; Norris et al., 2014; Giandolini et al., 2016; Gindre et al., 2016; Falbriard et al., 2018; Falbriard et al., 2020). Among them, when the IMU is positioned near the sacrum, that is, close to the center of mass, the vertical acceleration signal can be used to determine effective contact ( t c e ) and flight ( t f e ) times, instead of t c and t f (Flaction et al., 2013; Gindre et al., 2016). To delineate these effective timings, the vertical force is calculated based on Newton’s second law using the body mass (m) of individuals and the vertical acceleration data. Then, these effective timings are based on effective FS (eFS) and effective TO (eTO) events. More precisely, eFS and eTO correspond to the instants of time where the vertical force increases above and decreases below body weight (mg), respectively (Cavagna et al., 1988). The authors (Flaction et al., 2013; Gindre et al., 2016) did not mention why a 20 N threshold was not used to determine FS and TO events from their IMU data, even though this is the reference when using force plates data for event detection (Smith et al., 2015). However, the vertical acceleration recorded by an IMU during t f is usually negative (Gindre et al., 2016), while a force plate measure gives exactly zero. Therefore, it could be suspected that a 20 N threshold would not be reliable to obtain FS and TO events when dealing with IMU data, while the time at which the vertical force is equal to body weight would be equivalent between IMU and force plate data.

Using effective timings or t c and t f provide the same step duration, that is, it is given by either the sum of t c and t f or t c e and t f e . Thus, this temporal information is not lost. As for the effect of running speed, t c e and t c both decrease with increasing running speed, even though the decrease is much more important for t c than t c e (Cavagna et al., 2008; Da Rosa et al., 2019). Concerning t f e and t f , their trend with increasing running speed is not similar. Indeed, t f e tends to slightly decrease, while t f increases almost up to a plateau with increasing running speed (Cavagna et al., 2008; Da Rosa et al., 2019). In addition, t c e and t f e cannot directly be related to t c and t f , the reason being that the fraction of time spends below body weight during t c depends on the shape of the vertical ground reaction force, which is not precisely known when using IMUs (see above). Thus, t c and t f , parameters that are directly related to them, for example, duty factor (Minetti, 1998; Folland et al., 2017), as well as parameters that can be estimated from them, for example, vertical oscillation and vertical stiffness (Morin et al., 2005), cannot be obtained. Hence, the assessment of running biomechanics is restricted when using t c e and t f e .

Nonetheless, the vertical ground reaction force can be approximated using a sine wave as F z ( t ) = F z , max ⁡ sin ( π t / t c ) , where, based on momentum conservation law, F z , max = m g π ( t f / t c + 1 ) / 2 (Alexander, 1989; Kram and Dawson, 1998; Dalleau et al., 2004; Morin et al., 2005). In such case, the vertical ground reaction force is symmetric around t c / 2 , which means that the time duration between FS and eFS as well as between eTO and TO, called t g in what follows, are the same. Thereby, under the sine wave assumption, t c and t f can be obtained from t c e and t f e using t c = t c e + 2 t g and t f = t f e − 2 t g , if t g is known. These timings and the sine wave vertical ground reaction force are depicted in Figure 1 for a typical running stride. Recognizing that F z ( t g ) = m g = F z , max ⁡ sin ( π t g t c ) , and using the definition of F z , max given before, the following equation is obtained: csc ( π t g t c e + 2 t g ) = π 2 ( t f e − 2 t g t c e + 2 t g + 1 ) , which could not be solved analytically for t g (transcendental equation; Supplementary File) using Mathematica v12.1 (Wolfram, Oxford, UK), that is, no closed-form solution exists. Therefore, a numerical solution is required for any pair of t c e and t f e . Ultimately, a mathematical modeling of t g over the numerical t c e x t f e grid could be performed, and its accuracy could be evaluated using advanced data analysis techniques like machine learning. Indeed, supervised machine learning models like linear regressions have been used to model relationships between biomechanical measures and clinical outcomes (Halilaj et al., 2018; Backes et al., 2020; Alcantara et al., 2021). However, to the best of our knowledge, no attempt to provide such a model equation for t g has been made so far.

Vertical ground reaction force ( F z ) under the sine wave approximation, peak vertical force ( F z , max ) , foot-strike (FS) and toe-off (TO) events together with their corresponding effective events (eFS and eTO), as well as contact ( t c ) , flight ( t f ) , effective contact ( t c e ) , and effective flight ( t f e ) times, and time to reach body weight ( t g ) , for a typical running stride. Noteworthy, step duration is the same when using effective or usual timings.

Hence, the purpose of this study was to obtain a mathematical modeling of t g under the sine wave approximation of the vertical ground reaction force so that t c and t f can be reconstructed from t c e and t f e . As a direct experimental application, the proposed model was applied to experimentally measured t c e values. Then, the reconstructed t c values were compared to their corresponding experimental ground truth (gold standard).

 Materials and Methods Numerical Analysis

Brent’s method (also known as van Wijngaarden Dekker Brent method) (Brent, 1973; Press et al., 1992) was used to find the zeros of Eq. 1 for any pair of t c e and t f e . The zero of interest for a given t c e and t f e pair was considered to lie between 0 and the minimum of Eq. 1, which was minimized using the Broyden Fletcher Goldfarb Shanno method (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970). The numerical analysis was carried out using t c e and t f e values varying between 2.5 and 505 ms and using a grid spacing of 7.5 ms (4,624 grid points). The grid limits were chosen due to the fact that running requires 1) both a ground contact and a flight phase, that is, t c e and t f e cannot be 0 and 2) t c belongs to the interval [ 100   m s ,   400   m s ] and t f belongs to the interval [ 0   m s , 250   m s ] (Cavagna et al., 2008; Da Rosa et al., 2019; Lussiana et al., 2019), and to include any atypical t c e and t f e pair, that is, atypical runners. Noteworthy, the justification of the grid spacing is provided in Appendix. The grid spacing was dependent on the error threshold set to the mathematical modeling.

 Mathematical Modeling Boundary Relationship Between <inline-formula id="inf106"> <mml:math id="m107"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>c</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula> and <inline-formula id="inf107"> <mml:math id="m108"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>f</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula>

The numerical analysis showed that a linear boundary relationship is present between t c e and t f e (see Results Figure 2), that is, there is no solution for t g if t f e is higher than a certain percentage of t c e . This boundary relationship was computed by extracting the boundary points, that is, the smallest existing t f e values for every t c e grid point (68 pair of points). Then, a linear regression using ordinary least square was performed on a training set consisting of 85% of the entire set of boundary points. The y-intercept of the fitted linear model was held fixed at 0, the reason being that a null t c e necessarily ensures a null t f e . The linear model was tested on the remaining 15% points (testing set) and evaluated using the coefficient of determination ( R 2 ) and root-mean-square error (RMSE).

Contour plots depicting a) the numerically calculated time ( t g ) necessary to reach body weight and B) the corresponding percentage (p tg ) of time under body weight during ground contact time ( t c ) . The numerical simulation assumed a sine wave model for vertical ground reaction force and was performed over the t c e x t f e grid.

 Modeling a <inline-formula id="inf122"> <mml:math id="m123"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mi>g</mml:mi> </mml:msub> </mml:mrow> </mml:math> </inline-formula> Surface as Function of <inline-formula id="inf123"> <mml:math id="m124"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>c</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula> and <inline-formula id="inf124"> <mml:math id="m125"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>f</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula>

The numerical analysis showed that t g could be described by a smoothly increasing surface when increasing t c e and t f e (see Results Figure 2). Therefore, a multivariate polynomial regression using ordinary least square was performed on a training set consisting of t g values corresponding to 85% of the points within the boundary limits (i.e., the non-discarded grid points). The regression was performed using polynomials of order 1 to 15 and including intercept and interaction terms. RMSE on the remaining 15% points (testing set) was computed for each fitted polynomial.

 Experimental Application Participant Characteristics

One hundred recreational runners (Honert et al., 2020), 75 males (age: 31 ± 8 years, height: 180 ± 6 cm, body mass: 70 ± 7 kg, and weekly running distance: 37 ± 24 km) and 25 females (age: 30 ± 7 years, height: 169 ± 5 cm, body mass: 61 ± 6 kg, and weekly running distance: 20 ± 14 km), voluntarily participated in the present study. For study inclusion, participants were required to be in good self-reported general health with no current or recent lower extremity injury (≤1 month), to run at least once a week, and to have an estimated maximal aerobic speed ≥14 km/h. The study protocol was approved by the Ethics Committee (CER-VD 2020–00334) and adhered to the latest Declaration of Helsinki of the World Medical Association.

 Experimental Procedure

After providing written informed consent, each participant performed a 7-min warm-up run on an instrumented treadmill (Arsalis T150—FMT-MED, Louvain-la-Neuve, Belgium). Speed was set to 9 km/h for the first 3 min and was then increased by 0.5 km/h every 30 s. This was followed, after a short break (<5 min), by three 1-min runs (9, 11, and 13 km/h) performed in a randomized order (1-min recovery between each run). 3D kinetic data were collected during the first 10 strides following the 30-s mark of running trials. All participants were familiar with running on a treadmill as part of their usual training program and wore their habitual running shoes.

 Data Collection

3D kinetic data (1,000 Hz) were collected using the force plate embedded into the treadmill and using Vicon Nexus software v2.9.3 (Vicon, Oxford, UK). The laboratory coordinate system was oriented such that x-, y-, and z-axes denoted mediallateral (pointing toward the right side of the body), posterioranterior, and inferiorsuperior axis, respectively. Ground reaction force (analog signal) was exported in .c3d format and processed in Visual3D Professional software v6.01.12 (C-Motion Inc, Germantown, MD, United States). 3D ground reaction force signal was low-pass–filtered at 20 Hz using a fourth-order Butterworth filter and down-sampled to 200 Hz to represent a sampling frequency corresponding to typical measurements recorded using a central inertial unit.

 Data Analysis

For each running trial, eFS and eTO events were identified within Visual3D by applying a body weight threshold to the z-component of the ground reaction force (Cavagna et al., 1988). More explicitly, eFS was detected at the first data point greater or equal to mg within a running step, while eTO was detected at the last data point greater or equal to mg within the same running step. t c e and t f e were defined as the time from eFS to eTO of the same foot and from eTO of one foot to eFS of the contralateral foot, respectively.

In addition, FS and TO events were also identified within Visual3D. These events were detected by applying a 20 N threshold to the z-component of the ground reaction force (Smith et al., 2015). More explicitly, FS was detected at the first data point greater or equal to 20 N within a running step, while TO was detected at the last data point greater or equal to 20 N within the same running step. t c and t f were defined as the time from FS to TO of the same foot and from TO of one foot to FS of the contralateral foot, respectively.

The recorded vertical ground reaction force permitted to precisely measure t c and t f as well as t c e and t f e . Then, each t c e and t f e pair was fed to the best multivariate polynomial model to compute t g , which ultimately allowed to obtain t c . An instrumented treadmill was used to measure t c e and t f e (gold standard), instead of an IMU to remove any potential measurement error that would come from the IMU itself. Hence, the error obtained when comparing the reconstructed t c (obtained using the mathematical model and   t c e and t f e ) to its corresponding experimental ground truth (obtained from FS and TO events) could solely be coming from the sine wave assumption and the mathematical modeling but not from the measurement of t c e and t f e .

 Statistical Analysis

All data are presented as mean ± standard deviation. The reconstructed t c values were compared to corresponding experimental ground truth t c values using a BlandAltman plot (Bland and Altman, 1995; Atkinson and Nevill, 1998). Noteworthy, as step time is conserved, differences between measured and reconstructed t f values depicted the opposite behavior compared with the differences between measured and reconstructed t c values.

Systematic bias, lower and upper limit of agreements, and 95% confidence intervals (CI) were computed as well as RMSE. The difference between reconstructed and ground truth t c values was quantified using Cohen’s d effect size and interpreted as very small, small, moderate, and large when |d| values were close to 0.01, 0.2, 0.5, and 0.8, respectively (Cohen, 1988). Statistical analysis was performed using Jamovi (v1.2, retrieved from https://www.jamovi.org), with the level of significance set at p ≤ 0.05.

 Results Numerical Analysis

The numerically calculated t g values over the t c e x t f e grid are provided in Figure 2A, while Figure 2B depicts the corresponding percentage of time ( p t g ) spent under body weight during t c , [ p t g = 100 ∗ 2 t g / ( t c e + 2 t g ) ] .

 Mathematical Modeling Boundary Relationship Between <inline-formula id="inf159"> <mml:math id="m160"> <mml:mrow> <mml:msub> <mml:mtext>t</mml:mtext> <mml:mrow> <mml:mtext>ce</mml:mtext> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula> and <inline-formula id="inf160"> <mml:math id="m161"> <mml:mrow> <mml:msub> <mml:mtext>t</mml:mtext> <mml:mrow> <mml:mtext>fe</mml:mtext> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula>

The linear regression gave the model (Eq. 2): t f e = 0.795   t c e .

Applying this model to the testing set led to an R 2 of = 99.9 % and RMSE of 3.2 ms. The linear regression, training, and testing sets as well as predicted values are depicted in Figure 3.

Boundary relationship between t c e and t f e . A linear regression (solid line) was obtained using 85% of the entire boundary points (training set, small gray dots) and validated on the remaining 15% points (testing set, large black dots). Predictions are given by the black circles and led to a root-mean-square error of 3.2 ms ( R 2 = 99.9 % ) .

 Modeling a <inline-formula id="inf166"> <mml:math id="m168"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mi>g</mml:mi> </mml:msub> </mml:mrow> </mml:math> </inline-formula> Surface as Function of <inline-formula id="inf167"> <mml:math id="m169"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>c</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula> and <inline-formula id="inf168"> <mml:math id="m170"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>f</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula>

The grid points which did not satisfy the previously obtained boundary relationship (Eq. 2) were discarded (1814 discarded points). RMSE computed for each multivariate polynomial regression (order 1–15) is depicted in Figure 4. The polynomial which provided an RMSE smaller than 0.5 ms was kept as the final model of choice (RMSE = 0.43 ms; R 2 = 99.99 % ) and corresponded to a polynomial model including up to eighth-order terms [ P 8 ( t c e , t f e ) , Eq. 3]. The coefficients ( α i j , where 0 ≤ i + j ≤ 8 ) of the multivariate polynomial model are given in Table 1. P 8 ( t c e , t f e ) = ∑ i = 0 8 ∑ j = 0 8 − i α i , j t c e i t f e j

Root-mean-square error (RMSE) computed on the testing set (15% points) for polynomial fits of order 1 to 15 performed on the training set (85% points). The red circle denotes the final model of choice, an eighth-order polynomial model (RMSE = 0.43 ms; R 2 = 99.99 % ), and the gray line depicts the RMSE threshold of 0.5 ms.

Coefficients ( α i j , where 0 ≤ i + j ≤ 8 ) of the eighth-order multivariate polynomial model given by Eq. 3.

j ( exponent   of   t fe )
0 1 2 3 4 5 6 7 8
j ( exponent   of   t f e ) 0 −5.17E-5 −6.18E-2 2.73E0 −4.41E1 3.532 −1.55E3 3.783 −4.83E3 2.513
1 2.84E-1 −1.41E1 2.64E2 −2.45E3 1.234 −3.38E4 4.834 −2.78E4
2 1.17E1 −3.12E2 3.91E3 −2.49E4 8.434 −1.43E5 9.534
3 8.26E1 −2.25E3 2.20E4 −1.01E5 2.155 −1.72E5
4 5.13E2 −9.73E3 6.68E4 −1.90E5 1.865
5 1.63E3 −2.32E4 9.82E4 −1.22E5
6 3.41E3 −2.76E4 4.65E4
7 3.15E3 −8.66E3
8 4.62E2

Noteworthy, the threshold on RMSE ensured an error smaller than 1 ms on the reconstructed t c . The differences between t g values computed numerically and using the eighth-order polynomial model for the testing set (15% points) are depicted in Figure 5.

Differences between t g values ( Δ t g ) computed numerically (Section 2) and using the eighth-order polynomial model for the testing set (15% points). A difference larger than 2 ms was depicted for only two points (green and yellow circles) in the testing set, which were close to the boundary limit.

 Experimental Application

Reconstructed t c values were compared to corresponding experimental ground truth t c values using a BlandAltman plot, which is depicted in Figure 6. A systematic positive bias of 34.3 ms (95% CI [33.8 ms, 34.7 ms]) was obtained. The lower and upper limits of agreements were 0.0 ms (95% CI [−0.8 ms, 0.8 ms]) and 68.6 ms (95% CI [67.8 ms, 69.3 ms]), respectively. The RMSE between reconstructed and measured t c was 38.5 ms (7.6%), and Cohen’s d effect size was large (d = 1.1).

BlandAltman plot comparing experimentally measured and reconstructed t c using the multivariate polynomial model given by Eq. 3, which reports a systematic bias of 34.3 ms (95% confidence intervals [33.8 ms, 34.7 ms]). Δ t c : measured t c − reconstructed t c and t c ¯ : average of measured and reconstructed t c .

 Discussion

The proposed eighth-order multivariate polynomial model (Eq. 3) could be used to obtain t c and t f when an IMU is used to measure t c e and t f e . Thereby, important parameters to assess running biomechanics such as duty factor (Lussiana et al., 2019; Patoz et al., 2020), as well as vertical oscillation and vertical stiffness (Morin et al., 2005), could be calculated more precisely. Having these parameters would allow coaches and practitioners to better evaluate running performance outside the laboratory such as in a coaching environment and during an entire training or race, and to prevent running-related injuries.

In the case where an algorithm based on effective timings is running on the fly to provide live feedbacks, such as in sports watches, one could simply add the proposed model in the end of the algorithm chain, right before computing the biomechanical outcomes. However, many operations should be performed in a very small amount of time, where the number of operations is directly related to the order of the polynomial. Indeed, knowing that the number of terms in an n th -order polynomial composed of two variables is given by C 2 n + 2 , then C 2 n + 2 − 3 calculations are required to compute the polynomial features, that is, t c e i and t f e i , where 2 ≤ i ≤ n . In addition, C 2 n + 2 − 1 multiplications and C 2 n + 2 − 1 additions are necessary to calculate t g . Therefore, such a large number of operations could be problematic for the small computing power available in sports watches. If this is really an issue, the order of the polynomial could be decreased. For instance, a third-order polynomial model gave an RMSE of 2.5 ms (Figure 4), which, depending on the application, might already be sufficient. In this case, the number of operations would be reduced from 130 (eighth order) to 25 (third order), leading to a 5 times speedup, assuming sequential calculations (no parallelization).

The multivariate polynomial model (Eq. 3) was applied to experimentally measured t c e values. These results permitted us to show that the experimental ground truth t c was, on average, 34.3 ms higher than the reconstructed one. Since the multivariate polynomial regression reported an RMSE of 0.43 ms, the large systematic bias obtained here was inherently due to the sine wave approximation of the vertical ground reaction force. To further justify the previous statement, the polynomial depicting the smallest RMSE, that is, the 14th-order polynomial (RMSE = 0.12 ms; Figure 7), was used to compute t c based on t c e . Doing so, the following results were obtained: RMSE = 38.6 ms (7.6%), d = 1.1 (large effect size), and systematic bias = 34.2 ms [95% CI (33.7 ms, 34.6 ms)]. Therefore, to go beyond the scope of this study, future research should focus on defining a more accurate model of the vertical ground reaction force. Indeed, the sine wave approximation constituted the main limitation of the novel multivariate polynomial model proposed in this study.

Root-mean-square error as a function of grid size ranging from 36 to 40,804 total points and for each polynomial regression (1st to 15th order). The red circle denotes RMSE corresponding to a polynomial (eighth order) chosen in Section 3.2 (0.43 ms), and the gray line depicts an RMSE threshold of 0.5 ms.

 Conclusion

To conclude, in the present study, an eighth-order multivariate polynomial model was constructed based on the numerical solution of the transcendental equation given by Eq. 1. The proposed model permitted to compute t c and t f from effective timings ( t c e   and   t f e ) using the sine wave approximation of the vertical ground reaction force. The model was chosen so that RMSE was smaller than 0.5 ms. Therefore, the error in the computation of t c and t f was coming from the sine wave approximation, while the polynomial regression did not introduce further error.

 Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

 Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committee (CER-VD 2020–00334). The patients/participants provided their written informed consent to participate in this study.

 Author Contributions

Conceptualization, AP, TL, CG, and DM; methodology: AP, TL, CG, and DM; investigation: AP, TL, and BB; formal analysis: AP and BB; writing—original draft preparation: AP; writing—review and editing: AP, TL, BB, CG, and DM; supervision: AP, TL, CG, and DM

 Funding

This study was supported by Innosuisse (grant no. 35793.1 IP-LS).

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

The authors warmly thank the participants for their time and cooperation.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2021.687951/full#supplementary-material

References Abendroth-Smith J. (1996). Stride Adjustments during a Running Approach toward a Force Plate. Res. Q. Exerc. Sport 67, 97–101. 10.1080/02701367.1996.10607930 Alcantara R. S. Day E. M. Hahn M. E. Grabowski A. M. (2021). Sacral Acceleration Can Predict Whole-Body Kinetics and Stride Kinematics across Running Speeds. PeerJ 9, e11199. 10.7717/peerj.11199 Alexander R. M. (1989). On the Synchronization of Breathing with Running in Wallabies (Macropusspp.) and Horses (Equus caballus). J. Zoolog. 218, 69–85. 10.1111/j.1469-7998.1989.tb02526.x Ammann R. Taube W. Wyss T. (2016). Accuracy of PARTwear Inertial Sensor and Optojump Optical Measurement System for Measuring Ground Contact Time during Running. J. Strength Conditioning Res. 30, 2057–2063. 10.1519/jsc.0000000000001299 Atkinson G. Nevill A. M. (1998). Statistical Methods for Assessing Measurement Error (Reliability) in Variables Relevant to Sports Medicine. Sports Med. 26, 217–238. 10.2165/00007256-199826040-00002 Backes A. Skejø S. D. Gette P. Nielsen R. Ø. Sørensen H. Morio C. (2020). Predicting Cumulative Load during Running Using Field‐based Measures. Scand. J. Med. Sci. Sports 30, 2399–2407. 10.1111/sms.13796 Bland J. M. Altman D. G. (1995). Comparing Methods of Measurement: Why Plotting Difference against Standard Method Is Misleading. The Lancet 346, 1085–1087. 10.1016/s0140-6736(95)91748-9 Brent R. P. (1973). Algorithms for Minimization without Derivatives. Englewood Cliffs, NJ: Prentice-Hall. Broyden C. G. (1970). The Convergence of a Class of Double-Rank Minimization Algorithms 1. General Considerations. IMA J. Appl. Math. 6, 76–90. 10.1093/imamat/6.1.76 Camomilla V. Bergamini E. Fantozzi S. Vannozzi G. (2018). Trends Supporting the In-Field Use of Wearable Inertial Sensors for Sport Performance Evaluation: A Systematic Review. Sensors 18, 873. 10.3390/s18030873 Cavagna G. A. Franzetti P. Heglund N. C. Willems P. (1988). The Determinants of the Step Frequency in Running, Trotting and Hopping in Man and Other Vertebrates. J. Physiol. 399, 81–92. 10.1113/jphysiol.1988.sp017069 Cavagna G. A. Legramandi M. A. Peyré-Tartaruga L. A. (2008). Old Men Running: Mechanical Work and Elastic Bounce. Proc. R. Soc. B. 275, 411–418. 10.1098/rspb.2007.1288 Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. Oxfordshire, England, UK: Routledge. Da Rosa R. G. Oliveira H. B. Gomeñuka N. A. Masiero M. P. B. Da Silva E. S. Zanardi A. P. J. (2019). Landing-takeoff Asymmetries Applied to Running Mechanics: A New Perspective for Performance. Front. Physiol. 10, 415. 10.3389/fphys.2019.00415 Dalleau G. Belli A. Viale F. Lacour J. R. Bourdin M. (2004). A Simple Method for Field Measurements of Leg Stiffness in Hopping. Int. J. Sports Med. 25, 170–176. 10.1055/s-2003-45252 Debaere S. Jonkers I. Delecluse C. (2013). The Contribution of Step Characteristics to Sprint Running Performance in High-Level Male and Female Athletes. J. Strength Conditioning Res. 27, 116–124. 10.1519/jsc.0b013e31825183ef Falbriard M. Meyer F. Mariani B. Millet G. P. Aminian K. (2018). Accurate Estimation of Running Temporal Parameters Using Foot-Worn Inertial Sensors. Front. Physiol. 9, 610. 10.3389/fphys.2018.00610 Falbriard M. Meyer F. Mariani B. Millet G. P. Aminian K. (2020). Drift-Free Foot Orientation Estimation in Running Using Wearable IMU. Front. Bioeng. Biotechnol. 8, 65. 10.3389/fbioe.2020.00065 Flaction P. Quievre J. Morin J. B. (2013). An Athletic Performance Monitoring Device. Washington, DC: U.S. Patent and Tradematk Office patent application. Fletcher R. (1970). A New Approach to Variable Metric Algorithms. Comp. J. 13, 317–322. 10.1093/comjnl/13.3.317 Folland J. P. Allen S. J. Black M. I. Handsaker J. C. Forrester S. E. (2017). Running Technique Is an Important Component of Running Economy and Performance. Med. Sci. Sports Exerc. 49, 1412–1423. 10.1249/mss.0000000000001245 Giandolini M. Horvais N. Rossi J. Millet G. Y. Samozino P. Morin J.-B. (2016). Foot Strike Pattern Differently Affects the Axial and Transverse Components of Shock Acceleration and Attenuation in Downhill Trail Running. J. Biomech. 49, 1765–1771. 10.1016/j.jbiomech.2016.04.001 Giandolini M. Poupard T. Gimenez P. Horvais N. Millet G. Y. Morin J.-B. (2014). A Simple Field Method to Identify Foot Strike Pattern during Running. J. Biomech. 47, 1588–1593. 10.1016/j.jbiomech.2014.03.002 Gindre C. Lussiana T. Hebert-Losier K. Morin J.-B. (2016). Reliability and Validity of the Myotest for Measuring Running Stride Kinematics. J. Sports Sci. 34, 664–670. 10.1080/02640414.2015.1068436 Goldfarb D. (1970). A Family of Variable-Metric Methods Derived by Variational Means. Math. Comp. 24, 23. 10.1090/s0025-5718-1970-0258249-6 Halilaj E. Rajagopal A. Fiterau M. Hicks J. L. Hastie T. J. Delp S. L. (2018). Machine Learning in Human Movement Biomechanics: Best Practices, Common Pitfalls, and New Opportunities. J. Biomech. 81, 1–11. 10.1016/j.jbiomech.2018.09.009 Hébert-losier K. Mourot L. Holmberg H.-C. (2015). Elite and Amateur Orienteers' Running Biomechanics on Three Surfaces at Three Speeds. Med. Sci. Sports Exerc. 47, 381–389. 10.1249/mss.0000000000000413 Honert E. C. Mohr M. Lam W.-K. Nigg S. (2020). Shoe Feature Recommendations for Different Running Levels: A Delphi Study. PLOS ONE 15, e0236047. 10.1371/journal.pone.0236047 Kram R. Dawson T. J. (1998). Energetics and Biomechanics of Locomotion by Red Kangaroos (Macropus rufus). Comp. Biochem. Physiol. B: Biochem. Mol. Biol. 120, 41–49. 10.1016/s0305-0491(98)00022-4 Lee J. B. Mellifont R. B. Burkett B. J. (2010). The Use of a Single Inertial Sensor to Identify Stride, Step, and Stance Durations of Running Gait. J. Sci. Med. Sport 13, 270–273. 10.1016/j.jsams.2009.01.005 Lussiana T. Patoz A. Gindre C. Mourot L. Hébert-Losier K. (2019). The Implications of Time on the Ground on Running Economy: Less Is Not Always Better. J. Exp. Biol. 222, jeb192047. 10.1242/jeb.192047 Lussiana T. Gindre C. (2016). Feel Your Stride and Find Your Preferred Running Speed. Biol. Open 5, 45–48. 10.1242/bio.014886 Maiwald C. Sterzing T. Mayer T. A. Milani T. L. (2009). Detecting Foot-To-Ground Contact from Kinematic Data in Running. Footwear Sci. 1, 111–118. 10.1080/19424280903133938 Minetti A. E. (1998). A Model Equation for the Prediction of Mechanical Internal Work of Terrestrial Locomotion. J. Biomech. 31, 463–468. 10.1016/s0021-9290(98)00038-4 Moe-Nilssen R. (1998). A New Method for Evaluating Motor Control in Gait under Real-Life Environmental Conditions. Part 1: The Instrument. Clin. Biomech. 13, 320–327. 10.1016/s0268-0033(98)00089-8 Moore I. S. Ashford K. J. Cross C. Hope J. Jones H. S. R. Mccarthy-Ryan M. (2019). Humans Optimize Ground Contact Time and Leg Stiffness to Minimize the Metabolic Cost of Running. Front. Sports Act Living 1, 53. 10.3389/fspor.2019.00053 Morin J.-B. Dalleau G. Kyröläinen H. Jeannin T. Belli A. (2005). A Simple Method for Measuring Stiffness during Running. J. Appl. Biomech. 21, 167–180. 10.1123/jab.21.2.167 Norris M. Anderson R. Kenny I. C. (2014). Method Analysis of Accelerometers and Gyroscopes in Running Gait: A Systematic Review. Proc. Inst. Mech. Eng. P: J. Sports Eng. Tech. 228, 3–15. 10.1177/1754337113502472 Novacheck T. F. (1998). The Biomechanics of Running. Gait & Posture 7, 77–95. 10.1016/s0966-6362(97)00038-6 Patoz A. Lussiana T. Thouvenot A. Mourot L. Gindre C. (2020). Duty Factor Reflects Lower Limb Kinematics of Running. Appl. Sci. 10, 8818. 10.3390/app10248818 Press W. H. Teukolsky S. A. Vetterling W. T. (1992). Numerical Recipes in FORTRAN: The Art of Scientific Computing. Cambridge, England: Cambridge University Press. Purcell B. Channells J. James D. Barrett R. (2006). Use of Accelerometers for Detecting Foot-Ground Contact Time during Running. Proc. SPIE - Int. Soc. Opt. Eng. 6036, 292–299. 10.1117/12.638389 Shanno D. F. (1970). Conditioning of Quasi-Newton Methods for Function Minimization. Math. Comp. 24, 647. 10.1090/s0025-5718-1970-0274029-x Smith L. Preece S. Mason D. Bramah C. (2015). A Comparison of Kinematic Algorithms to Estimate Gait Events during Overground Running. Gait & Posture 41, 39–43. 10.1016/j.gaitpost.2014.08.009 Appendix: Justification of the Choice of the <inline-formula id="inf213"> <mml:math id="m216"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>c</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> </mml:mrow> </mml:math> </inline-formula> x <inline-formula id="inf214"> <mml:math id="m217"> <mml:mrow> <mml:msub> <mml:mi>t</mml:mi> <mml:mrow> <mml:mi>f</mml:mi> <mml:mi>e</mml:mi> </mml:mrow> </mml:msub> <mml:mo> </mml:mo> </mml:mrow> </mml:math> </inline-formula>Grid

To justify the grid choice, a similar numerical analysis was carried out but using different grid spacings (2.5, 5, 7.5, 10, 25, 50, 75, and 100 ms). t c e and t f e values were varied between 2.5 and 505 ms, which led to 6 to 202 points for both t c e and t f e and grid sizes ranging from 36 to 40,804 total grid points. The boundary relationship between t c e and t f e was computed on each grid. RMSE on the testing set (15% points) as a function of the number of points along t c e is depicted in Figure A1. Noteworthy, as for grid spacings of 75 and 100 ms, using a 15% size for the testing set did not provide at least two points in such set. Therefore, two random points were forced to be attributed to the testing test (29 and 33% points in the testing set). As expected, RMSE decreased with decreasing grid spacing. Besides, it can be noticed that using a grid spacing of 10 ms did not seem to impact RMSE for the boundary relationship compared to the 7.5-ms grid spacing used before (RMSE ∼3.5 ms). However, the polynomial regression should also be performed on these different grids to observe any additional features.

Root-mean-square error (RMSE) as a function of the number of points along tce and ranging from 6 to 202. The red circle denotes RMSE corresponding to the boundary relationship computed in Section 3.1 (3.2 ms).

For this reason, a multivariate polynomial regression (polynomial order from 1 to 15) was performed on 85% of the points composing these different grids, after having discarded the points which were not within the corresponding boundary relationship. RMSE on the testing set (15% points) as a function of grid size is depicted for each polynomial order in Figure A1. It can be noticed that the eighth-order polynomial is the lowest order polynomial, leading to an RMSE smaller than 0.5 ms on the testing set. In addition, the smallest grid to obtain such an RMSE threshold is given by a grid using a spacing of 7.5 ms, that is, 4,624 grid points. As for the grid spacing of 10 ms, it requires a polynomial of order 10 to achieve the requested RMSE threshold, which is less convenient as it requires 21 extra coefficients than the eighth-order polynomial. Therefore, these previous statements justify the grid choice used to construct the multivariate polynomial model (Eq. 3).