Manual material handling in the supermarket sector. Part 2 Knee, spine and shoulder joint reaction forces

Manual material handling is common in supermarkets and may be a contributing factor to the high prevalence of work-related musculoskeletal disorders, particularly to the lower back. This cross-sectional study applied state-of-the-art musculoskeletal models driven by kinematic data obtained in two supermarkets to estimate joint re- action forces in the knees, shoulders and lumbar spine under dynamic lifting conditions. Based on 1479 lifts from 15 workers, 8 tasks for which the compression or shear forces in the L5-S1 joint exceeded well-known biome- chanical tolerance limits were identified. High shoulder forces were associated with lifting relatively heavy merchandise to high shelves, while the weight of the handled merchandise was the main predictor of high knee forces. The study addressed well-known limitations associated with traditional lifting analysis tools and was the first to present a detailed analysis of the biomechanical loads during manual material handling tasks in the supermarket sector based on field measurements.


Introduction
Low back pain is the most common musculoskeletal disorder worldwide and the single largest cause of years lived with disability (Vos et al., 2016) with a substantial proportion of global cases attributable to occupational exposures (Punnett et al., 2005). The most well-documented occupational risk factor for low back pain is manual material handling (MMH) (Burdorf and Sorock, 1997;NIOSH, 1998;Hoogendoorn et al., 1999;da Costa and Vieira, 2010;Coenen et al., 2014;Lötters et al., 2003;Punnett and Wegman, 2004;Fernandes et al., 2016;Cole and Grimshaw, 2003;Swain et al., 2020;Kuiper et al., 1999); a term that typically includes both lifting, lowering, carrying, pushing and pulling. For this reason, the physical demands of MMH has been studied extensively using biomechanical, psychophysical and physiological approaches (Dempsey, 1998). Task-based analysis of MMH range from self-reports and observational methods to direct measurements (e. g. motion analysis or electromyography) and load estimations based on biomechanical models (Rajaee et al., 2015;Russell et al., 2007;David, 2005;Takala et al., 2010). The choice of assessment method is typically a trade-off between accuracy, flexibility and cost, as ergonomists are often expected to provide solutions to injury problems in a way that require minimal capital investment (Dempsey and Mathiassen, 2006).
The biomechanical approach has traditionally incorporated estimates of the joint moments and compressive forces in the lower back, typically at the L4-L5 and L5-S1 joints, to assess the risk of low back pain and injury during MMH (Dempsey, 1998). One of the most famous and widely used tools for this type of analysis is The Revised NIOSH Lifting Equation (Waters et al., 1993), which provides a risk estimate based on a set of multipliers, such as load weight, vertical and horizontal distance. Although this represents one of the most comprehensive and substantiated tools available, it has some fundamental limitations: firstly, it is based on a relatively simplistic 2D static biomechanical model for estimating the compressive force at the L5-S1 joint. It is well-known that these types of models underestimate the forces on the spine considerably compared with dynamic models (Gagnon and Smyth, 1992;Garg et al., 1982;Menzer, 2005), as the acceleration components of the body segments are not considered. Secondly, it can only be used in idealized lifting conditions, i.e. two-handed lifting in the sagittal plane with both feet on the ground, and has been shown to be insufficient for analyzing asymmetric lifting (Arjmand et al., 2012;Behjati and Arjmand, 2019). Thus, current lifting recommendations suffer from these limitations and may therefore, not provide the best protection for the workers.
One approach to overcome these limitations has involved the development of detailed, dynamic 3D biomechanical models (Kingma et al., 1996;De Zee et al., 2007;Han et al., 2012;Beaucage-Gauvreau et al., 2019). For example, de Zee et al. (2007) and Han et al. (2012) developed a 3D dynamic lumbar spine model with detailed representations of muscles, ligaments and intra-abdominal pressure. This model has shown higher accuracy than other quantitative tools for computing forces in the lumbar spine when compared with in vivo data (Rajaee et al., 2015;Bassani et al., 2017), but has only been used scarcely for the analysis of MMH (Behjati and Arjmand, 2019;Stambolian et al., 2016;Kobluach, 2016). Furthermore, as this spine model is integrated in a full-body musculoskeletal model, providing anatomically detailed representations of the other involved joints (e.g. the shoulders and knees), it enables a more complete analysis of the internal loads that workers are subjected to. For example, it can provide a better understanding of how loads are transferred throughout the body and how a reduction of forces in one joint may lead to higher forces in another.
One of the main drawbacks of applying 3D dynamic models is that they require detailed kinematics and external forces (e.g. ground reaction forces) as input data, which has traditionally been obtained using marker-based motion analysis and force plate measurements in a laboratory environment. However, recent advancements in inertial-based motion capture technology (Paulich et al., 2018;Schepers et al., 2018) and the availability of methods to predict ground reaction forces and moments from segment kinematics and dynamical properties only (Skals et al., 2017a;Fluit et al., 2014;Karatsidis et al., 2017) has provided new opportunities for obtaining these input data in the field. This combination of methods has recently been shown to produce reasonable estimates of the compression forces in the lumbar spine during lifting compared with the traditional laboratory-based approach (Larsen et al., 2020), providing a powerful new tool to better understand the dynamic loading of the lower back during MMH in real-life scenarios.
One of the industries with the highest prevalence of MMH is the retail industry (Heran-Le Roy et al., 1999); an umbrella category that includes the supermarket or grocery sector among others. Unsurprisingly, the supermarket sector also has a high prevalence of low back pain compared with other major industries (Guo et al., 1999). Furthermore, epidemiological research has shown that musculoskeletal disorders are highly prevalent in the supermarket sector in general (Ryan, 1989;Forcier et al., 2008;Anton and Weeks, 2016) with MMH being the most commonly identified occupational risk factor (Campany and Personick, 1992;Kraus et al., 1997;Clarke, 2003;Violante et al., 2005;Rahman and Zuhaidi, 2017). For example, supermarket workers also report a high prevalence of shoulder, neck and knee pain (Forcier et al., 2008;Anton and Weeks, 2016).
In view of the above, we conducted a two-part analysis of the working postures, muscular demands and biomechanical loads associated with a wide range of common MMH tasks in the supermarket sector based on the newest available technology for full-body inertial-based motion analysis and musculoskeletal modelling, as well as electromyographic measurements. The first paper (Skals et al., 2021) reports joint angle and muscle activity data, while the present paper reports the kinetic results from the musculoskeletal modelling. In the present paper, we applied the aforementioned methods for inertial motion capture (IMC) and ground reaction force prediction to estimate the dynamic loading of the knees, shoulders and lumbar spine using a highly detailed, 3D dynamic biomechanical model based on field measurements from two supermarkets. The aim of this approach was to assess the risk of injury during common MMH tasks in the supermarket sector as well as to identify current limitations of the methodology and obstacles for implementing this kind of technology for field-based analysis of MMH on a larger scale. The data presented here is the first of its kind and will help improve our understanding of the dynamic loading of the body during MMH, hereby, contributing to the ongoing efforts to reduce the incidence of work-related musculoskeletal disorders, particularly the incidence of low back pain and injury.

