We developed a three dimensional musculoskeletal model of an alpine skier with two skis and 53 degrees of freedom (19 for the skier and 17 for each ski) to simulate turning maneuvers in alpine skiing (Figure 1). The skeletal model of the skier consisted of 20 rigid segments and was derived from the full-body OpenSim model of Hamner et al. (2010). At each lower extremity the subtalar and mtp joints were locked because of the skier’s ski boot, which allows only plantarflexion and dorsiflexion at the ankle joint. The restraining effect of the ski boot was represented by a passive moment at the ankle joint incorporating the non-linear relation between the boot-induced moment and the ankle joint angle (Eberle et al., 2017). To increase computational speed the position of the arms of the skier was locked in a typical position and the mass of the ski poles was neglected. In total, the skier model had 19 degrees of freedom (6 between pelvis and ground; 3, 1, and 1 at each hip, knee and ankle, respectively; 3 at the lumbar joint between trunk and pelvis).

Each ski was discretized into 18 rigid segments (7 rear segments, 1 center segment and 10 front segments) connected by revolute joints. Mass and inertia properties of the ski segments were derived form measurement data of a competitive giant slalom ski. The length of the ski was 2.02 m with a sidecut radius of 32 m and a mass of 2.1 kg. The center segment was firmly affixed to the foot-ski boot segment of the skier model. Rotational spring-damper elements were attached to the revolute joints to incorporate stiffness and damping properties of the skis (Heinrich et al., 2014). Stiffness and damping parameters were derived from laboratory measurements of bending deflection and bending vibration (Mössner et al., 2014). Torsional twist of the skis was neglected because it was shown to be low during turning maneuvers (Yoneyama et al., 2008).

The motion of the skier was actuated by 94 muscles (43 per leg and eight actuating the lumbar joint). The muscle model was based on the OpenSim model of Catelli et al. (2019), since deep squatting and high hip flexion are often encountered during turning maneuvers. In the OpenSim model of Catelli et al. (2019), 80 muscles were includes (40 per leg) actuating the joints of both legs only. To actuate the lumbar joint additionally, we added eight muscles (2 × erector spinae, rectus abdominus, external obliques, internal obliques) to the skier model (Harris et al., 2017). Furthermore, we added three hip muscles (gemelli, pectineus, quadratus femoris) to each leg to increase hip muscle strength based on the musculoskeletal model presented in Harris et al. (2017) focusing on hip musculature.

Muscle activation dynamics was assumed as a first-order process (He, 1991) and models the change in muscle activation as a function of the current active state and the neural excitation of the muscle. The corresponding time constants for muscle activation and deactivation were set to 10 and 40 ms, respectively (Zajac, 1989; Nitschke et al., 2020). Muscle contraction dynamics was modeled in analogy to Nitschke et al. (2020) assuming a three-element Hill-type muscle model, which incorporates the force-length-velocity characteristics of the muscle. Contraction dynamics was formulated in implicit form, which was shown to result in better convergence and increased computational speed (van den Bogert et al., 2011; De Groote et al., 2016).

The dynamics of the whole musculoskeletal skier model was given by the multibody dynamics of the skier and skis, the muscle activation dynamics and the muscle contraction dynamics. The dynamics was formulated in implicit form (van den Bogert et al., 2011) as f x , x ̇ , u = 0 (1) where x denotes the states of the musculoskeletal skier model, x ̇ the time derivative of the states and u the controls of the musculoskeletal skier model. Specifically, the states x of the musculoskeletal skier model were represented by the degrees of freedom and their derivatives of the multibody model of the skier and skis, the projected length of the muscle fibers (=length of the contractile element in the Hill muscle model) and the muscle activations; the controls u of the musculoskeletal skier model were represented by the neural muscle excitations of the muscles (van den Bogert et al., 2011).

We modeled the ski-snow contact using three types of forces acting on each segment of both skis (Mössner et al., 2014). First, we applied a penetration force F _p acting normal to the snow surface. The penetration force was based on an elastic force penetration relation and depended on the penetration depth and speed of the edge of the ski segment orthogonal to the snow surface as well as the edging angle of the ski. The edging angle was defined as the angle between the base surface of the ski segment and the snow surface. Second, we applied a shear force F _s acting orthogonal to the ski edge and parallel to the snow surface. The shear force provided resistance against lateral shearing and depended on the penetration depth of the ski edge. Finally, we applied a frictional force F _f acting antiparallel to the segment’s velocity with a constant friction coefficient μ = 0.1 (Heinrich et al., 2014; Mössner et al., 2014).

