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

Effective contact ( tce ) and flight ( tfe ) times, instead of ground contact ( tc ) and flight ( tf ) times, are usually collected outside the laboratory using inertial sensors. Unfortunately, tce and tfe cannot be related to tc and tf because the exact shape of vertical ground reaction force is unknown. However, using a sine wave approximation for vertical force, tce and tc as well as tfe and tf could be related. Indeed, under this approximation, a transcendental equation was obtained and solved numerically over a tce x tfe 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 tce values. Then, reconstructed tc (using the model) was compared to corresponding experimental ground truth. A systematic bias of 35 ms was depicted, demonstrating that ground truth tc values were larger than reconstructed ones. Nonetheless, error in the reconstruction of tc from tce 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 tc and tf in real time, which would allow coaches and practitioners to better evaluate running performance and to prevent running-related injuries.


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 ce ) and flight (t fe ) 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 , 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 ce and t fe . Thus, this temporal information is not lost. As for the effect of running speed, t ce and t c both decrease with increasing running speed, even though the decrease is much more important for t c than t ce (Cavagna et al., 2008;Da Rosa et al., 2019). Concerning t fe and t f , their trend with increasing running speed is not similar. Indeed, t fe 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 ce and t fe 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 ce and t fe .
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 mgπ(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 ce and t fe using t c t ce + 2t g and t f t fe − 2t 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 ) mg F z, max sin( πtg tc ), and using the definition of F z, max given before, the following equation is obtained: 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 ce and t fe . Ultimately, a mathematical modeling of t g over the FIGURE 1 | 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 ce ), and effective flight (t fe ) 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.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org November 2021 | Volume 9 | Article 687951 numerical t ce x t fe 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. 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 ce and t fe . As a direct experimental application, the proposed model was applied to experimentally measured t ce values. Then, the reconstructed t c values were compared to their corresponding experimental ground truth (gold standard).

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 ce and t fe . The zero of interest for a given t ce and t fe 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 ce and t fe 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 ce and t fe cannot be 0 and 2) t c belongs to the interval [100 ms, 400 ms] and t f belongs to the interval [0 ms, 250 ms] (Cavagna et al., 2008;Da Rosa et al., 2019;Lussiana et al., 2019), and to include any atypical t ce and t fe 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 t ce and t fe The numerical analysis showed that a linear boundary relationship is present between t ce and t fe (see Results Figure 2), that is, there is no solution for t g if t fe is higher than a certain percentage of t ce . This boundary relationship was computed by extracting the boundary points, that is, the smallest existing t fe values for every t ce 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 ce necessarily ensures a null t fe . 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).

Modeling a t g Surface as Function of t ce and t fe
The numerical analysis showed that t g could be described by a smoothly increasing surface when increasing t ce and t fe (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.

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 selfreported 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 zcomponent 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 ce and t fe 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 ce and t fe . Then, each t ce and t fe 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 ce and t fe (gold FIGURE 3 | Boundary relationship between t ce and t fe . 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%) . Frontiers in Bioengineering and Biotechnology | www.frontiersin.org November 2021 | Volume 9 | Article 687951 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 ce and t fe ) 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 ce and t fe .

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.

Numerical Analysis
The numerically calculated t g values over the t ce x t fe 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 100p2t g /(t ce + 2t g )].

Boundary Relationship Between t ce and t fe
The linear regression gave the model (Eq. 2): t fe 0.795 t ce . (2) Applying this model to the testing set led to anR 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.

Modeling a t g Surface as Function of t ce and t fe
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 ce , t fe ), Eq. 3]. The coefficients (α ij , where 0 ≤ i + j ≤ 8) of the multivariate polynomial model are given in Table 1.
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.

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).

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 ce and t fe . 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 n+2 2 , then C n+2 2 − 3 calculations are required to compute the polynomial features, that is, t i ce and t i fe , where 2 ≤ i ≤ n. In addition, C n+2 2 − 1 multiplications and C n+2 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 ce 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 ce . 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.

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 ce and t fe ) 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 FIGURE 7 | 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.
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.
FIGURE 6 | 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 .
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org November 2021 | Volume 9 | Article 687951 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.

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