Statistical Modeling of Lower Limb Kinetics During Deep Squat and Forward Lunge

Purpose: Modern statistics and higher computational power have opened novel possibilities to complex data analysis. While gait has been the utmost described motion in quantitative human motion analysis, descriptions of more challenging movements like the squat or lunge are currently lacking in the literature. The hip and knee joints are exposed to high forces and cause high morbidity and costs. Pre-surgical kinetic data acquisition on a patient-specific anatomy is also scarce in the literature. Studying the normal inter-patient kinetic variability may lead to other comparable studies to initiate more personalized therapies within the orthopedics.

Methods: Trials are performed by 50 healthy young males who were not overweight and approximately of the same age and activity level. Spatial marker trajectories and ground reaction force registrations are imported into the Anybody Modeling System based on subject-specific geometry and the state-of-the-art TLEM 2.0 dataset. Hip and knee joint reaction forces were obtained by a simulation with an inverse dynamics approach. With these forces, a statistical model that accounts for inter-subject variability was created. For this, we applied a principal component analysis in order to enable variance decomposition. This way, noise can be rejected and we still contemplate all waveform data, instead of using deduced spatiotemporal parameters like peak flexion or stride length as done in many gait analyses. In addition, this current paper is, to the authors’ knowledge, the first to investigate the generalization of a kinetic model data toward the population.

Results: Average knee reaction forces range up to 7.16 times body weight for the forwarded leg during lunge. Conversely, during squat, the load is evenly distributed. For both motions, a reliable and compact statistical model was created. In the lunge model, the first 12 modes accounts for 95.26% of inter-individual population variance. For the maximal-depth squat, this was 95.69% for the first 14 modes. Model accuracies will increase when including more principal components.

Conclusion: Our model design was proved to be compact, accurate, and reliable. For models aimed at populations covering descriptive studies, the sample size must be at least 50.

Introduction

In biomechanics, the safety and efficiency of novel surgical techniques as well as the development of biocompatible products ultimately rely on its capability of being tested on humans through clinical trials. The complete development chain of a new surgical technique or implant and their introduction into clinic practice is both time-consuming and economically demanding. Next to it, it is known that patient-specific surgery planning or implant design can improve the long-time outcome of an implant (Pietsch et al., 2013; Spencer-Gardner et al., 2016). This fact is due to the high anatomical variability between individuals and the different functional activities, which have a significant effect in the ratio of the force components on the lower limb between subjects (Kutzner et al., 2010) and on the functional alignment of the prosthetic components of a lower limb implant (Smoger, 2016; Spencer-Gardner et al., 2016). Within this context, methodologies such as statistical models of the human anatomy as well as kinematics or kinetics that account for the anatomical inter-variability of the population combined with biomechanical simulation studies can provide non-invasive pre-surgical clinical output.

Lower limb kinetics can be estimated based on musculoskeletal models and ground force plates using inverse dynamics (Carbone et al., 2012; Galloway et al., 2012; Vaitkus and Várady, 2015; Bagwell et al., 2016). These techniques do not often account for patient-specific variability as they use scaled generic models (Worsley et al., 2011; Vaitkus and Várady, 2015), while it was already widely shown that the geometry of the musculoskeletal models is very sensitive to muscle force predictions (Carbone et al., 2012). In addition, and to the authors knowledge, the available studies merely consider very limited population samples which may not be representative of the total variability of the lower limb anatomy. Lastly, the available literature lacks completeness as, to date, no study has considered a statistical model of the full lower limb, namely, on demanding tasks such as the deep squat and the forward lunge.

Hence, in order to create the foundations for the development and optimization of the design or the durability of orthopedic implants, it is mandatory to generate appropriate loading conditions that represent inter-patient variability across the population (Honari and Taylor, 2013; Bischoff et al., 2014). Patient-specific finite element analyses are the state-of-the-art technique to infer quantitative information on a specific design or performance of an arthroscopic implant (Shu et al., 2018). Taylor et al. (2012) found most studies to be focusing on variations on the morphological and bone properties rather than the consequences of variability because of loading. Furthermore, it has been proved that the application of single-representative models can be extended to account for variability by either parametrically or probabilistically varying the loading/boundary conditions. These approaches allow model generation which can significantly extend the coverage of the anatomical variability and ultimately create a powerful tool to assess the performance of medical devices (Taylor et al., 2012).

Recent developments in medical imaging significantly increased the accuracy of the three-dimension computational anatomical representation, enhancing the anatomical differences within a determined population (Almeida et al., 2016; Audenaert et al., 2019). Hence, combining the use of magnetic resonance imaging (MRI) with musculoskeletal models will provide us an insight on lower limb kinematics on patient-specific anatomies. The statistical analysis of kinematic time series by means of dimensionality reduction techniques such as principal component analysis (PCA) or independent component analysis is not novel per se (Daffertshofer et al., 2004; Galloway et al., 2012), but the inclusion of patient-specific anatomies is believed to more accurately represent inter-patient kinetic variability. Such approach, hereby presented, will allow for a large population of kinetic data to be generated without the time and the expense of collecting the motion capture data of hundreds of patients. Simultaneously, it will open the door to the generation of large simulated populations for use in clinical outcome simulation studies, injury biomechanics, musculoskeletal disease models, or implant design optimization (Henak et al., 2013; Zhang et al., 2016).