Methods
The experimental procedures of this two-part cross-sectional study were described in detail in part 1 (Skals et al., 2021), including additional information regarding the subjects, work characteristics, data collection in the two supermarkets as well as a detailed description of the IMC and electromyographic measurements. In the following, these sections (except the electromyographic procedures) are briefly summarized, while the musculoskeletal modelling procedures are described in detail.

Subjects
Seventeen healthy full-time (min. 30 h/week) workers from two stores of a major Danish supermarket company were recruited for the study (8 male and 9 female, age: 27 ± 8 years, height: 174.4 ± 9.1 cm, weight: 76.6 ± 14.7 kg, experience: 9.4 ± 7.4 years, h/week: 37.6 ± 6.3). The study followed the guidelines of The North Denmark Region Committee on Health Research Ethics and all subjects provided written informed consent. Data were collected in October 2018.

Data collection
The work performed by the supermarket workers consisted almost exclusively of receiving, stocking and re-arranging merchandise. The MMH tasks in the stores could roughly be subdivided into the following categories: 1) fruit and vegetables (FV), 2) bread (BR), 3) meat and dairy (MD) and 4) colonial (CO), i.e. edible and inedible goods with long shelf lives. The workers would typically stock goods in all the aforementioned categories on any given day.
Based on observations and conversations with industry stakeholders, a subset of MMH tasks were selected for further analysis. As most of the merchandise could be placed on several shelf heights (H) and from different starting positions (SP), multiple scenarios with different start and end positions were included in the analysis, indicated with the numbers 1-4 (low to high). During all tasks, the subjects were informed about the start and end position, and were asked to lift the merchandise by the handles with both hands. In total, 50 different MMH tasks were included in the analysis with the subjects performing four consecutive repetitions of each task. However, due to a number of issues related to the musculoskeletal modelling, 24 of the 50 tasks had to be excluded from this part of the analysis (described in section 3). The characteristics of the 26 MMH tasks that were successfully modelled are described in Table 1. The MMH categories as well as the tasks within each category were counterbalanced to avoid any order effects.
Bodyweight and body dimensions were measured with a scale and caliper, respectively, before the measurement equipment was instrumented on the subjects. Subsequently, two investigators (SS and RB) accompanied the subjects in to the shopping area with all the merchandise assembled in a transport cage with two shelves (SP1/Low and SP2/High) and the measurement systems, a laptop computer and video camera positioned on a rolling table (see supplementary material). The transport cage was the typical assistive device used during stocking with the selected shelf heights closely corresponding to the top and bottom of a filled pallet. When initiating a measurement, the investigators placed the transport cage next to the shelf and asked the subject to stand with their left side to the cage, while the rolling table with the measurement equipment was placed opposite to the subject. From here, the subject would lift the specific merchandise from the transport cage onto the appropriate shelf to their right side, whereafter one of the investigators returned the merchandise to the starting position. There were no shelves or other merchandise obstructing the subjects when performing the lifts. This procedure was repeated four times before initiating the next series of four lifts from a different starting position and/or to a different shelf height. When all combinations of lifts for the specific merchandise had been performed, the merchandise was exchanged to another product within the same food category and the procedure repeated. When all tasks in a specific food category had been performed, the investigators moved the transport cage and rolling table to the next area of the store and completed a calibration of the IMC system before the next series of measurements were performed.