To analyze a turning maneuver with the present musculoskeletal model we took measurement data collected by our working group in a previous study (Filippi Oberegger, 2011) where a professional skiing instructor performed a turn to the right with the same giant slalom skies implemented in the musculoskeletal skier model. The movement of the skier was captured by a multi-camera system consisting of three cameras and a frequency of 50 Hz. In a post-processing step 23 landmarks of the skier were manually digitized and the three dimensional coordinates were reconstructed using the direct linear transformation (DLT) algorithm. Given the 23 landmarks we used OpenSim to scale the skier model and used the inverse kinematics tool to compute the kinematics of the skier (i.e., joint angles of skier at the inside and outside leg, joint angles at the lumbar joint and the translation and orientation of the pelvis segment of the skier).

Given the musculoskeletal skier model, we used a forward dynamics optimization framework to simulate the movement of the skier, track the experimental data of the skier during the turning maneuver and to compute the joint moments of the skier. Specifically, we formulated a corresponding optimal control problem (i.e., tracking problem). The task of the optimal control problem was to find the states x and controls u of the musculoskeletal skier model such that a given objective function J is minimized (van den Bogert et al., 2011). Specifically, we used the following objective function J = 1 T ∫ 0 T w 1 ‖ e r r q ‖ 2 2 ︸ tracking error + w 2 ∑ i = 1 nmus a i 2 ︸ muscle effort + w 3 ‖ x ̇ ‖ 2 2 + ‖ u ̇ ‖ 2 2 ︸ regularization including a tracking error term, a muscle effort term and a regularization term. Simulation time is denoted by T, ‖ − ‖₂ denotes the Euclidean norm, nmus the number of muscles and w ₁, w ₂ and w ₃ are weighting factors.

The first term in the objective function corresponded to the tracking error where err _q denotes the deviation of the degrees of freedom of the skier model (i.e., pelvic translation and rotation, joint angles at the lumbar, hip, knee and ankle joints) and the corresponding measurement data. The second term in the objective function corresponded to muscle effort and was used to resolve muscle redundancy (having more muscles than degrees of freedom). In the literature several criteria have been suggested (Erdemir et al., 2007). One common criterion is to use muscle activation a squared as a surrogate for muscle effort. This criterion has been used in a number of studies involving dynamic movement tasks such as jump landing (Laughlin et al., 2011), squatting (Catelli et al., 2019) or cutting (Weinhandl and O’Connor, 2017). Finally, a small regularization term was added with a small weight factor w ₃ to enhance convergence by minimizing the derivatives of the states x and controls u (Nitschke et al., 2020).

The optimal control problem was subjected to constraints due to the dynamics of the musculoskeletal skier model (i.e., muscle activation and contraction dynamics and multibody dynamics of the skier model) as well as lower and upper bounds on the states x and controls u (van den Bogert et al., 2011; Nitschke et al., 2020).

Solving an optimal control problem is computationally challenging. Recently, however, several efficient computational frameworks have been developed for solving dynamic optimization problems (e.g., Falisse et al., 2019; Nitschke et al., 2020). Similar to the approach of Nitschke et al. (2020) we implemented and solved the optimal control problem in an efficient way. We used MotionGenesis (Motion Genesis LLC, Menlo Park, CA, United States) to generate the equations of motion for the multibody dynamics of the skier model. The equations of motion were exported as C-code and imported in MATLAB via the MEX interface. In MATLAB we fused the equations for the multibody dynamics and the muscle contraction and activation dynamics and formulated the optimal control problem (i.e., tracking problem). To solve the optimal control problem, we transformed it into a constrained nonlinear programming problem (NLP) using direct collocation and the implicit Euler formula (van den Bogert et al., 2011). To increase computational speed, we provided analytical derivatives of the objective function and constraints to the NLP solver IPOPT (Nitschke et al., 2020).