While the gait cycle has been the most researched activity in the current literature, it is not particularly demanding for the lower limb joints. For the purpose of implant wear testing, implant fixation, and joint stability analysis, there are other more challenging activities commonly performed in daily living that might be of particular interest (Hartmann et al., 2013). Clinical, experimental, and computational studies have clearly reported increased complication risk and wear rate under high contact stress conditions (Kang et al., 2008; O’Brien et al., 2015; de Ruiter et al., 2017).

In sum, the purpose of this study is to build statistical models of deep squatting and forward lunging for applications in pre-clinical testing of orthopedic implants and surgery in an asymptomatic adult population and ultimately to analyze and validate the inter-individual variations in lower limb kinetics.

Materials and Methods

Participants

Fifty-three asymptomatic volunteers participated in the study. In order to eliminate sex and race differences and reduce the potential influence of age and body mass index (BMI), only healthy Caucasian men who were not overweight and aged between 17 and 25 years are included. The admission requirement is practicing sports for at least 2 h a week. The subjects were asked to perform five times a smooth maximal-depth squat and a right forward lunge step with a predetermined frequency and fluency after a short training. In addition, the volunteers underwent full lower limbs MRI. An ethics committee (Ghent University Hospital, Belgium) approved these investigations (EC2014/0286). The characteristics of our study population are listed in Table 1. Because of missing data, there was no complete data acquisition for the squat among the three subjects.

TABLE 1

Table 1. Demographic and anthropometric characteristics of the study population.

In both examinations, 28 reflective markers are stuck on the skin on palpable anatomical landmarks. The application of skin markers to investigate kinetics is obvious but rather inaccurate. By contrast, using more accurate measurements with implants would raise ethical concerns.

Instrumentation

Our motion capture acquisition strategy was based on a similar study by Deluzio et al. (Deluzio and Astephen, 2007). Spatial marker trajectory data and the corresponding force registrations are imported into the Anybody Modeling System (AMS version 7.1.0, Anybody Technology, Aalborg, Denmark) (Damsgaard et al., 2006) as well as geometric data from a 3-Tesla MAGNETOM Trio-Tim System MRI device (Siemens AG, Erlangen, Germany).

Musculoskeletal Modeling

Motion capture musculoskeletal models were personalized with subject-specific bone geometry that was incorporated in a simulation model from the Twente Lower Extremity Model (TLEM 2.0) dataset (Carbone et al., 2015). An overview of the musculoskeletal model input is presented in Figure 1. In the simulation output, the forces are described in three fixed perpendicular planes.

FIGURE 1

Figure 1. Overview of data input for the motion capture musculoskeletal simulation model. (A) Motion is performed when standing on two force plates. Motion capture data synchronized with ground reaction forces are exported as .c3d file. (B) Twenty-eight reflective markers are placed on anatomical bony landmarks. A MRI scan of the full lower limb is performed. Segmentation of pelvis, thigh, and shank with corresponding positions of marker landscape. (C) Motion capture squat model. Anybody squat (D) and lunge (E) model.

Data Processing

The output data from musculoskeletal models are numerous, multivariate, and multidimensional (Deluzio and Astephen, 2007; Lai et al., 2009). In contrast to some gait studies that modeled kinematic and kinetic data together, we used only kinetic data (Deluzio and Astephen, 2007; Reid et al., 2010; Galloway et al., 2012). We think that integrating linear quantities (forces) and rotation quantities (angles) is like comparing apples and oranges. On top of that, the kinetic data in Anybody is generated by an inverse dynamics approach starting from the kinematic data.

The beginning and end frames of all motion lab recordings are not useful due to irrelevant transients. Analogously, the peak evolution will vary from the center of the recorded data. Hence, data alignment and trimming are essential prior to incorporating the subjects’ motion recordings into a statistical model. These operations are executed using standard implementations in MATLAB (MathWorks, Natick, MA, United States).

The frame recorded with the peak knee flexion angle is defined as 50% progress of the motion. Trimming is based on knee flexion. For the lunge, the best is to consider only the closed chain part. As such, recordings where the right foot is not on the right force plate are left aside. Several arbitrary ways to execute an open-chain motion could be an important source of noise η. Noise is defined as artifacts when processing the input data to the output data (Lai et al., 2009). On top of data, we used only information from the leg that was the most loaded. So, in contrast to the squat data, a lot of waveform data are not used for the lunge.

Interpolation is performed to ensure that the measurements are running synchronized in real time. All trimmed measurements are subdivided into 0–50–100 proceedings, corresponding to the onset, the middle, and the finish of motion, respectively. Each set of kinetic data is arranged in a feature vector and concatenated into a training matrix. The training data matrix X contains observations in the rows and subjects in the columns as described in Eqs. [1] to [3].

$$X=[x_1,x_2,x_3,\mathrm{\dots },x_i,\mathrm{\dots },x_{p-1},x_p](1)$$

An observation expresses several dynamic parameters on a certain progress of the aligned lunge or squat motion from 0 to 100%. For each participant i, the kinetic model input data are taken from the musculoskeletal model output. The kinetic variables are implemented into a subject vector x_i for the ith subject (out of p).p represents the number of training samples, being 53 for the lunge and 50 for the squat.