Instrumentation
Full-body kinematics were measured using the Xsens MVN Awinda wireless motion-tracker (Xsens Technologies BV, Enschede, The Netherlands). The IMC system consists of 17 inertial measurement units (IMUs), sampling at 60 Hz, which were attached to the subjects with the accompanying velcro straps, headband and a tight-fitting customized tshirt. Each IMU contains a 3D accelerometer, gyroscope and magnetometer, and transmit data in real-time to a master receiver (Paulich et al., 2018). The IMU data is used to drive a 23 segment kinematic model, which is scaled to the subjects based on the manually measured body dimensions. The collected data were reprocessed using the embedded tool in Xsens MVN Analyze v. 2019.0.0 and subsequently exported as Biovision Hierarchy files (BVH-files).

Musculoskeletal modelling
Musculoskeletal models were developed in the AnyBody Modeling System v. 7.2 (AMS) (Anybody Technology A/S, Aalborg, Denmark) based on the BVH_Xsens template from the AnyBody Managed Model Repository v. 2.2.3. The model template included a lower extremity model based on the cadaver study of Carbone et al. (2015), a lumbar spine model based on the work of Hansen et al. (2006), de Zee et al. (2007 and Han et al. (2012), as well as shoulder and arm model based on the work of Van Der Helm et al. (1992) and Veeger et al. (1991Veeger et al. ( , 1997. The lumbar spine was modelled with non-linear disc stiffness and included 188 muscle elements with the representation of ligaments and intra-abdominal pressure similar to Han et al. (2012). Body segments were represented as rigid bodies with the spine model consisting of 14 rigid segments representing the cervical vertebras, thorax, and lumbar vertebras as well as the sacrum (De Zee et al., 2007;Hansen et al., 2006). The model followed a spine rhythm that distributes the trunk motion over the vertebral bodies using a coupled-mechanism (Hansen et al., 2006). Each shoulder and arm model had 140 muscle elements, including 12 for deltoideus, 6 for subscapularis, 6 for infraspinatus, 6 for supraspinatus, 6 for teres minor and 10 for pectoralis major. Muscles with complex wrapping behavior, e.g. deltoideus, have insertion points on the bones and wrap over analytical surfaces, e.g. ellipsoids or spheres. The knee was modelled as a hinge joint with a fixed rotation center and axis, and the patellar tendon defined as a non-deformable element connecting the patella to the tibia. Each leg in the lower extremity model include 169 muscle elements with several spanning the knee joint, e.g. gastrocnemius and biceps femoris. Further details about the model architecture can be found in the online documentation of the Table 1 Abbreviations and description of the selected manual material handling tasks, including the merchandise handled, the starting position, shelf height, weight and dimensions (length (L) x width (W) x height (H)). Note that there are two listed shelf heights, corresponding to the first and second store where data were collected.

Abbreviation (tables)
Abbreviation ( model repository (Lund et al., 2020). The models had a total of 80 degrees-of-freedom (DOF), including 2 x 3 DOF at the ankle joints, 2 x 1 DOF at the knee joints, 2 x 3 DOF at the hip joints, 6 DOF at the pelvis, 3 DOF between pelvis and thorax, 10 x 3 DOF for the lumbar and cervical vertebras, 2 x 2 DOF at the elbow joints, 2 x 3 at the glenohumeral joints, 2 x 3 DOF at the sternoclavicular joints, 2 x 2 DOF at the wrist joints, 1 DOF at the neck joint and 6 DOF between the hands and box. Examples of the models during 3 MMH tasks are illustrated in Fig. 1. 3D computer-aided design models of the handled merchandise were created in SolidWorks v. 2017.5.0 (Dassault Systems, Vélizy-Villacoublay Cedex, France). The geometry, mass and inertial properties of the models were based on measurements made during data collection (see Table 1). For all the merchandise, a cuboid with the material properties of the merchandise was modelled and placed within a shell with the properties of cardboard. For the milk crates however, the 15 milks in the crates were modelled individually in a similar fashion, and placed in a lidless shell with the material properties of high-density polyethylene.

Model scaling and kinematics
The musculoskeletal models were scaled according to the manually measured segment dimensions, which were inputted to the IMC software before initiating the measurements. The BVH-files that are exported from the IMC software uses a hierarchical structure that contain a description of the kinematic model in a static pose, and the absolute position and orientation of the root pelvis segment as well as the joint angles between the segments for each time frame (Karatsidis et al., 2019). To enable scaling and marker tracking between the stick figure and the musculoskeletal model, we used the framework first presented in Skals et al. (2017b) and specifically applied for processing IMC data in Karatsidis et al. (2019). This framework was added as the standard tool for processing BVH-files in AMS v. 7.2.3.
In short, the segments of the musculoskeletal model were scaled according to the joint-to-joint distances of the stick figure. However, as the dimensions of the pelvis, foot and trunk were not directly scalable from the stick figure, additional nodes were added to the unscaled segments of the musculoskeletal model at points that were identifiable on the stick figure. The distance between the added nodes were then computed, whereafter the ratio between the unscaled segments and the distance between the added nodes were multiplied onto the segment lengths of the stick figure and then applied to scale the musculoskeletal model. Lastly, virtual markers were introduced on both the stick figure and the musculoskeletal model to enable marker tracking. The kinematics were solved using a non-linear least-square optimization problem, which minimizes the least-square difference between the virtual markers of the two models (Andersen et al., 2010). The mass of the body segments were determined by distributing the body mass of the subjects according to the regression equations presented in Winter et al. (2009).
As the motion of the handled merchandise was not measured, the kinematics were solved by creating spherical joints between the hands and boxes. The joint constraints were placed at the handles on either side of the boxes, where the subjects had been instructed to place and keep their hands when performing the tasks. In addition, a kinematic constraint was created to control the box rotation in the sagittal plane by adding a point at the proximal and distal end on the right side of the boxes, which had to remain at the same height in relation to the ground plane. This constraint ensured that the boxes would remain horizontal in the sagittal plane, while the other DOF were controlled by the movement of the hands.

Prediction of external forces
The ground reaction forces and moments were predicted using a method first presented by Fluit et al. (2014) and further developed and evaluated by Skals et al. (2017a) and Karatsidis et al. (2019), which has shown comparable accuracy to force plate measurements during activities of daily living (Fluit et al., 2014), sports-related movements (Skals et al., 2017a) and IMC-based gait analysis (Karatsidis et al., 2019). The method employs a set of 25 dynamic contact elements that were defined under each foot of the musculoskeletal models, each consisting of five uniaxial force actuators able to generate a positive normal force, as well as positive and negative static friction forces in the anteroposterior and mediolateral directions (with a friction coefficient of 0.8). A non-linear strength function ensures that the contact elements only generated forces when they were close to the ground and almost stationary (Skals et al., 2017a). The activation threshold distance and velocity in relation to the ground plane were set to 35 mm and 0.8 m/s, respectively. Lastly, small residual forces and moments were placed at the pelvis with a strength of 10 N and 10 Nm to improve numerical stability of the prediction algorithm. The actuation of each contact element was computed as part of the muscle recruitment (see section 2.3.3).
Twelve muscle-like contact elements with a high strength (400,000 N and Nm) were also defined at the joints between the hands and boxes to estimate the external forces and moments of the handled merchandise, hereby creating actuators for the added DOF and enabling the inverse dynamic analysis to be solved.

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

Data analysis
L5-S1 peak compression, anteroposterior shear and mediolateral shear forces, as well as the peak resultant JRFs at the left and right knee and glenohumeral joints were extracted from the musculoskeletal models and normalized to percentage of body weight (%BW). In addition, the peak compression and anteroposterior shear forces at the L5-S1 joint were compared with the widely used compression and anteroposterior shear tolerance limits of 3400 N (Waters et al., 1993;Nelson et al., 1981) and 1000 N (McGill et al., 1998;Gallagher and Marras, 2012), respectively, in order to evaluate whether the analyzed MMH tasks posed a risk for developing low back disorders. The start and end points of the lifting cycles were determined as the instant where the subjects lifted the handled merchandise off the ground with a secure grip and the instant just before the box was securely placed on the shelf, respectively. The data were resampled to 101 data points (one lifting cycle) for illustrative purposes.

Statistical analysis
Repeated measures linear mixed models (Proc Mixed, SAS) were used to test if any differences existed between the MMH tasks for the selected variables. The peak JRFs were the dependent variables, while the MMH tasks were included as fixed effects and the subjects as random effects to take repeated measures into account. The main purpose of the linear mixed model analyses were to determine least square means with 95% confidence intervals for each MMH task in order to rank the tasks from highest to lowest for each outcome variable to facilitate the risk assessment. Residual diagnostics plots were visually inspected to ensure a normal distribution of the residuals as well as homogeneity of variance, while within-subject correlation was assumed. The covariance structure was set to Variance Components, while the model was fit using restricted maximum likelihood estimation (REML). The results are presented as least square means with 95% confidence intervals based on a Satterthwaite approximation. All statistical analyses were performed in SAS v. 9.4 (SAS Institute Inc., Cary, NC, USA).

Results
Of the 17 subjects from which data were collected, 15 were included in the final analysis (see Part 1 for further details). From these 15 subjects, 78 trials were either missing or excluded: 48 trials were missing as the tasks FV2 and FV3 were not placed on both shelves for subject 2 and 3, while FV3 was not lifted from both starting positions for subject 12 and 17. Five additional trials of the tasks BR1 and MD2 were missing for subject 2. Sixteen trials were excluded for subject 17 as the wrong hand was used during two of the one-handed lifts (the one-handed tasks were eventually excluded as described below), 4 trials excluded due to errors in the kinematic data, while 5 trials were excluded due to extreme values attributed to errors in the system's calculation of a joint angle. A total of 2922 IMC trials containing all 50 MMH tasks were initially used to drive the musculoskeletal models. However, there were a number of issues during the modelling procedures that resulted in the exclusion of a large proportion of the analyzed trials: firstly, as the handled merchandise followed the movement of the hands, inaccuracies in the IMC systems determination of the hand placement led to large errors in the orientation of the boxes in some cases. This error occurred almost exclusively for the MMH tasks where narrow, small boxes were handled, but to an extend that necessitated the exclusion of all trials related to the handling of cold cuts, yoghurts and fresh herbs. Secondly, during a large proportion of the one-handed lifts as well as some individual trials for the other tasks, a high degree of palmar dorsiflexion lead to errors in the muscle wrapping of one or more wrist flexors, which caused the inverse dynamic analysis to fail. Correcting the muscle wrapping for each individual trial was a time-consuming and complicated task, and we, therefore, decided to exclude all the one-handed lifts, specifically the handling of minced beef, single yoghurts and vegetable oil. Seventy-five additional trials of the other MMH tasks were also excluded due to this error. Finally, 6 trials were excluded due to random errors occurring when attempting to solve the inverse dynamics problem, which did not have any immediate solution. In total, 1479 trials of the remaining 26 MMH tasks (see Table 1) were included in the final analysis, comprising all the relatively heavy, two-handed lifts.
The results of the linear mixed model analyses showed significant differences for the main effect (i.e. the MMH tasks) for each outcome variable (p < 0.001). The least square means with 95% confidence intervals for the L5-S1 JRFs are listed in Table 2, and the results for the knee and glenohumeral resultant JRFs listed in Table 3. The 26 MMH tasks are ranked from highest to lowest for each selected variable. Fig. 2 illustrate the JRFs over the complete lifting cycles for the tasks FV1-SP1-H1/Bananas-LowToLow, FV3-SP1-H3/Cucumbers-Low-ToHigh and MD3-SP1-H3/Milk-LowToHigh, respectively, while additional figures for the other analyzed tasks can be found in a supplementary database .

L5-S1 JRFs
The MMH tasks exhibiting the highest compression forces at the L5-S1 joint were the handling of bananas (553 and 539 %BW) and milk (from 424 to 506 %BW). Furthermore, the handling of cucumbers (from 413 to 449 %BW) and bread placed on the lowest shelf (422 and 425 % BW) also showed relatively high compression forces. These overall results were more or less mirrored for the L5-S1 anteroposterior shear force, where the handling of bananas (142 and 155 %BW) and milk lifted from or to a low position (from 134 to 144 %BW) showed the highest values, followed by the handling of cucumbers and bread, particular when lifted from or to a low position. The L5-S1 mediolateral shear forces were in general very low compared with the compression and anteroposterior shear forces, ranging from 5 to 16 %BW across all the analyzed tasks.

Knee resultant JRFs
The tasks showing the highest knee resultant JRFs were overall different for the left and right leg with the handling of bananas (761 and 848 %BW) and cucumbers (from 681 to 751 %BW) showing the highest forces in the left knee, while the handling of milk showed the highest values for the right knee (from 744 to 799 %BW). Overall, the data exhibited a trend towards the heaviest merchandise resulting in the highest knee forces, while there was no clear trend associated with either the starting positon or shelf height.

Glenohumeral resultant JRFs
For the glenohumeral resultant JRFs, the highest ranked tasks for both the left and right side were Cucumbers-LowToHigh (225 and 227 % BW), Bread-HighToHigh (223 and 225 %BW), Cucumbers-HighToHigh (222 and 220 %BW) and Bread-LowToHigh (217 %BW for both sides). Furthermore, the handling of bananas also showed relatively high glenohumeral forces (from 140 to 167 %BW), although these forces were significantly lower. Hence, the data showed that the heaviest merchandise lifted to the highest shelf heights resulted in the highest forces in the glenohumeral joint.

Table 2
Peak axial compression, mediolateral shear and anteroposterior shear forces presented in percentage of bodyweight (%BW) with 95% confidence intervals. The 26 manual material handling tasks are ranked from highest to lowest. L5-S1 axial compression force L5-S1 mediolateral force L5-S1 anteroposterior force

Biomechanical tolerance limits
When viewing the data in relation to the biomechanical tolerance limits of 3400 N for compression and 1000 N for anteroposterior shear force, there were 8 and 6 tasks where the upper confidence limit exceeded the compression and shear tolerance limits, respectively, namely Bananas-LowToLow (4188 N for compression and 1191 N for shear), Bananas-HighToLow (4088 N and 1097 N), Milk-LowToLow (3854 N and1113 N), Milk-LowToMid (3811 N and1069 N), Milk-LowToHigh (3775 N and1062 N), Milk-HighToLow (3742 N and1041 N), Milk-HighToMid (3577 N) and . Several other tasks showed spinal forces that were near the tolerance limits, particularly Cucumbers-LowToMid (3398 N and 978 N).

Discussion
We applied IMC-based kinematic data obtained in two supermarkets to drive musculoskeletal models and estimate the biomechanical loads during common MMH tasks. Based on the dynamic loading of the L5-S1 joint, we identified 8 and 6 MMH tasks that exceeded well-known biomechanical tolerance limits for compression and anteroposterior shear force, respectively, and hence, could pose a risk for the development of low back pain and injury. Furthermore, relatively high glenohumeral resultant JRFs were found when heavier merchandise were moved to high shelf heights, while the forces in the knees were highest for the heaviest merchandise, i.e. handling of bananas and milk, regardless of start and end position.
The main finding of the study was the identification of several MMH tasks that posed a risk for developing back pain and injury, namely the handling of bananas as well as milk crates lifted from or to a low position (13.5-15 cm above the floor). These tasks were the heaviest merchandise (20.2 and 17.3 kg, respectively) handled on a daily basis in the two supermarkets where data were obtained. Furthermore, the handling of cucumbers and bread, particularly from a low to a high position, also resulted in spinal loads that were close to the biomechanical tolerance limits, despite these merchandise having a much lower weight (10.2 and 7.9 kg, respectively). Hence, while the weight of the handled merchandise was a strong predictor of the compressive forces in the lumbar spine, the height of the starting positon and shelf height also influenced these forces considerably. This was also evident from the instant in the lifting cycle where peak compression occurred, which was either at the initiation of the lift or when placing the merchandise on the shelf (see Fig. 2 and supplementary database for further details).
For the knee resultant JRFs, the merchandise with the highest weight resulted in the highest knee forces regardless of start and end position. This was not surprising, as the resultant forces in the knees are highly dependent on the total mass during a slow, controlled movement. Furthermore, during the handling of milk, which showed the highest forces in the right knee, the subjects had to turn and walk one step to place the merchandise on the shelf. Hence, the peak knee forces during this task occurred when the subjects were standing with the merchandise on one leg at the end of the lifting cycle. For the handling of bananas, which showed the highest forces in the left knee, the peak force occurred during the initiation of the lift when the subjects were squatting down to pick up the merchandise. It is difficult to compare these forces to previous studies, as 3D dynamic joint reaction forces in the knees have not previously been reported during manual handling. A study by Skals et al. (2017a), which reported knee resultant JRFs for sports-related movements using a similar musculoskeletal model, may provide some context: for example, the handling of bananas resulted in knee forces of around 8 times the subjects' bodyweights, which was equivalent to performing counter-movement jumps and running (Skals et al., 2017a). However, knee flexion-extension moments based on dynamic biomechanical models has been estimated in previous research on lifting, as for instance, by De Looze et al. (1993). In this study, the authors found peak knee flexion-extension moments of approximately 170 and 35 Nm during squat-and straight leg lifting of a 15.3 kg barbell from a starting position of 10% of the subjects' body height, respectively. In the present study, the peak knee flexion-extension moments ranged from 54 (Bread-MidToMid) to 105 Nm (Bananas-MidToLow) with a median of Fig. 2. Time-series curves of the L5-S1 mediolateral (top left), axial compression (top center) and anteroposterior (top right) joint reaction forces, as well as the glenohumeral (middle) and knee (bottom) resultant joint reaction forces for the right (R) and left (L) side during the tasks FV1-SP1-H1/Bananas-Low-ToLow (red), MD3-SP1-H3/Milk-LowToHigh (blue) and FV3-SP1-H3/Cucumbers-LowToHigh (green). The forces are presented as the mean (solid line) with one standard deviation (shaded area). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) 80 Nm .
The glenohumeral resultant JRFs were by far the highest when cucumbers and bread were placed at the highest shelf height (from 217 to 227 %BW). These forces were approximately 36% higher than the second highest ranked task, Bananas-HighToLow (160 and 167 %BW). Hence, it seems reasonable to recommend that these boxes should only be placed on shelves below shoulder height or that the merchandise is stocked individually. A study by Faber et al. (2009) used the 3D electromyographic-assisted trunk model of Kingma et al. (1996) to estimate the net total joint moment in the right shoulder during masonry work. Peak moments ranged from approximately 30 to 50 Nm when lifting building blocks of 6-16 kg from a table to shoulder level with two hands. In the present study, the peak total joint moments in the shoulders ranged from 20 (TomatoCans-HighToLow) to 73 Nm (Cucumbers-LowToHigh) with a median of 45 Nm (see supplementary database for further details )). In comparison, the median total joint moment of 45 Nm for the analyzed MMH tasks roughly corresponded to lifting a 14 kg building block to shoulder level with two hands, meaning that the handling of most variations of milk (from 40 to 68 Nm), bananas (from 58 to 70 Nm), cucumbers (from 46 to 73 Nm) and bread (25-69 Nm) shows higher total moments in the shoulders.
We chose to use the well-known spinal compression tolerance limit of 3400 N proposed by The National Institute of Occupational Safety and Health (Waters et al., 1993;Nelson et al., 1981) as well as the shear tolerance limit of 1000 N recommended by McGill et al. (1998) and Gallagher and Marras (2012). The compression limit is based on the results of field studies that provided some quantitative data linking compressive forces with the incidence of low back disorders, and represent the maximum load acceptable for the majority of workers for infrequent exertions (Waters et al., 1993). However, this limit should be used with caution, as it is highly generalized and do not control for important individual factors, such as the workers age, sex and bodyweight. In a review of in vitro studies of the compression tolerance of spinal units, Genaidy et al. (1993) proposed a regression equation to calculate the maximum permissible compression forces during MMH for different percentiles of the population, which adjusts for sex and age. The results for the 10th percentile in the 20-29 age group shows large differences between men (4480 N) and women (3431 N), and a substantial influence of age: for instance, the maximum permissible compressive forces for women in the ages 20-29, 30-39, 40-49 and 50-59 were 3431 N, 2572 N, 1712 N and 853 N, respectively. These results exemplify the importance of considering individual factors in relation to job demands as well as carefully planning the work to accommodate the inherent loss of physical capacity (e.g. muscle strength) with age. Another important factor to consider when assessing the injury risk to the spine during lifting is the influence of spine posture in relation to the fatigue failure tolerance during compressive loading. For instance, an in vitro study by Gallagher et al. (2005) found a reduction in number of compressive load cycles to failure for lumbar spinal units of 61% and 97% when the specimens were positioned to replicate 22.5 • and 45 • torso flexion, respectively, compared with a neutral posture. Hence, frequent lifting that require high degrees of torso flexion may substantially reduce the spinal structures tolerance to compression forces.
The use of 3D dynamic biomechanical models to estimate L5-S1 compression forces during lifting has traditionally been based on input data from laboratory settings. For example, Marras and Davis (1998) found mean peak L5-S1 axial compression forces of approximately 3700 N when a 13.6 kg load was lifted from knee height and a horizontal distance of 53.3 cm by male subjects. When the load was placed to the left of the subjects with an asymmetry angle of 60 • , the compression force increased to around 4200 N. If we compare these data to the estimates for the seven male subjects in our sample, the symmetrical lift corresponded roughly to the seventh ranked task, Milk-HighToMid (3693 N), while the asymmetrical lift corresponded to the second ranked tasks, Bananas-HighToLow (4252 N). The highest compression force found for the male subjects in this study was for .
When evaluating the results in relation to the shear tolerance limits of 1000 N, it should be noted that this threshold is for infrequent loading (<100 loadings per day), while a limit of 700 N has been proposed for frequent loading (100-1000 loadings per day) (McGill et al., 1998). This is important as previous research has estimated that Danish supermarket workers handle 1212 kg on average during a workday, and that an exposure-response relationship exists between their cumulative workload and increased low back pain intensity . In addition, the evaluation of the methodological framework used to estimate these forces in the present study showed that the anteroposterior shear force was underestimated due to an underestimation of the trunk flexion angle (Larsen et al., 2020). Therefore, other tasks may also have actually exceeded the tolerance limit, particularly the handling of bread and cucumber to or from a low position.
In conclusion, we recommend that the handling of bananas, milk, cucumbers and bread in the participating supermarket company should be reconsidered in order to reduce the risk of injury to the workers. More specifically, full boxes of bananas and milk crates should optimally not be handled manually as they produce spinal loads that are considered hazardous based on the best available criteria. Instead, individual stocking could be considered, but limiting the amount of merchandise in each box or the size of the boxes may be warranted. In addition, boxes of cucumbers and bread should preferably not be lifted to the highest shelf heights (108-146.5 cm), as these tasks produce substantially higher shoulder JRFs compared with the other analyzed merchandise. Finally, our results indicate that deep squatting or single leg stances (e.g. during walking) when handling the relatively heavy merchandise (i.e. bananas, milk and cucumbers) should be minimized to reduce the peak forces in the knees.

Limitations and strengths
Besides the possible underestimation of the anteroposterior shear force, the study had a number of other limitations. First, as described in section 3, there were a number of issues with the musculoskeletal modelling procedures, which were mostly related to errors in the kinematic data. Specifically, the errors in the IMC systems estimation of the hand placement lead to errors in the orientation of all the smaller, narrow boxes, leading to the exclusion of a large proportion of the trials. Furthermore, the kinematic data showed excessive palmar dorsiflexion for multiple subjects and trials, particularly during one-handed lifting, which led to errors in the muscle wrapping of the wrist flexors. These issues may be overcome over time, as these novel methods are enhanced due to developments of the hardware (Xsens) and software (AMS) in response to limitations being identified. However, in its present state of development, the methodology was not applicable for the analysis of smaller, narrow boxes and one-handed lifts, but proved effective for the analysis of the relatively heavy merchandise in larger boxes. Second, as the position of the hands in relation to the box has to be accurately defined when performing the inverse dynamic analysis, we instructed the subjects to lift with both hands using the handles on either side of the merchandise. In doing so, we imposed a lifting technique on the workers and did therefore, not accurately capture the natural intra-and intersubject variability. In its present state of development, the methodology requires standardizing the procedures in the field to some degree, which limits its applicability. Third, a well-known issue when using IMUs is the effect of magnetic distortions from the surrounding environment, causing orientation and positional drift over the course of a measurement. To what extent these distortions influenced the data is hard to tell, but we were well-aware of this issue during data collection and performed frequent calibrations of the system to minimize its influence. Furthermore, significant efforts have been made by the developers of the IMC system to continuously correct for drift using an advanced Kalman filter as well as a post-processing tool (HD re-process) in the accompanying software that improves the consistency and precision of the kinematics and global position estimates. The overall function and effect of the Kalman filter was described in Paulich et al. (2018). Fourth, the methodology is very time-consuming and require highly specialized skills, which presently makes its use as a standard industrial ergonomic assessment tool infeasible. Finally, the study only analyzed the acute loading of the joints and did not assess cumulative loading. This is problematic for the analysis of MMH in the supermarket sector, as most of the work tasks are highly repetitive. Hence, while many of the handling tasks do not appear hazardous in regards to the acute loading of the involved joints, they may well be hazardous if handled frequently over the course of a working day, week or year.
Based on the identified limitations, there are still challenges to be overcome to improve the versatility and applicability of this methodology for risk assessment of MMH. Most importantly, the methodology required imposing handling techniques on the workers to some extent as well as isolating each work task in a standardized manner to facilitate a somewhat efficient modelling of the tasks. Furthermore, as many of the trials had to be excluded for reasons outlined above, the methodology in its present state of development is best used for in-depth biomechanical analysis of heavier, two-handed lifts, which could possibly be identified as potentially hazardous using observational methods.
The study also contains several strengths. It was the first to apply state-of-the-art musculoskeletal models to estimate the dynamic loading of the knees, shoulders and L5-S1 joint during MMH using IMC-based kinematic data obtained in the field. The methods employed address well-known limitations of traditional quantitative lifting analysis tools, e.g. The Revised NIOSH Lifting Equation, which are typically limited to static or quasi-static analysis of sagittal plane lifting with both feet on the ground. Hence, many of these tools do not sufficiently replicate the real-life conditions of performing work-related MMH, which is a dynamic movement, involving both lifting, lowering, carrying and walking with the handled materials, while often being highly asymmetric. Despite the current limitations, the framework for estimating dynamic joint loads presented here shows great potential for risk assessment of MMH, and could help enhance our understanding of the dynamic loading of the involved joints in real-life working 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 Modeling System, which was used for the simulations. Mark de Zee is also a minority shareholder in the company.