To evaluate the simulation of the turning maneuvers and the associated tracking error, we first calculated the root mean squared difference (RMSD) between the joint angles of the skier derived from the measurement data and the corresponding joint angles of the skier in the simulation. Second, we compared the track of the skier where we used the ankle joint centers of the outside and inside leg as reference points (Supej et al., 2020). Finally, we computed the RMSD between the measured and simulated speed of the skier during the turning maneuver. In the data analysis, we focused on the steering phase of the turning maneuvers, where the skier is subjected to the highest loads (Klous et al., 2012, 2014). Similar to the recent study of Supej et al. (2020), we defined the beginning and end of steering phase when the turn radius of the skier was below the side cut radius of the ski (R = 32 m). After the evaluation we computed the joint moments at the hip, knee and ankle joint of the outside and inside leg as well as the lumbar joint of the skier. Joint moments were represented as internal joint moments and hip flexion, adduction and internal rotation, knee extension and ankle dorsiflexion moments were denoted as positive. To analyze the loading of the knee joints in more detail, we further computed the full 6-DOF intersegmental joint moments and forces at the knee of the outside and inside leg (van den Bogert et al., 2013), respectively, solving the Newton-Euler equations consecutively starting at the ski segments to the shank segment. Since intersegmental joint forces are not the total forces at joint and have limited utility on their own (Derrick et al., 2020), we focused on the analysis of the intersegmental knee joint moments. The intersegmental knee joint moments and forces were represented in the local coordinate system of the shank; the x-axis, y-axis and z-axis referred to the anterior-posterior, superior-inferior and medial-lateral direction.

In the forward dynamics optimization framework, the turning maneuver could be successfully simulated (an animation of the simulated turning maneuver is provided as Supplementary Material) in about 35 min of computational time on a single core of a workstation (Thinkstation 330, 3.5 GHz E-2146 CPU). In the simulation, the musculoskeletal skier model was able to track the measured kinematic data closely (Figure 2). Specifically, the RMSD between the measured joint angles and the joint angles obtained by the musculoskeletal skier model were in the range from 0.50 to 2.72° (Table 1). The lowest differences were observed for lumbar bending and knee flexion at the outside left leg; the highest differences were observed for pelvis list and hip adduction at the inside right leg.

In addition, the simulated track of the skier and the speed of the skier were in good agreement with the measurement data At the inside and outside leg, the mean deviation between the measured and simulated track was 0.025 and 0.018 m, respectively (Figure 3A). The speed of the skier increased from about 13.5 m/s to 14.5 m/s during the simulated turning maneuver and matched the measured speed with a root mean squared difference of 0.12 m/s (Figure 3B). For the speed comparison, the midpoint between the right and left hip joint center was chosen as the reference point. The turn radius of the center of mass of the skier dropped at the beginning of the steering phase to a minimum of about 16 m and remained almost constant afterwards in the range from 18 to 19 m (Figure 4A).

In the simulation of the turning maneuver, the ground reaction forces were higher on the outside leg compared to the inside leg (Figure 4B). Computing the force distribution between the inside and outside leg, about 60% of the total ground reaction force was acting on average on the outside leg. Consequently, the load on the outside leg was on average 50% higher. Peak forces reached 1.00 BW and 0.94 BW on the outside and inside leg, respectively. The local maximum at the beginning of the steering phase was induced by the skier performing an unloading-loading motion after the phase of edge change and the beginning of the steering phase (see animation provided online as Supplementary Material).

The highest internal joint moments were observed at the lumbar joint with a maximum value of 1.88 Nm/kg for lumbar extension. This was about 2.5 times larger compared to the maximum lumbar bending moment rising to 0.75 Nm/kg and about 12 times larger compared to maximum lumbar rotation moment (Table 2; Figure 5). At the outside leg, the highest internal joint moments corresponded to the hip extension moment with 1.27 Nm/kg, the knee extension moment with 1.02 Nm/kg and the ankle plantarflexion moment with 0.85 Nm/kg. At the inside leg, peak knee and hip extension moments were of similar order compared to the inside leg. The ankle plantarflexion moment and the passive boot moment, however, were 45 and 60% lower, respectively (Table 2; Figure 5).