$$x_i={[HJR{F}_x,HJR{F}_y,HJR{F}_z,KJR{F}_x,KJR{F}_y,KJR{F}_z,]}^T(2)$$

The input matrices JRF_axis consist all of 101 observations o.

$$JR{F}_{\text{axis}}=\left[\begin{array}{c}\hfill \begin{array}{c}\hfill \text{measurement}\mathrm{\hspace{0.33em}1}\hfill \\ \hfill \text{measurement}\mathrm{\hspace{0.33em}2}\hfill \\ \hfill \mathrm{\dots }\hfill \end{array}\hfill \\ \hfill \mathrm{\dots }\hfill \\ \hfill \begin{array}{c}\hfill \text{measurement}51(\mathrm{at}\mathrm{maximal}\mathrm{right}\mathrm{knee}\mathrm{flexion})\hfill \\ \hfill \mathrm{\dots }\hfill \\ \hfill \text{measurement}101\hfill \end{array}\hfill \end{array}\right](3)$$

D serves as a diagonal matrix with row-wise standard deviations d_o for each observation o. The total number of observations is the multiplication of the number of dynamic variables and aligned time instances.

$$D=\left[\begin{array}{cc}\hfill \begin{array}{cc}\hfill d_1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill d_2\hfill \end{array}\hfill & \hfill \begin{array}{cc}\hfill 0\hfill & \hfill \mathrm{\dots }\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{array}\hfill \\ \hfill \begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \\ \hfill \mathrm{\dots }\hfill & \hfill 0\hfill \end{array}\hfill & \hfill \begin{array}{cc}\hfill \mathrm{\dots }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill d_{606}\hfill \end{array}\hfill \end{array}\right](4)$$

After normalization by row-wise standard deviation in [4] and [5] as well as mean centering in [6], a residual matrix R is created. R comprises the entry data for the model M as a measure of dispersion.

$$\tilde{X}=\left[\begin{array}{ccc}\hfill \mathrm{\dots }\hfill & \hfill \begin{array}{c}\hfill \mathrm{\dots }\hfill \\ \hfill {\tilde{x}}_i\hfill \\ \hfill \mathrm{\dots }\hfill \end{array}\hfill & \hfill \mathrm{\dots }\hfill \end{array}\right]=D^{-1}X(5)$$

$$R=\left[\begin{array}{ccc}\hfill \mathrm{\dots }\hfill & \hfill \begin{array}{c}\hfill \mathrm{\dots }\hfill \\ \hfill {\widehat{x}}_i\hfill \\ \hfill \mathrm{\dots }\hfill \end{array}\hfill & \hfill \mathrm{\dots }\hfill \end{array}\right]with{\widehat{x}}_i={\tilde{x}}_i-\overline{x}(6)$$

PCA is a powerful dimensionality reduction technique developed by Karl Pearson. It is not a method to investigate the center size of the data but the common variability. PCA is mathematically defined as an orthogonal linear transformation. PCA transforms the data; as such, most of the variance of the data will come to lie in the first components. This allows us to create statistical models. Altogether PCA is a reliable tool in capturing the salient features of waveform data (Robbins et al., 2013; Jolliffe and Cadima, 2016).

Using this for a statistical model, it enables to generate population data from a small set of clinical data. The kinetic model should represent waveform data as a linear combination of vectors, representing the primary modes of variation in experimental data (Jolliffe et al., 2002; Saliba et al., 2018). Eigenvalues and eigenvectors have been created by singular value decomposition.

$$R=UxLx{A}^T(7)$$

In Eq. [7], U and A are the left and right singular vectors, so U^T.U = I and A^T.A = I because of orthogonality. I refers to the unity matrix. L is a diagonal matrix that contains the square roots of the eigenvalues $\sqrt{{\mathrm{\lambda }}_k}$ belonging to R^T x R, as k{1,2,3,…,p}. A^T contains the eigenvectors of R^TxR, whereas U has the scaled versions of the principal component (PC) scores u_ok. Here, o corresponds to the observation and k to the PC number. The PC scores are mentioned in Eq. [8].

$$z_{ok}=u_{ok}\sqrt{{\mathrm{\lambda }}_k}(8)$$

The PC scores from a single waveform quantify the contribution of each feature. The variance of the scores for the kth eigenvalue of R^T x R amounts to $\frac{{\mathrm{\lambda }}_k}{n-1}$, as λ_k represents the variance of the kth PC, whereas n is the number of observations.

The cumulative variance of each mode k is defined as

$$\text{Compactness}\left(M\right)=\frac{1}{p-1}\sum _{m=1}^k{\mathrm{\lambda }}_m(9)$$

The PC weight matrix W in Eq. [10] involves the correlation coefficients between components and test subject data.

$$W=Lx{A}^T(10)$$

A set of patient data can be approximately reconstructed by using t selected PCs in Eq. [11].

$${\widehat{x}}_{ij}=D\sqrt{p-1}\sum _{i=1}^t{u}_{oi}\sqrt{l_i}{a}_{ki}+\overline{x}(11)$$

As mentioned before, Drepresents the diagonal matrix of row-wise standard deviations and p stands for the subject count. $\sqrt{l_i}$ is the ith diagonal element of L, also from the singular value decomposition (Jolliffe et al., 2002; Galloway et al., 2012; Jolliffe and Cadima, 2016).

