Effects of load mass and position on the dynamic loading of the knees, shoulders and lumbar spine during lifting a musculoskeletal modelling

Musculoskeletal models may enhance our understanding of the dynamic loading of the joints during manual material handling. This study used state-of-the-art musculoskeletal models to determine the effects of load mass, asymmetry angle, horizontal location and deposit height on the dynamic loading of the knees, shoulders and lumbar spine during lifting. Recommended weight limits and lifting indices were also calculated using the NIOSH lifting equation. Based on 1832 lifts from 22 subjects, we found that load mass had the most substantial effect on L5-S1 compression. Increments in asymmetry led to large increases in mediolateral shear, while load mass and asymmetry had significant effects on anteroposterior shear. Increased deposit height led to higher shoulder forces, while the horizontal location mostly affected the forces in the knees and shoulders. These results generally support the findings of previous research, but notable differences in the trends and magnitudes of the estimated forces were observed.


Introduction
Manual material handling (MMH) constitute a substantial proportion of the work performed in many industries, such as manufacturing, construction and retail (Heran-Le Roy et al., 1999). The term typically encompass the acts of lifting, lowering, carrying, pushing and pulling various materials (Dempsey, 1998), but can also be an underlying aspect of other commonly used terms to describe physical workplace factors, such as heavy physical work, forceful exertions and awkward postures. Numerous epidemiological studies have identified MMH as a risk factor for developing work-related musculoskeletal disorders: for example, MMH has been associated with musculoskeletal disorders in the lower back (da Costa and Vieira, 2010;Coenen et al., 2014;Lötters et al., 2003;Punnett and Wegman, 2004;NIOSH, 1998;Burdorf and Sorock, 1997;Kuiper et al., 1999;Hoogendoorn et al., 1999), knees (da Costa and Vieira, 2010; Reid et al., 2010;Fransen et al., 2011) and shoulders (Van Rijn et al., 2010;Mayer et al., 2012). Hence, the disorders resulting from MMH pose a substantial burden on society, both in terms of the human and economic costs (Dempsey, 1998;Dempsey and Hashemi, 1999).
For these reasons, substantial efforts have been made to understand how the body is loaded during MMH and how some loads may increase the risk of injury to the involved joints. Biomechanical models have played a major role in this regard. In the 1960s, Morris et al. (1961) made one of the earlier attempts to calculate the forces on the spine using a simple static, sagittal-plane biomechanical model (Potvin, 2008). In the following decade, Chaffin et al. (1975) further developed this model, which would later be used in the NIOSH Work Practices Guidelines (Nelson et al., 1981) -one of the most cited and important works in the field of occupational biomechanics (Chaffin, 2009). In the NIOSH guidelines, the biomechanical model was used to estimate the compressive force on the L5-S1 joint, which was evaluated against failure tolerance data of cadaver discs (Nelson et al., 1981). The NIOSH guidelines were later revised to account for hand coupling and asymmetrical lifting (Waters et al., 1993(Waters et al., , 1994. This version has since been shown to be a reasonable predictor of injury risk to the lower back during lifting (Marras et al., 1999a;Waters et al., 2011).
In the years following the publication and implementation of static, sagittal-plane biomechanical models, it became clear that these models generally tended to underestimate the compressive forces in the lumbar spine as they ignore the acceleration components and inertial loads of the body segments as well as the co-contraction of the trunk muscles (Garg and Kapellusch, 2009). For example, studies by Garg et al. (1982), McGill and Norman (1985) and de Looze et al. (1994) all showed that either the compressive forces or joint moments in the lumbar spine were significantly higher when employing 2-D dynamic models, even without the representation of trunk muscles. In 1986, McGill andNorman (1986) made a major advancement with the publication of a 3-D EMG-driven dynamic model with representations of 7 ligaments and 48 muscles (Potvin, 2008). Since then, several other detailed 3-D dynamic models have been published, e.g. by Granata and Marras (1995) and Kingma et al. (1996).
Today, scientific and technological developments, particularly with respect to computer processing power, has enabled the application of even more detailed and versatile biomechanical models. In particular, modelling and simulation software such as OpenSim (Delp et al., 2007) and The AnyBody Modelling System (AMS) (Damsgaard et al., 2006) provides a continuously growing library of complex musculoskeletal models. In AMS, musculoskeletal models of the various body parts, developed by many different research groups independently over the last three decades, are integrated in a full-body model with upwards of a 1000 muscle elements to provide a comprehensive and extensively validated tool to estimate dynamic muscle and joint reaction forces (Lund et al., 2020). For example, the lumbar spine model (de Zee et al., 2007;Han et al., 2012) has been shown to produce spinal forces in close agreement with in vivo data (Bassani et al., 2017;Rajaee et al., 2015). Furthermore, as the full-body model also includes detailed cadaver-based models of the shoulders and lower extremities, it enables an accurate estimation of the forces in all the major joints involved in MMH. However, few studies to date have used the AMS for the analysis of MMH (Larsen et al., 2020;Behjati and Arjmand, 2019;Stambolian et al., 2016;Koblauch, 2016;Wagner et al., 2007;Skals et al., 2021); two of these examples were validation studies (Larsen et al., 2020;Stambolian et al., 2016), while three others were based on only one subject (Behjati and Arjmand, 2019;Koblauch, 2016;Wagner et al., 2007). A recent field study by the authors (Skals et al., 2021) is the only example in the literature, where musculoskeletal models were applied for a comprehensive analysis of MMH. However, these models have yet to be implemented for a systematic analysis of MMH in a laboratory setting, which could provide valuable and accurate reference data of the dynamic forces in multiple joints simultaneously during different lifting conditions. Therefore, the aim of the present study was to determine the effects of load mass (LM), asymmetry angle (AA), horizontal location (HL) and deposit height (DH) on the dynamic loading of the knees, shoulders and lumbar spine during lifting using state-of-the-art, full-body musculoskeletal models. In addition, The Revised NIOSH Lifting Equation (Waters et al., 1993(Waters et al., , 1994 was used to calculate recommended weight limits (RWLs) and lifting indices for all tested conditions, hereby providing an opportunity to evaluate the appropriateness of the equation to assess the risk of low back injury during lifting. The main objectives of this approach were to provide a comprehensive analysis of the loads in the involved joints during different lifting conditions using one of the most advanced musculoskeletal models available and facilitate the evaluation of injury risk during occupational lifting based on a collective assessment of multiple joint loads.

Subjects
A convenience sample of 22 healthy, voluntary subjects (16 male and 6 female, age: 30 ± 10, height: 178.1 ± 11.3 cm, mass: 80.6 ± 12.1 kg) participated in the study. The study followed the guidelines of The North Denmark Region Committee on Health Research Ethics and the subjects provided written informed consent. Data were collected in January and May 2019.

Instrumentation
Full-body kinematics and ground reaction forces and moments were obtained by measuring the trajectories of 43 passive reflective markers using eight infrared Oqus cameras (Qualisys AB, Göteborg, Sweden) and from two forces plates (one under each foot) embedded in the laboratory floor (AMTI, Watertown, MA, USA), respectively. The marker and force plate data were sampled at 120 Hz and digitally low-pass filtered with cut-off frequencies of 6 and 15 Hz, respectively, using a second-order, zero-phase Butterworth filter. Thirty-nine markers were placed on the subjects (see supplementary material) and one marker placed on each corner of a custom-built box used during the experiment (see Fig. 1). The box was made of a combination of steel and aluminum with a steel pole screwed into the base plate at the center of the box, where iron weight plates were placed. As can be seen in Fig. 1, force transducers were placed under the handles with the cables connecting the transducers to the amplifier held by a string in the ceiling to avoid the cables affecting the execution of the lifts. However, due to a systematic error in the measurements, the resulting force data could not be used. Instead, the external forces and moments of the box were estimated using contact elements (see section 2.4).  -55) and 60 cm (HL-60), and five lifts with DHs of 30 (DH-30), 60 (DH-60), 90 (DH-90), 120 (DH-120) and 150 cm (DH-150) above floor level. Note that the LM-10 and HL-45 conditions were identical and therefore, based on the same measurements. Except for the DH condition, the chosen lifting factors were inspired by the multipliers used in The Revised NIOSH Lifting Equation to determine the risk of developing lifting-related low back pain (Waters et al., 1993(Waters et al., , 1994. Hence, the definition of the distances and angles between the subjects and the starting locations of the lifted load resembled the specifications outlined in the NIOSH lifting equation. The variations in LM are self-explanatory, while the definitions of the load positions relative to the subjects are described in the following.

Data collection
The AAs were defined as the angle between the sagittal line of the Percentage of bodyweight per second RWL Recommended weight limit body at the midpoint between the malleoli and the sagittal line of the lifted load at the midpoint between the hands. The HLs were defined as the distance between the midpoint of the line between the malleoli and the midpoint between the hands. Finally, the DHs were defined as the vertical height of the shelves above the floor. The distance from the floor to the midpoint between the hands were approximately 10 cm higher (distance from the center of the hands to the bottom edge of the load). For the LM and HL conditions, the box was placed on a EUR-pallet with a height of 15 cm, meaning that the distance between the hands and floor were approximately 25 cm at the initiation of the lifts. For all lifts related to the LM condition, the horizontal location was 45 cm. For the AA and DH conditions, the load was placed on a custom-made wooden bench, also with a height of 15 cm. Both series of lifts were initiated with a HL of approximately 35 cm, while the HL was approximately 68 cm when the box was placed on the shelves during the DH condition. For the AA, HL and DH conditions, the LM was 10 kg. The lifting procedures for AA, HL and DH are illustrated in Fig. 1. For the AA, HL and LM conditions, the subjects were instructed to lift the box to an upright position with their hands slightly above waist height and subsequently lower the box down to the starting position. For the DH condition, the subjects were asked to lift the box from the starting position and place it on the appropriate shelf. No instructions were given regarding lifting technique or pace, but the subjects were encouraged to lift in a controlled fashion. Four repetitions were performed for each lifting condition, totaling 84 lifting trials for each subject. The order of the lifting conditions were counter-balanced to avoid any order effects, e.g. due to physical or mental fatigue.

Musculoskeletal modelling
Musculoskeletal models were developed in The AMS v. 7.3 (AnyBody Technology A/S, Aalborg, Denmark) using the Plug-in-gait-MultiTrial-StandingRef template from the AnyBody Managed Model Repository v. 2.3 (Lund et al., 2020). The model included 947 muscle elements and 86 degrees-of-freedom (DOF): 2 × 3 DOF at the ankle joints, 2 × 1 DOF at the knee joints, 2 × 3 DOF at the hip joints, 6 DOF at the pelvis, 3 DOF between thorax and pelvis, 12 × 3 DOF for the cervical, thoracic and lumbar spine, 2 × 2 DOF at the elbow joints, 2 × 3 DOF at the glenohumeral joints, 2 × 3 DOF at the sternoclavicular joints, 2 × 2 DOF at the wrist joints, 1 DOF at the neck joint and 6 DOF between the hands and box. An example of the modelling procedure is illustrated in Fig. 2. The lumbar spine model was based on the work of Hansen et al. (2006), de Zee et al. (2007) and Han et al. (2012), and consisted of seven rigid segments, representing the lumbar vertebrae, thoracic spine and sacrum. It was modelled with non-linear disc stiffness in the lumbar region with representations of seven ligaments and intra-abdominal pressure (Han et al., 2012). The thoracic part of the spine was modelled as a rigid segment with spherical joints between T1 and C7 as well as between T12 and L1. The model followed a kinematic rhythm that distributes the trunk motion over the vertebral bodies through a coupled mechanism and was actuated by 188 muscle elements.
The shoulder and arm model was based on the work of van der Helm et al. (1992) and Veeger et al. (1991;. The shoulder and arm model included 146 muscle elements: the shoulder model included representations of the deltoideus, subscapularis, infraspinatus, supraspinatus, teres minor, pectoralis major among others, while the arm model further included representations of for instance, the brachialis, biceps and triceps brachii, as well as multiple wrist flexor, extensor, supinator and pronator muscles. Muscles spanning the shoulder joint with complex wrapping behaviour (e.g. deltoideus) have insertion points on the bones and wrap over analytical geometric shapes, such as ellipsoids or spheres. The lower extremity model was based on the cadaver study of Carbone et al. (2015) as well as the study of De Pieri et al. (2018). The model includes representations of 169 muscle-tendon elements in each leg with multiple muscles spanning the knee joints, e.g. the gastrocnemius and biceps femoris. In the configuration used in the present study, the knee was modelled as a hinge joint with a fixed rotation center and axis, and the patella tendon defined as a non-deformable element connecting the patella to the tibia.
A 3D computer-aided design model of the box was created in Solid-Works v. 2017.5.0 (Dassault Systems, Vélizy-Villacoublay Cedex, France). The mass of the different box configurations were based on measurements made during data collection and used along with the geometry to estimate the inertial properties in SolidWorks.
Twelve contact elements with a high strength (400,000 N and Nm) were defined between the hands and box to estimate the external forces and moments (Larsen et al., 2020;Skals et al., 2021), while small residual forces and moments with a strength of 10 N and 10 Nm, respectively, were placed on the pelvis to improve numerical stability. The actuation of the contact and residual elements were solved as part of the muscle recruitment algorithm (see section 2.6).

Model scaling and kinematics
The geometric and inertial parameters of the musculoskeletal models were scaled to the subjects by applying a length-mass-fat scaling law (Rasmussen et al., 2005), and the total body mass distributed to the individual segments using the regression equations presented in Winter et al. (2009).
For each subject, a single trial (DH-120) was initially used to determine segment lengths and marker positions using the optimization method of Andersen et al. (2009). The optimized segment dimensions and marker positions for each subject were then saved and used for the analysis of all other trials for that subject. Kinematic analysis was performed by minimizing the least-square difference between model and experimental markers over the whole trial duration (Andersen et al., 2010).

Muscle recruitment
The muscle, joint, contact and residual forces were distributed by solving a second-order quadratic optimization problem (Damsgaard et al., 2006): G is the second-order quadratic objective cost function, n ( i is the strength of the residual force. C is the coefficient matrix for the dynamic equilibrium equations, f is all the unknown muscle and joint reaction forces (JRFs) and d contains all the external and inertial forces. The muscle strengths were determined from the physiological cross-sectional area and the length-mass-fat scaling law (Rasmussen et al., 2005;Frankenfield et al., 2001) with the non-negativity constraints dictating that the muscles can only pull. The muscles were modelled without contraction dynamics.

Data analysis
Peak L5-S1 axial compression (A-C) force and impulse, anteroposterior (A-P) and mediolateral shear (M-L) forces as well as the peak resultant JRFs for the left and right knee and shoulder (glenohumeral) joints were selected for further analysis. The forces were normalized to percentage of bodyweight (%BW) and impulse to percentage of bodyweight per second (%BWs). Trunk flexion, lateral bending and rotation as well as knee and shoulder flexion joint angles were also extracted from the musculoskeletal models and are presented in the supplementary material (Sup. Figs. 1 and 2). In addition, The Revised NIOSH Lifting Equation was used to calculate recommended weight limits (RWLs) and lifting indices for all conditions following the procedures outlined in Waters et al. (1993;. For LM, HL and AA, the lifting cycles were defined as the instant when the subjects lifted the box from its base to the instant it was returned to the starting position. For DH, the lifting cycle ended at the instant the box was placed on the appropriate shelf. The peak JRFs were analyzed for the second half of the lifting cycle during the DH conditions to determine the effect of the end location rather than the starting location, except for calculating the lifting indices. The selected variables were resampled to 101 data points (one lifting cycle).

Statistical analysis
The statistical analyses were performed in SAS v. 9.4 (SAS Institute Inc., Cary, NC, USA): repeated measures linear mixed models (Proc Mixed, SAS) were used to test if any significant differences existed between the different levels for each condition separately with the peak JRFs and L5-S1 A-C impulse as the dependent variables and the condition levels included as fixed effects. The subjects were included as random effects to consider repeated measures. The covariance structure was set to "Variance Components" and the model fit using restricted maximum likelihood estimation. Normal distribution of the residuals and homogeneity of variance were ensured by visually inspecting residual diagnostic plots. The data are presented as least square means with 95% confidence intervals based on a Satterthwaite approximation, and differences of least square means reported for p < 0.05.

Results
From the 22 subjects, a total of 1832 lifting trials were included in the analysis. Sixteen trials were either missing or excluded: four trials of HL-35 and one trial of DH-30 were missing for subject 22 and 7, respectively, while four trials of DH-150 were missing for subject 12, as the subject was not able to reach this shelf height. All four trials of HL-35 for subject 16, two trials of AA-30 for subject 21 and one trial of LM-10 for subject 1 were excluded due to random errors in the force plate data.
Time-series curves of the mean JRFs for all lifting conditions are illustrated in Fig. 3a,b, 4a and b, while least square means with confidence intervals for each condition and indications of significant differences between levels listed in Tables 1 and 2. The NIOSH lifting indices and recommended weight limits for all tested conditions are presented in Table 3. The linear mixed model analysis for each lifting condition showed significant differences for all outcome variables overall (p < 0.0001). The differences of least square means analyses showed further significant differences between individual levels for all conditions (see Tables 1 and 2), which are specified in the following.

Horizontal location
The L5-S1 A-C force ranged from 502 to 549 %BW and significant differences were found between all levels with the exception of HL-40 vs. HL-45. The L5-S1 A-P forces ranged from 49.6 to 54.5 %BW with no significant differences between the three highest HLs (i.e. HL-50, HL-55 and HL-60) as well as between the two lowest HLs (i.e. HL-30 and HL- 35), but more notable differences between the lowest and highest levels. The L5-S1 M-L forces were low in general (from 6.5 to 7.5 %BW), but significant differences were found between the four lowest levels (i.e. HL-30, HL-35, HL-40 and HL-45) and HL-50 and HL-60 (p < 0.05). The L5-S1 A-C impulse was not substantially affected by the increments in HL, ranging from 1376 (HL-45) to 1490 %BWs (HL-55). For the resultant JRFs in the knee (left: 362-466 %BW, right: 388-503 %BW) and shoulder joints (left: 52.6-103.6 %BW, right: 53.4-102.2 %BW), the data generally showed higher forces with increments in the HL with a notable difference between the magnitude of forces in the left and right knee. Furthermore, a substantial increase in the differences between levels was found for the shoulder forces when the HL increased to 45 cm and above (see Table 2).

Deposit height
For the L5-S1 JRFs, the most notable trends in the data were that elevations in the DHs resulted in significant reductions in the A-C (from 515 to 387 %BW) and A-P shear forces (from 42.2 to 34.0 %BW), particularly when the DH was increased to 90 cm and above, while the changes in M-L shear were negligible (from 3.3 to 6.2 %BW). The L5-S1 A-C impulse significantly increased with all increments in DH (from 782 to 1123 %BWs), particularly when the DH was increased to 120 and 150 cm, with the exception of DH-60 vs. DH-90. The largest differences were found for the shoulder forces, which more than doubled from the lowest (left: 115 %BW, right: 111 %BW) to the highest shelf level (left: 240 % BW, right: 245 %BW). The knee JRFs (left: 162-243 %BW, right: 181-263 %BW) were significantly reduced from DH-30 (253 %BW on average) to DH-90 (172 %BW on average), but increased slightly from DH-90 to DH-120 (193 %BW on average) and DH-150 (202 %BW on average).

Load mass
For all outcome variables, increments in LM resulted in significantly higher JRFs for nearly all levels. Specifically, all increments in LM resulted in significantly higher L5-S1 A-C (from 472 to 672 %BW) and A-P shear forces (from 47.5 to 67.9 %BW). Similar tendencies were found for the M-L shear force (from 6.2 to 9.5 %BW), although no significant differences were established for LM-5 vs. LM-10, LM-10 vs. LM-15 and LM-10 vs. LM-20. The L5-S1 A-C impulse also significantly increased for each increment in LM, ranging from 1216 (LM-5) to 2148 %BWs

Asymmetry angle
For the L5-S1 A-C force (from 513 to 543 %BW), the lowest (AA-15) and highest level (AA-75) were significantly different from all other levels, while no significant differences were established between AA-30 (530 %BW), AA-45 (531 %BW) and AA-60 (531 %BW). For the A-P (from 51.9 to 67.8 %BW) and M-L shear forces (from 14.1 to 40.0 %BW), every increment in AA resulted in significantly higher forces (p < 0.0001). In particular, the M-L shear forces were substantially higher compared with the other lifting conditions with the force increasing approximately 6 %BW for every increment in AA. The L5-S1 A-C impulse was identical for AA-15 and AA-30 (1478 %BWs), but significantly higher when the AA was increased to 45 cm and above (from 1553 to 1747 %BWs). As the box was moved to the right of the subjects during the AA lifts, the left knee JRFs expectedly decreased with increments in the AA (from 348 to 263 %BW), while the forces in the right knee increased substantially (from 439 to 679 %BW). The forces in the right knee significantly increased for every increment in AA, reaching the highest forces in the right knee across all conditions for AA-75. Finally, every increment in AA resulted in significantly higher bilateral shoulder JRFs (left: 63.1-106.6 %BW, right: 57.1-113.2 %BW) of 12.5 %BW on average. Fig. 4b. Time-series curves of the mean knee and shoulder (glenohumeral) resultant joint reaction forces (JRF) for the left (L) and right (R) side of the 22 subjects, normalized to percentage of bodyweight (%BW), during the lifts with a load mass (left) of 5 (red), 10 (green), 15 (blue), 20 (cyan) and 25 kg (magenta) as well as asymmetry angles (right) of 15 (red), 30 (green), 45 (blue), 60 (cyan) and 75 • (magenta). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

NIOSH lifting indices
For the HL conditions, the RWLs ranged from 7.6 (HL-60) to 15.3 kg (HL-30), while the lifting indices ranged from 0.66 to 1.31. The L5-S1 A-C forces were estimated to between 3969 and 4346 N with HL-45 producing the lifting index closest to the action limit of 1.0, which corresponded to an A-C force of 4132 N. The RWLs and lifting indices calculated over the whole lifting cycle for the DH conditions ranged from 12.8 to 14.7 kg and 0.68 to 0.78 with the A-C forces ranging from 3928 (DH-120) to 4105 N (DH-30). For the variations in AA, the RWLs ranged from 9.9 (AA-75) to 12.5 kg (AA-15) with lifting indices and A-C forces of 0.80-1.00 and 4060-4295 N. AA-75 was the only condition for which the lifting index (1.00) equaled the action limit, which resulted in an A-C force of 4295 N. Finally, the LM conditions produced lifting indices from 0.49 to 2.46 and a RWL of 10.2 kg for all load instances, which resulted in A-C forces from 3734 (LM-5) to 5316 N (LM-25).

Table 1
Least square means of the peak L5-S1 axial compression force (A-C) and impulse as well as anteroposterior (A-P) and mediolateral (M-L) shear forces*. The forces are normalized to percentage of bodyweight with 95% confidence intervals in parentheses, while impulse is presented as percentage of bodyweight per second. Significant differences (p < 0.05) between condition levels are indicated with numbers in superscript.

Table 2
Least square means of the peak knee and shoulder (glenohumeral) resultant joint reaction forces (JRFs) for the left (L) and right side (R)* normalized to percentage of bodyweight with 95% confidence intervals in parentheses. Significant differences (p < 0.05) between condition levels are indicated with numbers in superscript.

Discussion
In this study, we used state-of-the-art musculoskeletal models to determine the effects of well-known lifting factors on the dynamic loading of the knees, shoulders and lumbar spine. The main findings were that all the tested conditions had a significant overall effect on the peak loading of the involved joints with particular lifting factors exhibiting a substantial influence on one or more joint forces. Increments in the HL had minor effect on the L5-S1 forces, but a substantial effect on the peak forces in the knee and shoulder joints. The DH had a significant effect on the peak forces in the knees, shoulders and lumbar spine. The increments in LM had the most clear positive linear relationship across all outcome variables, most notably with respect to the L5-S1 A-C and A-P shear forces. Finally, each increment in the AA led to significantly higher peak L5-S1 A-P and M-L shear forces, while also showing significant increases of the peak forces in the right knee as well as both shoulders.
Variations in LM had the most substantial influence on the peak L5-S1 A-C and A-P shear forces. Each 5 kg increment in LM resulted in an average increase in A-C and A-P shear force of 50.0 %BW (396 N) and 4.9 %BW (38.8 N), respectively, with significant differences between all condition levels. Furthermore, the L5-S1 A-C impulse also increased substantially with each increment in LM. This is in line with the results of previous research showing increased loads in the lower back with increments in LM (Granata et al., 1999;Davis et al., 2010;Lavender et al., 2003;Plamondon et al., 2012). Similarly, the forces in the knees and shoulders increased significantly with every increment in LM reaching average peak forces of 547 %BW (4329 N) and 127 %BW (1005 N) for LM-25, respectively. Interestingly, Schipplein et al. (1990) found an inverse relationship between the dynamic peak flexion-extension moment in the right knee (from 53 to 13 Nm) and L5-S1 joint (from 225 to 344 Nm) with similar increments in LM (approximately 5-25 kg). In the present study, the peak resultant JRFs in the right knee (from 388 to 503 %BW), L5-S1 A-C and A-P shear forces all increased with increments in LM. Schipplein et al. (1990) attributed the decreased moment in the knee to the subjects changing to a so-called stoop-lifting technique with increments in load. This change in lifting technique was not observed in the present study, as the subjects maintained a squat-lifting technique with a similar degree of knee and trunk flexion for all load instances (see Sup. Fig. 2a and b). This surely contributed to the high degree of linearity between the increments in LM and the JRFs, as the subjects did not attempt to alter their technique in response to the increased load.
The variations in HL had a much smaller effect on the L5-S1 forces than LM with a difference between the closest (HL-30) and farthest location (HL-60) of 47.6 %BW (377 N) and 4.9 %BW (39 N) for the A-C and A-P shear forces, respectively, corresponding to the average of the differences observed between each increment in LM. Hence, although most of the 5 cm increments in HL were statistically significant with respect to the A-C force, these differences have minor practical significance when compared with the effect of LM. However, this may only be the case at relatively light loads (e.g. 10 kg), as an increased LM would presumably amplify the effects of the HL on the L5-S1 A-C and A-P shear forces. The increments in HL had a more pronounced effect on the forces in the knees and shoulders: for the knees, the forces increased with 109 %BW (866 N) between HL-30 and HL-60, while the forces in the shoulders almost doubled from 53.0 %BW (420 N) to 103 %BW (814 N). Interestingly, instead of straightening the knees and leaning more forward in response to the increased HL, the subjects chose to further flex their knees and shoulders, hereby maintaining a similar degree of trunk flexion for all condition levels (see Sup. Fig. 1a and b). This choice of lifting technique may explain why the increments in HL had such a pronounced effect on the knee and shoulder JRFs, but only a minor effect on the L5-S1 JRFs. Furthermore, this may also have led to the subjects increasing the horizontal force posteriorly when picking up the box, as they did not move their center of gravity forward in response to the increased HLs, hereby potentially reducing the effect of the HL increments on the A-C and A-P shear forces. However, this was not possible when repositioning the load on the pallet, which led to notably larger differences in A-C force between condition levels at the deposit phase compared with the pick-up phase (see Fig. 3a).
Several previous modelling studies have found increased low back Table 3 Recommended weight limits (RWL) and lifting indices (LI) calculated using The Revised NIOSH Lifting Equation with the equation multipliers listed, specifically the initial horizontal (HL) and vertical location of the midpoint between the hands (VL), vertical travel distance (D) and asymmetry angle (AA). The lifting frequency was set to one lift per minute for up to 2 h of a workday, while the highest point of the load was set to 60% of the subjects' average body height for calculating D during the HL, LW and AA conditions. In addition, the table shows the load mass (LM) and least square means of the peak axial compression force (A-C) in newton (N)*. *Note that for the DH condition, the peak A-C forces were determined for the whole lifting cycle unlike the results in Tables 1 and 2. loading with increments in the HL (Schipplein et al., 1995;Faber et al., 2011;Lavender et al., 1999), while other studies did not show this effect (Marras et al., 1999b;Faber et al., 2007). Faber et al. (2007) suggested that the different findings may be a result of the subjects being permitted to more freely adapt their lifting technique in response to the increased HLs in the studies showing no significant change in low back loading (Marras et al., 1999b;Faber et al., 2007). One example from the literature that analyzed knee loads in addition to low back loads for different HLs (from 20 to 60 cm) is the study by Schipplein et al. (1995). They found a nonlinear increase in the L5-S1 flexion-extension moment with increased HL with the largest differences occurring when the HL was over 40 cm. When lifting approximately 10 kg, the knee flexion-extension moment at the time of peak L5-S1 moment shifted from extension to flexion between 20 and 40 cm and then remained fairly constant with further increments. It was suggested that the subjects squatted down and shifted their weight posteriorly as the HL increased (Schipplein et al., 1995). In the present study, the subjects also squatted down in addition to flexing their shoulders further in response to the increased HLs (from 30 to 60 cm), while maintaining a similar degree of trunk flexion. This resulted in slight increases in the peak L5-S1 A-C and A-P shear forces, while both the left (from 56 to 84 Nm) and right peak knee flexion-extension moments (from 58 to 85 Nm) increased linearly, similar to the knee resultant JRFs. Increments in the DH had a significant effect on the L5-S1 A-C and A-P shear forces, but it should be noted that the peak forces over the whole lifting cycle mostly occurred at the initiation of the lifts rather than during the placement of the box on the shelves. This was most likely due to the initial lifting height (25 cm from hands to floor) being lower than the shelf heights, as previous research has shown that a higher lifting height may significantly reduce the peak load on the lower back (Lavender et al., 2003;Plamondon et al., 2012). Interestingly, while a higher DH reduced the peak A-C and A-P shear forces, the A-C impulse increased, particularly when the DH was 120 and 150 cm. Hence, while a higher DH may alleviate the acute loads on the lower back to some extent, these benefits may be counteracted by a longer duration of exposure if the DH is 120 cm or above. The DH-30 condition showed the highest knee JRFs, while the lowest forces were found when the DH was 90 cm. It appears that a 90 cm DH was the most optimal for reducing the forces in the knees, as the subjects were able to place the box with the knees fairly extended without having to accelerate the box excessively during pick up to compensate the shoulders. When viewing Fig. 3b and Sup. Fig. 1b, it is reasonably clear that the DH-120 and DH-150 conditions showed the highest knee flexion angles and resultant JRFs at the initiation of the lifts, indicating that the subjects attempted to accelerate the box with their legs more than for the lower DHs. On the other hand, during the DH-30 and DH-60 conditions, the subjects maintained a higher degree of knee flexion throughout the lifting cycle, hereby also maintaining a relatively high load. Finally, the increments in DH had the most substantial influence on the forces in the shoulders, reaching 245 % BW (1942 N) and240 %BW (1903 N) for the left and right side, respectively, when lifting to a shelf height of 150 cm. This incremental increase in peak resultant JRFs is consistent with previous studies employing 2-D static (Pekkarinen and Anttonen, 1988) and 3-D dynamic biomechanical models (Faber et al., 2009) to study the effect of DH on the shoulder flexion-extension and total moment, respectively. Collectively, these results indicate that lifting to near waist height (e.g. DH-90) seems considerably less strenuous compared with both lower and higher DHs based on a simultaneous assessment of the knee, shoulder and L5-S1 JRFs as well as the L5-S1 A-C impulse.
The increments in AA had a negligible effect on the L5-S1 A-C force, but a substantial influence on both the A-P and M-L shear forces. When compared to DH-30, which was initiated from a symmetrical but otherwise identical position, the A-C, A-P and M-L forces were 4.6, 52.7 and 525% higher when the AA was 75 • . The M-L force increased by 6.5 %BW (51.4 N) on average for every increment in AA, reaching a peak across all conditions of 40.0 %BW (317 N) when the AA was 75 • . This can partly be explained by the high degree of trunk flexion, lateral bending and rotation as well as glenohumeral flexion that were associated with the high AAs (see Sup. Fig. 2a and b). In comparison, Marras and Davis (1998) found a 3 and 5% increase in A-C, 4.5 and 8% decrease in A-P and 35 and 58% increase in M-L force when lifting from AAs of 30 and 60 • to the right of their subjects, respectively, compared with a symmetrical lift. Hence, the relative difference in A-C force was very similar between studies, while substantial differences can be seen for the shear forces. These notable differences likely stem from a combination of differences in lifting techniques, measurement tools and biomechanical models; Marras and Davis (1998) used a back electrogoniometer and force plate measurements as input to a dynamic, EMG-assisted biomechanical model (Granata and Marras, 1995). Besides the effect on the L5-S1 shear forces, the AAs also had a significant effect on the knee and shoulder forces. As could be expected, the load on the knees shifted towards the side of the load with increments in the AA, resulting in a substantial difference of 416 %BW (3289 N) between the left and right knee for AA-75. Although this tendency is unavoidable, the fact that the subjects were prohibited from moving or rotating their feet probably exacerbated this effect. Furthermore, the forces in both shoulders increased incrementally with increased AA, reaching forces at AA-75 equal to a symmetrical lift with double the LM (20 kg). Collectively, these results exemplify the negative consequences of performing asymmetrical lifts, as it not only leads to increased shear forces in the lumbar spine, but also significantly higher peak forces in the knees and shoulders.
A common observation across the LM, HL and AA conditions were that the deposit phase showed similar JRFs on average, but the duration of the largest forces were longer, while some of the differences between condition levels were increased. This was particularly visible for the L5-S1 A-C and A-P shear forces during the LM condition, where especially the durations of exposure were longer, while the differences in A-P shear were slightly larger. These overall tendencies were probably a result of the subjects having to place the box as close to its original location as possible, while additionally being encouraged to lift in a controlled fashion in general. In an occupational context, workers will often be able to drop boxes from higher heights instead of carefully placing them, hereby reducing the cumulative load. Another interesting difference between the pick-up and deposit phases were the degree of trunk rotation during the AA condition (See Sup. Fig. 2a). On average, the degree of trunk rotation as well as duration of exposure increased notably during deposit compared with pick up, indicating that carefully placing the box on the wooden bench prolonged the subjects' exposure to unfavorable postures. However, the L5-S1 and knee resultant JRFs were generally slightly lower during deposit, but the exposure to large forces was more prolonged. In combination, these findings could indicate that the subjects accelerated the load during pick up from a slightly more optimal posture, which resulted in higher peak JRFs, but a shorter duration of exposure.
The calculation of the NIOSH lifting index showed that LM was a particularly strong predictor of an increased risk of low back disorders, as the load increments increased the lifting index linearly. The LM-10 condition resulted in a lifting index near the action limit (0.98), which produced an A-C force of 4132 N. This force estimate is considerably higher than the 3400 N A-C tolerance limit proposed by NIOSH (Nelson et al., 1981;Waters et al., 1993), which should correspond to the action limit. However, this discrepancy was expected, as static, sagittal-plane models like those used to develop the NIOSH equation will produce lower estimates of low back loading compared with dynamic models (Garg et al., 1982;McGill and Norman, 1985;de Looze et al., 1994). Similar to the LM conditions, the lifting index also appeared to be a reasonable predictor of low back loading during the HL conditions, but an inverse relationship was found between the lifting indices and the increments in DH as well as the peak L5-S1 A-C force. This can be explained by the travel distance multiplier (D), which increased the lifting index, but the peak A-C forces typically occurred during the pick-up phase with DH-30 being the only exception. However, as the results for the L5-S1 A-C impulse showed a significantly higher cumulative load with every increment in DH, the NIOSH equation is not necessarily inappropriate for evaluating the risk of injury to the spine during these types of conditions. Finally, the NIOSH equation did produce higher lifting indices with increments in AA, but the only level that resulted in a lifting index equal to the action limit was AA-75 (1.00), which corresponded to an A-C force of 4295 N. Previous research has indicated that the lifting equation may produce RWLs resulting in L5-S1 A-C and A-P shear forces exceeding their proposed tolerance limits of 3400 and 1000 N (McGill et al., 1998;Gallagher and Marras, 2012), respectively, during lifts involving a high degree of asymmetry and trunk flexion (Behjati and Arjmand, 2019). In the present study, the A-P shear forces did not exceed 1000 N for any of the tested conditions, but large M-L shear forces were found for the highest AAs. As no tolerance limit for M-L shear have been proposed in the literature, it is difficult to evaluate the injury risk associated with these forces, but one could argue that the NIOSH equation does not adequately adjust for increased M-L shear forces based on the results of this study. Therefore, caution should be taken when applying the equation during highly asymmetric lifts involving high degrees of trunk flexion, which has also been suggested in previous research (Behjati and Arjmand, 2019).
The study contains a few limitations that should be noted. First, as the subjects had to stand still on the force plates, it restricted their ability to change technique in response to the changing conditions. This included rotating their feet to account for the increments in asymmetry, which could potentially have reduced the joint loads. Second, the use of marker-based motion analysis as well as additional measurement equipment (see supplementary material) may further have restricted their natural execution of the lifts. Third, instructing the subjects to lift in a controlled fashion might also have affected the procedures, as for instance, by incentivizing the subjects to not accelerate or drop the box during the pick-up and deposit phase, respectively. In an occupational setting, workers may use their legs to accelerate loads during pick up at low initial lifting heights to a larger extent or drop the boxes instead of placing them to avoid the additional effort. Fourth, the subjects were barefoot during the experiment, which does not mimic real-life working conditions, where workers often wear protective footwear. Hence, it would have been more representative of real-life working conditions if the subjects had worn protective footwear as well as decreased the risk of injury if they had dropped the box on their feet. We did, however, fit the box with protective foam padding. Fifth, as the forces between the hands and box were not measured (as mentioned in section 2.2), we used contact elements to estimate these forces, similar to previous studies (Larsen et al., 2020;Koblauch, 2016;Skals et al., 2021;Muller and Corbeil, 2020). However, the accuracy of the estimated external hand forces and moments is unclear. Finally, including a convenience sample may be problematic in relation to the generalization of the results. Taking individual factors into accounte.g. age and individual physical capacitycould be considered in future studies.

Conclusions
This study provided a detailed analysis of the dynamic peak forces in the major joints during MMH with varying loads, asymmetry, start and end locations based on state-of-the-art musculoskeletal models. The results mostly support the findings of previous research employing dynamic biomechanical models, showing that LM and AA had a substantial effect on the L5-S1 JRFs, the DH had the most pronounced effect on the shoulder JRFs, while the forces in the knees were significantly affected by all the lifting factors. However, notable differences between studies were identified in regards to the magnitude and trends of the estimated forces, most probably due to differences in lifting techniques and biomechanical models. Finally, calculations of RWLs and lifting indices using the NIOSH lifting equation showed that this method was mostly appropriate during lifts with varying LM, HL and DH, but did not adequately adjust for increased M-L shear forces during highly asymmetrical lifts. Our findings provide a detailed insight into the influence of ergonomic lifting factors on the peak dynamic loading of the major joints and can be used to evaluate the injury risk to multiple joints simultaneously in response to a wide range of lifting conditions.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Mark de Zee is co-founder of the company AnyBody Technology A/S that owns and sells the AnyBody Modelling System, which was used for the simulations. Mark de Zee is also a minority shareholder in the company.