The intersegmental knee joint moments in the frontal plane showed that primarily an internal adduction moment was acting on the knee joint of the outside leg during the turning maneuver (Figure 5), which was mainly induced by the ground reaction force passing laterally to the knee. Correspondingly, the intersegmental mediolateral force at the knee joint pointed medially to counteract the lateral component of the ground reaction force. Contrary, primarily an internal abduction moment acted on the knee joint of the inside leg during the turning maneuver (Figure 6) caused by the ground reaction force passing medially to the knee. Correspondingly, the intersegmental mediolateral force at the knee joint pointed laterally to counteract the medial component of the ground reaction force. Peak internal knee adduction and abduction moments reached 0.23 and 0.31 Nm/kg at the outside and inside leg, respectively (Table 3). In the transverse plane, the rotation moments at the knee ranged between -0.08 and 0.22 Nm/kg at the outside leg with alternating phases of internal and external rotation (Figure 6). At the inside leg, mainly an internal rotation moment was present throughout the turning maneuver with values ranging between 0.05 and 0.28 Nm/kg (Table 3).

We developed a novel three dimensional musculoskeletal model of an alpine skier with two skis and 53 kinematic degrees of freedom and applied it successfully to simulate and analyze a turning maneuver. To the authors’ knowledge, the present study incorporates the first three dimensional musculoskeletal model for analyzing turning maneuvers in alpine skiing. In previous musculoskeletal simulation studies in alpine skiing only two dimensional models were developed (Gerritsen et al., 1996; Heinrich et al., 2014, 2018). These models were applied to analyze jump landing maneuvers in downhill skiing in the sagittal plane and possible risk factor for an injury of the anterior cruciate ligament. Developing a three dimensional musculoskeletal simulation model is computational challenging due to the increased model complexity and consequently the high computational cost involved for the optimization. To reduce the computational cost, we employed recent advances in musculoskeletal modeling and simulation. In particular, we formulated the system dynamics in implicit form (van den Bogert et al., 2011) and used direct collocation with analytical derivatives of the system dynamics to solve the underlying dynamic optimization problem (Nitschke et al., 2020).

The simulation of the turning maneuver was based on a forward dynamics optimization framework (i.e., forward assisted data tracking), where measured kinematic data of an alpine skier performing a turning maneuver were tracked by the musculoskeletal skier model. In the simulation of the turning maneuver, the experimental data of the skier could be tracked closely with a RMSD below 3° at all joints. This RMSD is considered to be low, because it is well within the precision of current mobile measurement devices such as inertial measurement units (IMU) based systems (Fasel et al., 2018b) or machine learning techniques (Ostrek et al., 2019). Additionally, the track and the speed of the skier could be tracked well. Comparing the inside and outside leg, higher errors were detected at the inside leg. The higher errors could have been caused by the higher hip and knee flexion on the inside leg during the turning. This might have complicated the manual digitization process due to occlusions and consequently reduced the accuracy of the corresponding measurement data.

In the forward dynamics optimization framework, kinematic data obtained by video-based stereophotogrammetry and a multi-camera setup were used as input data. IMU based systems in combination with a Global Navigation Satellite System (GNSS) or computer vision and human pose estimation have been shown to be promising approaches for capturing the kinematics of a skier during turning maneuvers on the ski slope (e.g., Fasel et al., 2018b; Ostrek et al., 2019). Since only measured kinematic data are mandatory in the forward dynamics optimization framework (see also next section), the present musculoskeletal simulation model might also be combined with either of these approaches in future research without the need to provide ground reaction force data captured by mobile force platforms or pressure insoles. This might open up new opportunities regarding the analysis of the loading of the skier during turning maneuvers in the natural environment if only kinematic data are available. If measured ground reaction force, however, are additionally available as input data, these data might also be tracked in the optimization framework. Specifically, in the objective function of the optimization framework the deviation of the simulated and measured ground reaction forces might be taken into account to track the measured ground reaction forces as good as possible (van den Bogert et al., 2011; Nitschke et al., 2020).

In the present study we used a forward dynamics optimization framework to compute the joint moments during a tuning maneuver. In the literature there are only a few studies analyzing the loading of the skier during turning maneuvers dynamically (Klous et al., 2012, 2014; Hirose et al., 2013; Lee et al., 2017). However, in all of these studies the authors used an inverse dynamics approach to compute the joint moments at the lower extremities. Compared to the classical inverse dynamics analysis, the present framework has four major advantages.