Model Validation

Validation is defined as the process of ensuring that the dimensionality-reduced PCA model accurately represents real-world kinetics. Probably the most important problem arising with this process is the choice of the optimal number of the principal components to be retained. PCA projects the input data from a high dimensional space into a subspace of lower dimension, which can then further be divided into two separate subspaces: the kinetic data subspace, preserving the essence of the original kinetic data as lossless as possible, and the noise subspace, corresponding to the remaining tail of principal components associated with the smallest eigenvalues. Given the complexity of the problem of optimally defining the threshold between signal and noise principal components, the literature on the topic is overwhelming and beyond the scope of this work. Reliable results in distinguishing components that express meaningful correlations among variables as opposed to trivial components, explaining noise, have been provided using the Monte Carlo permutation test (Peres-Neto et al., 2005). The principal components were tested for representing valid correlations as opposed to residual error using the following two criteria: rank of roots and equality of roots (Jackson, 2005; Vasco, 2012).

Further, four quantitative model parameters are investigated. “Goodness” measures are chosen according to the statistical shape modeling study of Styner et al. (2003) in which there is also a PCA dimensionality reduction algorithm. This study is, to the authors’ best knowledge, the first to provide such an approach, implemented for a kinetic model.

Model Accuracy

$$\text{Accuracy}\left(M\right)=\frac{1}{p}\sum _{i=1}^p{||{\widehat{x}}_i(M)-x_i||}^2(12)$$

The first validation test that analyzes relevant information is retained by the model or otherwise states how well the original data can be reconstructed from the model given the number of principal components retained. Here, the root-mean-square error (RMSE) is computed in Eq. [12] as the average absolute difference between the original training data and the reconstructed data for models with 80, 90, 95, and 98% variance of the original data.

Model Compactness

The model will be compact enough if it can describe the variance in kinetic measurements with a minimal number of modes. Eq. [9] is used to describe the compactness with the cumulative variance for a certain number of modes.

Model Generalization

$$\text{Generalization}\left(M\right)=\frac{1}{T_g}\sum _{i=1}^{T_g}{||{\widehat{x}}_i(M_i^{\prime })-x_i||}^2(13)$$

The model generalization quantifies the ability of models to represent new instances. The generalization ability is evaluated by performing a series of leave-one-out tests on the training data. The question here is: how many training samples are necessary to approach the population precisely? The generalization ability is therefore a means for post hoc sample size evaluation. If having enough training samples, we expect the model to be able to describe unseen data quite accurately (Wang and Shi, 2017). The generalization value can be interpreted as the median out-of-sample accuracy value.

The generalization evolution gives the RMSE between the excluded subject data and the best-matched 95% variance model M’ values of randomly selected training data by ascending number of training samples in the model M’. The higher the T_g test value in Eq. [13], the higher the precision of the median generalization value. Here the number of models created for each number of training samples amounts to T_g = 10,000.

Model Specificity

$$\text{Specificity}\left(M\right)=\frac{1}{T_s}\sum _{i=1}^{T_s}{||{\widehat{x}}_i^{\prime }(M)-x_i||}^2(14)$$

A population model is able to generate new data. The model specificity measures the soundness of new instances randomly generated by the developed model M. Models with 80, 90, 95, and 98% of variance are tested. ${\widehat{x}}_i^{\prime }(M)$ refers to a randomly generated subject.

We assume that the PCs of the model are normally distributed (Jolliffe et al., 2002; Jackson, 2005; Galloway et al., 2012). The specificity estimator is defined in Eq. [14]. For each observation o, an imaginary subject i is defined by choosing random normal distributed values n ∈ N(0,1) for each mode m in the model M as in Eq. [15].

$${\widehat{x}}_{i,o}^{\prime }\left(M\right)={\overline{x}}_o+d_o\sum _{m=1}^{p}n.\sqrt{{\mathrm{\lambda }}_m}.z_m=\overline{x}+d_o\sum _{m=1}^{p}n.u_m(15)$$

The RMSE is defined as the error between the virtually subject data and the most similar sample in the training dataset. The specificity value can be interpreted as the median approximation error of T_s generated subjects. The higher the T_s test value, the higher the precision of the specificity. Here the number of models created is set to T_s = 1,000,000.

Results

The average hip and knee peak flexion angles are, respectively 95° and 104° for the lunge and 107° and 112° for the squat motion, respectively. The average peak hip joint reaction force (HJRF) amounts to 3.08 times body weight (BW) for the maximum-depth squat and 4.76 BW for lunging. The means for peak knee joint reaction force (KJRF) are still higher: 4.52 BW for squat motion and 7.16 BW for the lunge. The trimmed original waveform data from HJRF and KJRF of our musculoskeletal model are represented by gray curves in Figures 3–8.

A statistical model of kinetic output data from the AMS was made for deep squatting and another one for lunging. Figure 2 displays the cumulative variance of modes in the statistical model. The variances of the first three modes in the squat model are illustrated in Figures 3–5. Together they represent 59.73% of the population variance. For the lunge, the first three modes accounts for 66.40% of population variance. More about these modes are detailed in Figures 6–8. In Table 2, the in-sample model accuracy and the specificity median RMSE are described for each model. The boxplots in Figure 9 illustrate the in-sample squat model accuracy for the ascending number of the components included. There is a boxplot for every variable in the model. The median and the interquartile range of the RMSE, when compared to the initial data, decrease as more PCs are included in the reconstructed data, as expected. The model accuracies from the lunge model are quite similar but around one IQR RMSE higher. In contrast, the lunge model is a bit more compact. For calculating the out-of-sample accuracy based on leave-one-out tests, we based on the lunge data because it has the most test data. The results of the model compactness and the statistical findings of the permutation testing related to the number of the principal components used are demonstrated in Table 3.

FIGURE 2

Figure 2. Scree plot with the cumulative variance of the modes (or principal components) in the lunge (orange) and squat (purple) kinetic model.

FIGURE 3

Figure 3. Relation between the kinetic waveform simulation output and the squat progress for each individual sample in gray. Mean values of the measurements in green ±2 standard deviations of the first mode in red and blue. The first mode accounts for 33.80% of the inter-subject population variance. Note the different y axis calibrations.

FIGURE 4

Figure 4. Mean values of joint reaction forces during deep squatting in green ±2 standard deviations of the second mode in red and blue. The second mode accounts for 14.05% of the inter-subject population variance.

FIGURE 5

Figure 5. Mean values of joint reaction forces during deep squatting in green ±2 standard deviations of the third mode in red and blue. The third mode accounts for 11.88% of the inter-subject population variance.

FIGURE 6

Figure 6. Mean values of joint reaction forces during lunging in green ±2 standard deviations of the first mode in red and blue. The first mode accounts for 40.87% of the inter-subject population variance.

FIGURE 7

Figure 7. Mean values of joint reaction forces during lunging in green ±2 standard deviations of the second mode in red and blue. The second mode accounts for 15.07% of the inter-subject population variance.

FIGURE 8

Figure 8. Mean values of joint reaction forces during lunging in green ±2 standard deviations of the third mode in red and blue. The third mode accounts for 10.46% of the inter-subject population variance.

FIGURE 9

Figure 9. RMSE for the original squat training data versus reconstructed squat data with an increasing number of principal components on the x axis.

TABLE 2

Table 2. Validation analyses of the squat and lunge statistical kinetic model. We consider squat and lunge models that capture 80, 90, 95, and 98% of inter-individual population variance.

TABLE 3

Table 3. Choosing the optimal amount of principal components for the squat kinetic datasets.

Regarding Figure 10, for each training data input amount going from 4 to 52, 10,000 models were created, including 95% population variance, to reconstruct an excluded subject. Out-of-sample accuracy RMSE from the reconstructed data versus the original excluded data are given on the y axis in box-and-whisker diagrams. The boxplots are log–log-scaled in order to visualize the downward trend of the out-of-sample accuracy. Also plotted is a horizontal line of the in-sample model accuracy of our 95% model. The out-of-sample accuracies are less than 0.1 BW, except for the KJRF in the transverse plane. From 50 test subjects up, the out-of-sample accuracies are clearly stagnating for the HJRF in the frontal and the transverse plane as well as for the KJRF in the sagittal plane.

FIGURE 10

Figure 10. Accuracy evolution of kinetic lunge data with log–log scaling (boxplot with root-mean-square error of the reconstructed data with 95% variance versus the original training data) for different levels of prior knowledge expressed as amounts of training data in a kinetic model. The green horizontal line indicates the in-sample target accuracy.

Discussion

The validation analysis confirmed that our models have a high degree of compactness and accuracy. Many types of noise are in the higher components. The PCA technique has adequately allowed rejection of the error variance from the model. The meaningful variance is obviously divided over the first 12 or 14 components. This multidimensionality describes the silent features in the data and, eventually, they could be linked to the varying characteristics of the study population. A common source of meaningless variance originates from data alignment. It is impossible to avoid this because we do not want to introduce supplementary noise in the data by aligning them more stringently. Since all subjects have a BMI lower than 25, skin shift errors during movements are limited (Cappozzo et al., 1996).

According to the lunge, the model only describes the closed-chain part of motion for two reasons. First, femoroacetabular impingement and joint reaction forces are more pronounced at higher flexion (Audenaert et al., 2012). Secondly, while creating a model from the onset of the lunge back to the original position, the model would be no longer compact enough because there is too much degree of freedom when moving a leg in the air.

The dominant mode is supposed to describe the overall variance (Jolliffe et al., 2002), as is clearly apparent in the lunge model. In the squat model, the overall variance is limited for the HJRF in the frontal and the sagittal plane. This is due to low hip flexion and rotation moments because the center of gravity will lie almost perfectly between the hip joints (and not the knee joints during squat), in contrast to the lunge case. For this statement, we based on Schwab et al. (2006). They found that, for young adults, the femoral head position appears to be a reliable indicator for the gravity line in the sagittal plane during stance (Schwab et al., 2006).

The second mode of the squat model indicates that a high HJRF component in the transverse plane results in high KJRF components in the frontal and sagittal plane in order to counterbalance the downward force at the hip. The third mode correlates the depth of squatting with the joint reaction force components in the frontal plane. For the lunge, the association of the frontal joint reaction force components is mainly summarized in the second mode. Finally, according to our interpretation, the third mode of the lunge model may take alignment errors into account.