First, using a forward dynamics optimization framework the obtained kinematics of the skier as well as the corresponding ground reaction forces are dynamically consistent. Consequently, no residual forces and torque have to be added such that the equations of motion of the skier are satisfied (Bailly et al., 2021). This is important since residual forces and moments do not exist in reality and affect the computation of kinetic variables such as joint moments, power and work (Hatze, 2002; Faber et al., 2018). Contrary, using an inverse dynamics approach the validity of the computed joint moments is impaired depending on the size of the introduced residual forces and moments (Faber et al., 2018).

Second, optimization of a forward dynamic model with data tracking and effort minimization also has the advantage that the data being tracked can be any number of variables. The number of measurements can be overdetermined (more measurements than kinematic degrees of freedom and external loads), or underdetermined (fewer measurements than kinematic degrees of freedom and external loads). This makes it possible to perform a full dynamic analysis without external loads (i.e., ground reaction forces) provided for example by instrumented skis, as demonstrated previously in a planar analysis of jump landing in skiing (van den Bogert et al., 2011) or in the present three dimensional analysis of a turning maneuver in skiing. In contrast, classical inverse dynamic analysis requires the measurement of all kinematic degrees of freedom and external loads, not more and not less.

Third, the computed joint moments are less sensitive to errors in the measurement data, since only the measured kinematic data (i.e., joint angles, pelvic translation and orientation) are tracked in the simulation. First or second derivatives of the measured kinematic data are not required. Contrary, inverse dynamics requires the second derivative of the measured kinematic data, which amplifies measurement errors (Cahouët et al., 2002; Pàmies-Vilà et al., 2012).

Fourth, the computed joint moments are constrained to stay within physiological limits. The physiological limits were induced by the musculoskeletal model, which included a three element Hill-type muscle model with activation and contraction dynamics. Contraction dynamics incorporates the force-length-velocity characteristics of the muscle, the active state as well as the maximum isometric force (van den Bogert et al., 2011) and limits the acting muscle forces and the corresponding joint moments. In addition, the activation dynamics limits the change of the corresponding joint moments, incorporating time constants for muscle activation and deactivation in the dynamics equation (Bailly et al., 2021).

Based on the musculoskeletal model and the forward dynamics optimization framework we computed the joint moments at the lumbar joint as well as the hip, knee and ankle joints of the outside and inside leg of the skier during the turning maneuver. At the outside leg, highest lower-limb joint moments were identified at the hip joint (1.27 Nm/kg, hip extension), followed by the knee joint (1.02 Nm/kg, knee extension) and ankle joint (0.85 Nm/kg, ankle plantarflexion). Knee and hip extension moments were similar at the outside and inside leg, although the ground reaction forces were on average about 50% higher on the outside leg. This can be explained by the increased knee and hip flexion on the inside leg, which required higher activation of the knee and hip extensor muscles. The ankle plantarflexion moment and also the passive boot moment were lower on the inside leg, which indicated that the skier was pushing more against the shaft of the ski boot at the outside leg.

In accordance with the present study, the hip extension moment was reported as the highest joint moment in the kinetic studies of Hirose et al. (2013) and Lee et al. (2017) analyzing turning maneuvers in skiing based. In both studies, the authors used on an inverse dynamics approach to compute the joint moments given kinematic data from an IMU based system and measured ground reaction forces. Specifically, in the study of Lee et al. (2017), the peak knee extension moment (0.5 Nm/kg) and ankle plantarflexion moment (1.1 Nm/kg) were roughly of the same order of magnitude as in the present study. Strikingly, however, the peak hip extension moment reached 6 Nm/kg. In the study of Hirose et al. (2013), the hip extension moment reached 12 Nm/kg assuming a mass of the skier of 75 kg. Compared to the present study, these reported peak hip extension moments are a factor of 5 and 10, respectively, higher and are likely to be unrealistic high. In particular, these values exceed reported maximum voluntary joint moments by about 70 and 240%, respectively, if we take the data of the age group 19–25 in the study of Anderson et al. (2007) as reference data and add one standard deviation to the mean value. Unfortunately, Hirose et al. (2013) and Lee et al. (2017) did not present any information about the external forces acting on the skier (i.e., ground reaction forces or speed and turn radius of the skier) or the kinematics of the skier, which makes a more detailed comparison impossible. On the other hand, similar hip extension moments as in the present study were reported in the study of van den Bogert et al. (1999). Nine male subjects were instrumented with a 12-channel accelerometer system mounted on the upper body and the hip extension moment was computed through an inverse dynamics analysis. The authors reported an average peak hip extension moment of 1.75 Nm/kg during long turns on a flat slope assuming a mean mass of the skiers of 75 kg. The skiers were asked to load the outer ski only and lift the inner ski, which shifted the entire load to the outer ski. In contrast, only 60% of the external load was acting on the outside leg of the skier during the turning maneuver in the present study.