The RMSE for model accuracy are far below 0.05 BW, as opposed to similar studies. The specificity was almost equal for models with 80, 90, 95, and 98% of variance. It questions the relevance of taking the model specificity into account in this setting. According to the generalization evolution, we could conclude that, minimally, 50 samples are enough to provide reliable models at 0.1 BW precision for both squat and lunge motion. Nevertheless, we recommend exceeding this threshold number because the in-sample accuracy is still lower, especially for the squat. Note that gender, age group, BMI group, and race differences are not included here. Therefore, it is very likely that, in more heterogeneous populations like the elderly, 50 samples will be too low to ensure reliable models.

Unfortunately, electromyography data are not collected during this study. This could give information about muscle activation and muscle strength. Motor unit action potentials could be registered non-invasively by using surface electromyography. It has been stated several times that the muscle activation patterns depend on several aspects like training level and osteoarthritis (Benedetti et al., 2003; Knoop et al., 2012; Mei et al., 2017). The integration of electromyography and kinetic data could help to declare aberrant kinetic patterns.

By applying correlation matrix PCA to obtain uncorrelated maximum-variance linear combinations and given that there is only kinetic data with limited scaling differences, some more PCs are required to account for the same amount of covariance compared to classical covariance matrix PCA (Jolliffe et al., 2002; Jolliffe and Cadima, 2016). This makes the selection of PCs in the kinetic data subspace even more crucial to ensure model validation properties like accuracy, compactness, generalization, and specificity, which is the major drawback of PCA (Jolliffe et al., 2002; Peres-Neto et al., 2005; Vasco, 2012). To handle this, there are numerous methods described in the literature, but there is no consensus yet. We objectified our selection strategy based on eigenvalues by considering the validation measures for different cutoffs. On top of that, for the model generalization and specificity abilities, we assume multivariate normal distribution which is seldom true (Vasco, 2012).

The most important limitation of the present work, however, relates to the population under investigation, namely, young male, Belgian adolescents and the unknown extent of which findings can be extrapolated to other populations. Nevertheless, in general terms, we expect our results to be representative by extension for a Western European population.

Conclusion

We created two models that describe kinetics from both hip and knee joint, contrary to the limited number of studies available with PCA analyses of waveform data considering the knee only (Deluzio and Astephen, 2007; Reid et al., 2010; Galloway et al., 2012). Since all muscles from the knee, except from the M. popliteus (Paulsen and Waschke, 2011), are biarticular and the body should be seen as a whole, a model with the HJRF as well as KJRF is preferable. We proved that such a model for 95% of population variance was compact and very accurate (<0.015 BW). To describe the population at <0.1 BW precision, our small sample size was still sufficient. Using t-tests to investigate differences in PC scores in comparing studies will enable the creation of personalized hip and tibial implants with specific weight-bearing properties, resulting in prolonged longevity. In our opinion, this is a very important feature since total knee arthroplasty and total hip arthroplasty are increasingly utilized to treat more physically active patients.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by the Commission for Medical Ethics, UZ Gent. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

JV was financially supported by Ph.D. grant 11V2215N from the Research Foundation Flanders. EA was financially supported by Senior Clinical Fellowship from the Research Foundation Flanders.

Conflict of Interest

PG is an employee of AnyBody Technology. No financial benefits have been received or will be received from any commercial party related directly or indirectly to the subject of this article.

The remaining 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.

Acknowledgments

The authors would like to thank Ashwin Schouten for his contribution to the musculoskeletal modeling in the AMS.

References

Almeida, D. F., Ruben, R. B., Folgado, J., Fernandes, P. R., Audenaert, E., Verhegghe, B., et al. (2016). Fully automatic segmentation of femurs with medullary canal definition in high and in low resolution CT scans. Med. Eng. Phys. 38, 1474–1480. doi: 10.1016/j.medengphy.2016.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Audenaert, E. A., Pattyn, C., Steenackers, G., De Roeck, J., Vandermeulen, D., and Claes, P. (2019). Statistical shape modelling of skeletal anatomy for sex discrimination: their training size, sexual dimorphism and asymmetry. Front. Bioeng. Biotechnol. 7:302. doi: 10.3389/fbioe.2019.00302

PubMed Abstract | CrossRef Full Text | Google Scholar

Audenaert, E. A., Peeters, I., Vigneron, L., Baelde, N., and Pattyn, C. (2012). Hip morphological characteristics and range of internal rotation in femoroacetabular impingement. Am. J. Sports Med. 40, 1329–1336. doi: 10.1177/0363546512441328

PubMed Abstract | CrossRef Full Text | Google Scholar

Bagwell, J. J., Snibbe, J., Gerhardt, M., and Powers, C. M. (2016). Hip kinematics and kinetics in persons with and without cam femoroacetabular impingement during a deep squat task. Clin. Biomech. 31, 87–92. doi: 10.1016/j.clinbiomech.2015.09.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Benedetti, M. G., Catani, F., Bilotta, T. W., Marcacci, M., Mariani, E., and Giannini, S. (2003). Muscle activation pattern and gait biomechanics after total knee replacement. Clin. Biomech. 18, 871–876. doi: 10.1016/s0268-0033(03)00146-3

CrossRef Full Text | Google Scholar

Bischoff, J. E., Dai, Y., Goodlett, C., Davis, B., and Bandi, M. (2014). Incorporating population-level variability in orthopedic biomechanical analysis: a review. J. Biomech. Eng. 136:021004. doi: 10.1115/1.4026258

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappozzo, A., Catani, F., Leardini, A., Benedetti, M. G., and Croce, U. D. (1996). Position and orientation in space of bones during movement: experimental artefacts. Clin. Biomech. 11, 90–100. doi: 10.1016/0268-0033(95)00046-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbone, V., Fluit, R., Pellikaan, P., van der Krogt, M. M., Janssen, D., Damsgaard, M., et al. (2015). TLEM 2.0 – A comprehensive musculoskeletal geometry dataset for subject-specific modeling of lower extremity. J. Biomech. 48, 734–741. doi: 10.1016/j.jbiomech.2014.12.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbone, V., van der Krogt, M. M., Koopman, H. F., and Verdonschot, N. (2012). Sensitivity of subject-specific models to errors in musculo-skeletal geometry. J. Biomech. 45, 2476–2480. doi: 10.1016/j.jbiomech.2012.06.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Daffertshofer, A., Lamoth, C. J., Meijer, O. G., and Beek, P. J. (2004). PCA in studying coordination and variability: a tutorial. Clin. Biomech. 19, 415–428. doi: 10.1016/j.clinbiomech.2004.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Damsgaard, M., Rasmussen, J., Christensen, S. T., Surma, E., and De Zee, M. (2006). Analysis of musculoskeletal systems in the anybody modeling system. Simul. Model. Pract. Theory 14, 1100–1111. doi: 10.1016/j.simpat.2006.09.001

CrossRef Full Text | Google Scholar

de Ruiter, L., Janssen, D., Briscoe, A., and Verdonschot, N. (2017). The mechanical response of a polyetheretherketone femoral knee implant under a deep squatting loading condition. Proc. Instit. Mech. Eng. H J. Eng. Med. 231, 1204–1212. doi: 10.1177/0954411917738805

PubMed Abstract | CrossRef Full Text | Google Scholar

Deluzio, K. J., and Astephen, J. L. (2007). Biomechanical features of gait waveform data associated with knee osteoarthritis: an application of principal component analysis. Gait Posture 25, 86–93. doi: 10.1016/j.gaitpost.2006.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Galloway, F., Worsley, P., Stokes, M., Nair, P., and Taylor, M. (2012). Development of a statistical model of knee kinetics for applications in pre-clinical testing. J. Biomech. 45, 191–195. doi: 10.1016/j.jbiomech.2011.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghaffari, A., Davis, I., Storey, T., and Moser, M. (2018). Current concepts of femoroacetabular impingement. Radiol. Clin. North Am. 56, 965–982.

Google Scholar

Hartmann, H., Wirth, K., and Klusemann, M. (2013). Analysis of the load on the knee joint and vertebral column with changes in squatting depth and weight load. Sports Med. 43, 993–1008. doi: 10.1007/s40279-013-0073-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Henak, C. R., Anderson, A. E., and Weiss, J. A. (2013). Subject-specific analysis of joint contact mechanics: application to the study of osteoarthritis and surgical planning. J. Biomech. Eng. 135:021003. doi: 10.1115/1.4023386

PubMed Abstract | CrossRef Full Text | Google Scholar

Honari, G., and Taylor, J. S. (2013). Commentary on Crawford, et al., The role of patch testing in the evaluation of orthopedic implant-related adverse effects: current evidence does not support broad use. Dermatitis 24, 108–111. doi: 10.1097/der.0b013e3182979722

PubMed Abstract | CrossRef Full Text | Google Scholar

Jackson, J. E. (2005). A user’s Guide to Principal Components. Hoboken, NJ: John Wiley & Sons.

Google Scholar

Jolliffe, I. T., Bickel, P., Diggle, P., Fienberg, S., Krickeberg, K., Olkin, I., et al. (2002). Principal Component Analysis, 2nd Edn. Cham: Springer.

Google Scholar

Jolliffe, I. T., and Cadima, J. (2016). Principal component analysis: a review and recent developments. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 374:20150202. doi: 10.1098/rsta.2015.0202

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, L., Galvin, A. L., Brown, T. D., Jin, Z., and Fisher, J. (2008). Quantification of the effect of cross-shear on the wear of conventional and highly cross-linked UHMWPE. J. Biomech. 41, 340–346. doi: 10.1016/j.jbiomech.2007.09.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Knoop, J., Dekker, J., Klein, J. P., van der Leeden, M., van der Esch, M., Reiding, D., et al. (2012). Biomechanical factors and physical examination findings in osteoarthritis of the knee: associations with tissue abnormalities assessed by conventional radiography and high-resolution 3.0 Tesla magnetic resonance imaging. Arthritis Res. Ther. 14:R212. doi: 10.1186/ar4050

PubMed Abstract | CrossRef Full Text | Google Scholar

Kutzner, I., Heinlein, B., Graichen, F., Bender, A., Rohlmann, A., Halder, A., et al. (2010). Loading of the knee joint during activities of daily living measured in vivo in five subjects. J. Biomech. 43, 2164–2173. doi: 10.1016/j.jbiomech.2010.03.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Lai, D. T., Begg, R. K., and Palaniswami, M. (2009). Computational intelligence in gait research: a perspective on current applications and future challenges. IEEE Trans. Inf. Technol. Biomed. 13, 687–702. doi: 10.1109/TITB.2009.2022913