The internal knee joint moments showed that in the frontal plane primarily an adduction moment was acting on the knee joint of the outside leg in combination with alternating phases of an internal and external rotation moment in the transverse plane as well as a knee extension moment. Klous et al. (2012) studied the loading of the knee joint during a carved and skidded turn. The turn radius of the skier was about 10 m in both turns, while the mean speed during the carving turn was higher than during the skidded turn (13.9 m/s vs 10.4 m/s). The authors determined peak knee extension moments up to 8.35 and 4.07 Nm/kg for the skidded and carved turn, respectively, based on an inverse dynamics approach. In addition, they computed peak adduction moments up to 5.70 and 5.75 Nm/kg in the frontal plane and external rotation moments up to 6.85 and 2.75 Nm/kg in the transverse plane for the skidded and carved turn. All of these values are significantly higher compared to the values of the present study. The differences might be explained in part by the decreased turn radius of the skier in the study of Klous et al. (2012), which induced a higher external loading onto the skier. However, due to measurement noise and the limitations of the inverse dynamics approach, the computed joint moments might have been overestimated exceeding physiologically plausible values. For example, if we compare the peak knee extension moments with literature data such as maximum voluntary knee extension moments (Anderson et al., 2007) or maximum isometric knee extension moments (Pincivero et al., 2004; Domire et al., 2011), these values are unrealistic high. In addition, the reported knee abduction moments are about three times higher than the assumed injury threshold of 125 Nm valgus moment (= 1.67 Nm/kg assuming a mass of 75 kg) in the study of McLean et al. (2008), although they did not investigate an injury prone situation.

Interestingly, in the present study the peak joint moment at the lumbar joint exceeded the peak values at the lower limbs during the turning maneuver. In particular, the lumbar extension moment rose up to 1.88 Nm/kg, while the lateral bending and rotation moments were remarkably lower. Consistent with the results of the present study, the highest joint moment was observed at the lumbar joint for lumbar extension in the study of Hirose et al. (2013). High values at the lumbar joint imply that the lower back of the skier is subjected to high loads during turning maneuvers. These high values may be linked to lower back pain, which is a common overuse injury in alpine skiing (Spörri et al., 2018). Furthermore, it has been reported that the combination of lumbar flexion, lateral bending and axial rotation amplifies the loading at the lower back (Spörri et al., 2018). Further studies, however, and a more sophisticated model of the lower back are necessary to quantify the internal loading at the lower back in more detail.

Some limitations of the present study have to be mentioned. First, in the simulation of the turning maneuver we did not track the movement of the arms of the skier. The reason was to reduce the complexity of the model and to decrease computational time, which is one of the big challenges in three dimensional musculoskeletal simulations. However, since we included the mass and inertia properties of the arm in the model and assumed a mean posture of the arms in front of the skier, the impact on the computed joint moments is expected to be low.

Second, we did not implement a detailed spine model, but used a single lumbar joint at the lower back. Consequently, the present simulation model is expected to provide only basic features regarding the loading of the lower back of the skier turning maneuvers. While these basic features might contribute to the understanding of lower back pain, which is a common overuse injury in alpine skiing (Spörri et al., 2018), a more detailed spine model might provide further insight.

Third, in the present simulation study we analyzed data of a professional skiing instructor performing a giant slalom turning maneuver. Changing the characteristics of equipment, the present simulation model might also be used to analyze turning maneuvers in other disciplines such as slalom, super-G (super giant slalom), or downhill skiing. Furthermore, the present simulation model might also be used to analyze jump landing maneuvers in super-G and downhill, which have been identified as a common situation leading to injury (Gilgien et al., 2014). Regarding jump landing maneuvers the present simulation model might extend the current knowledge derived from two dimensional simulation models (Heinrich et al., 2014) regarding the loading of the skier in the frontal and transverse plane.