PubMed Abstract | CrossRef Full Text | Google Scholar

Mei, Q., Gu, Y., Fu, F., and Fernandez, J. (2017). A biomechanical investigation of right-forward lunging step among badminton players. J. Sports Sci. 35, 457–462. doi: 10.1080/02640414.2016.1172723

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Brien, S. T., Luo, Y., and Brandt, J.-M. (2015). In-vitro and in-silico investigations on the influence of contact pressure on cross-linked polyethylene wear in total knee replacements. Wear 332, 687–693. doi: 10.1016/j.wear.2015.02.048

CrossRef Full Text | Google Scholar

Paulsen, F., and Waschke, J. (2011). Sobotta Atlas of Human Anatomy, 15th Edn. Amsterdam: Elsevier.

Google Scholar

Peres-Neto, P. R., Jackson, D. A., and Somers, K. M. (2005). How many principal components? Stopping rules for determining the number of non-trivial axes revisited. Comput. Stat. Data Anal. 49, 974–997. doi: 10.1016/j.csda.2004.06.015

CrossRef Full Text | Google Scholar

Pietsch, M., Djahani, O., Hochegger, M., Plattner, F., and Hofmann, S. (2013). Patient-specific total knee arthroplasty: the importance of planning by the surgeon. Knee Surg. Sports Traumatol. Arthrosc. 21, 2220–2226. doi: 10.1007/s00167-013-2624-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Reid, S. M., Graham, R. B., and Costigan, P. A. (2010). Differentiation of young and older adult stair climbing gait using principal component analysis. Gait Posture 31, 197–203. doi: 10.1016/j.gaitpost.2009.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Robbins, S. M., Astephen Wilson, J. L., Rutherford, D. J., and Hubley-Kozey, C. L. (2013). Reliability of principal components and discrete parameters of knee angle and moment gait waveforms in individuals with moderate knee osteoarthritis. Gait Posture 38, 421–427. doi: 10.1016/j.gaitpost.2013.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Saliba, C. M., Clouthier, A. L., Brandon, S. C. E., Rainbow, M. J., and Deluzio, K. J. (2018). Prediction of knee joint contact forces from external measures using principal component prediction and reconstruction. J. Appl. Biomech. 34, 419–423. doi: 10.1123/jab.2017-0262

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwab, F., Lafage, V., Boyce, R., Skalli, W., and Farcy, J. P. (2006). Gravity line analysis in adult volunteers: age-related correlation with spinal parameters, pelvic parameters, and foot position. Spine 31, E959–E967.

PubMed Abstract | Google Scholar

Shu, L., Yamamoto, K., Yao, J., Saraswat, P., Liu, Y., Mitsuishi, M., et al. (2018). A subject-specific finite element musculoskeletal framework for mechanics analysis of a total knee replacement. J. Biomech. 77, 146–154. doi: 10.1016/j.jbiomech.2018.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Smoger, L. M. (2016). Statistical Modeling to Investigate Anatomy and Function of the Knee. Denver: Thesis.

Google Scholar

Spencer-Gardner, L., Pierrepont, J., Topham, M., Baré, J., McMahon, S., and Shimmin, A. (2016). Patient-specific instrumentation improves the accuracy of acetabular component placement in total hip arthroplasty. Bone Joint J. 98, 1342–1346. doi: 10.1302/0301-620x.98b10.37808

PubMed Abstract | CrossRef Full Text | Google Scholar

Styner, M. A., Rajamani, K. T., Nolte, L. P., Zsemlye, G., Szekely, G., Taylor, C. J., et al. (2003). Evaluation of 3D correspondence methods for model building. Inf. Process. Med. Imaging 18, 63–75. doi: 10.1007/978-3-540-45087-0_6

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, M., Barrett, D. S., and Deffenbaugh, D. (2012). Influence of loading and activity on the primary stability of cementless tibial trays. J. Orthopaedic Res. 30, 1362–1368. doi: 10.1002/jor.22056

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaitkus, M., and Várady, T. (eds) (2015). “A general framework for constrained mesh parameterization,” in Proceedings of the 31st Spring Conference on Computer Graphics, (New York, NY: ACM).

Google Scholar

Vasco, M. (2012). Permutation tests to estimate significances on Principal Components Analysis. Comput. Ecol. Softw. 2:103.

Google Scholar

Wang, J., and Shi, C. (2017). Automatic construction of statistical shape models using deformable simplex meshes with vector field convolution energy. Biomed. Eng. Online 16:49. doi: 10.1186/s12938-017-0340-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Waxman, A. (2004). WHO Global Strategy on Diet, Physical Activity and Health. Geneva: World Health Organization.

Google Scholar

Worsley, P., Stokes, M., and Taylor, M. (2011). Predicted knee kinematics and kinetics during functional activities using motion capture and musculoskeletal modelling in healthy older people. Gait Posture 33, 268–273. doi: 10.1016/j.gaitpost.2010.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Fernandez, J., Hislop-Jambrich, J., and Besier, T. F. (2016). Lower limb estimation from sparse landmarks using an articulated shape model. J. Biomech. 49, 3875–3881. doi: 10.1016/j.jbiomech.2016.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar