Experimental procedure

Ten subjects (6 males, 4 females, age: 22 ± 1 years old, height: 176 ± 8 cm, mass: 72 ± 7 kg) participated in the experiment. The study was approved by the local institution’s Research Ethics Committee and each subject signed an informed consent form prior to the experiment. Subjects were young university students with very little manual material handling experience. They had had no injuries in the last six months. Five two-person teams were created, including mixed, women’s and men’s teams, with various height differences between members of each team (Table 1).

Table 1

Individual data on each team.

	Team member ♯1			Team member ♯2
	Height [cm]	Mass [kg]	Sex	Height [cm]	Mass [kg]	Sex
Team 1	165	59	F	169	63	F
Team 2	173	69	M	171	68	F
Team 3	178	69	M	191	86	M
Team 4	172	65	F	179	83	M
Team 5	186	74	M	179	82	M

Each trial was composed of two handling tasks (Fig 1): transfer a box from a lift location centered in relation to their body (facing the load to be lifted) to a deposit location 36 cm off-centered to the right of one team member and to the left of the other, and then bring it back again. Both lift and deposit locations were 4 cm above the ground and were marked with tape. Each team member began the trial in a neutral upright standing position with arms on each side of the body. After the experimenter gave the start signal, team member ♯1 was instructed to signal the other member of the tandem to begin the lift. After the first deposit, a short pause of approximately 2 seconds was taken in the neutral upright position without any load in their hands. Then team member ♯1 gave the signal to start the second box transfer (return). After the task, participants were asked to stand in a neutral position. Each handling trial thus comprised 9 successive phases: standing in a neutral position without a load in the hands; pre-maneuver 1 (starting at the first touch on the box and ending at the instant of lift); transfer 1; post-maneuver 1 (beginning at the instant of deposit and ending when the participants were no longer touching the box); pause; pre-maneuver 2; transfer 2; post-maneuver 2; and standing in a neutral position without a load in the hands.

Fig 1

Experimental protocol with two team members lifting the box at the beginning of transfer 1.

One box was adapted to have four different configurations: B17_E: 17 kg evenly distributed; B51_E: 51 kg evenly distributed; B51_♯1: 51 kg where the center of mass had been shifted toward team member ♯1; and B51_♯2: 51 kg where the center of mass had been shifted toward team member ♯2. Several weights were added to a rigid base structure to obtain the B51_E, B51_♯1 and B51_♯2 configurations. For each one, the same weights were added to ensure the same total mass. The box’s dimensions were 80.5 × 44.5 × 32 cm. Four repetitions of each tasks (including an outward and a return phase) were performed back to back with each configuration. Subjects were given a two-minute familiarization period to practice box transfers with configuration B51_E. After this familiarization period, trials with the B17_E configuration were performed. Then, the order of B51_E, B51_♯1 and B51_♯2 was randomized. Weight changes between conditions were executed out of sight of the participants.

A motion capture system (10 cameras, Vicon, Peak, UK) recorded at 100 Hz the 3D coordinates of markers located on both participants (see S1 File) and on the box. At the beginning of each experiment, the locations of 82 markers on each participant were recorded (Fig 2). Some of them were located on anatomical landmarks following the ISB recommendations [21, 22] and the others constituted clusters. Fifty-eight markers were left on each participant during handling trials. Markers on the participants were used to estimate anatomical landmarks position during tasks, then used as input to assess joint coordinates (more details about this step are provided under “Inverse kinematics”). In addition, 10 markers were located on the box. Four of them were used to construct a coordinate system associated with the box and the others were used to reconstruct these 4 markers in case of occlusion.

Fig 2

Experimental markers initially located on each participant.

GRF were measured at 2000 Hz using three force platforms (Advanced Mechanical Technology Inc., Watertown, MA). One subject stood on one of the force platforms and the second subject stood on two, with one foot on each platform. Subjects were instructed not to move their feet during the box transfers.

Biomechanical model

Back loading was estimated based on a whole-body osteoarticular model, composed of 16 rigid segments (pelvis, lower trunk, upper trunk, head, upper arms, lower arms, hands, upper legs, lower legs and feet) linked by 15 joints corresponding to 35 degrees of freedom (3 for the pelvis / lower trunk joint, 3 for the lower trunk / upper trunk joint, 3 for the neck, 3 for each shoulder, 2 for each elbow, 2 for each wrist, 3 for each hip, 1 for each knee and 2 for each ankle). The geometrical parameters were subject-specific, calibrated using motion capture data and an optimization-based method [23–25]. Body segment inertial parameters (BSIP) were extracted from anthropometric tables [26].

For the prediction methods, discrete contact points were defined on the model, corresponding to the potential contact points with the ground or with the box. Fourteen points under each foot and 11 points under each hand were defined to map the contact area [15, 16].

Inverse kinematics

From the positions of 36 anatomical landmarks (estimated from the positions of the rigid clusters), the joint coordinates were computed with a multibody optimization [27] and then filtered with a 4th-order Butterworth low-pass filter with a cut-off frequency of 6 Hz and no phase shift [28].

External forces

GRF were estimated by three different prediction methods (M1, M2 and M3), each compared to force platform measurement, considered as the gold standard (Fig 3). The three prediction methods differed in the type of experimental input data used: M1 used only motion data, M2 used motion data and the traction/compression force applied to the box and M3 used motion data and one team member’s GRF.

Fig 3

Setup of the gold standard (GS) and the three prediction methods (M1, M2, and M3) used to estimate GRF.

As the traction/compression force was not measured during these experiments, this force was calculated from the HF estimates. HF applied by each team member were estimated from the GRF and motion measured, by subtracting each body segment’s contribution to GRF [29]. The traction/compression force represents the minimum common HF applied by both team members in the longitudinal direction of the box (Eq 1) (Fig 4). The longitudinal direction was estimated using the box markers. At each sample time, it corresponded to a scalar value.

Fig 4

HF applied by both team members and their respective projections on the longitudinal axis of the box $$ \overrightarrow{y_B}$$.

Considering that the amplitude of the traction/compression force applied by each team member was equal and opposite, no resultant acceleration of the box emerged due to these forces. The acceleration of the box was directly proportional to the HF applied on the box by both team members along the three orthogonal directions. This point is discussed further under “Limits and perspectives”.

$$ \overrightarrow{H{F}_1}$$ and $$ \overrightarrow{H{F}_2}$$ are the HF applied by team member ♯1 and team member ♯2, respectively. $$ \overrightarrow{y_B}$$ represents the longitudinal axis of the box.

The rationale behind the addition of the traction/compression force applied to the box (Method M2) is that this is a component whose effect on the box cannot be captured from the whole-body and box kinematics. Thus, it is important information to ensure a good estimate of back loading.

The next three subsections describe the gold standard and the three prediction methods.

L5/S1 joint moments

L5/S1 joint moments were estimated with a recursive Newton-Euler algorithm [30], using a bottom-up approach. Measured or predicted GRF were used depending the method in question.

To estimate back loading during team handling, marker-based motion capture and GRF measurements are the most widely used methods in previous studies [31–33] and will be referred to as the gold standard in this study.

Analysis software

All prediction methods and L5/S1 joint moment estimation methods were implemented and processed with CusToM (Customizable Toolbox for Musculoskeletal simulation [34]). CusToM is a toolbox developed in Matlab^® that enables musculoskeletal analysis based on inverse dynamics approaches with a high level of customization.

Data analysis

Prediction methods were evaluated by comparing the predicted and measured GRF and the L5/S1 joint moments computed with the predicted and measured GRF during the transfer phases. The asymmetrical moment was defined as the vector sum of the lateral bending and torsion moment components. For each comparison, a Pearson correlation coefficient, a root mean square error (RMSE) and a relative RMSE (rRMSE) (Eq 2) [35] were computed. For the L5/S1 joint moments, the mean error was also computed. A positive value corresponded to an underestimate and a negative value corresponded to an overestimate.

A two-way repeated measures ANOVA was conducted on the Pearson correlation coefficient, the RMSE and the rRMSE by combining the three prediction methods (M1, M2 and M3) and the four box configurations (B17_E, B51_E, B51_♯1 and B51_♯2).

Method M3 can be used in two different ways: either using team member ♯1’s GRF or using team member ♯2’s GRF. Analyzing team member ♯1 using team member ♯1’s GRF is equivalent to using the gold standard method. Analyzing team member ♯2 using team member ♯2’s GRF is also equivalent to using the gold standard method. Thus, method M3 was evaluated by analyzing team member ♯1 using team member ♯2’s GRF and team member ♯2 using team member ♯1’s GRF. Comparison measures of each repetitive trial were averaged. Significant effects were analyzed using a Bonferroni post hoc test. The sphericity assumption was tested with Mauchly’s test and, when this assumption was violated, the Huynh–Feldt correction was applied. The significance level was set a priori at 0.05. All statistical analyses were performed with IBM SPSS^® software (version 26.0).

GRF estimation

For the vertical axis, no significant differences in any of the dependent variables were observed between Methods or Box configurations. The correlations were 0.99 and the errors approximately 2%. The vertical axis contained the greatest forces and amplitudes since it contained the weights of the subjects and the box.

Significant differences appeared for the antero-posterior axis. M1 had a lower correlation (mean: 0.31) and higher errors (mean: 54%) than the other methods. In the antero-posterior axis, the predicted forces were close to zero with M1, whereas each team member applied a non-zero opposite force (Fig 5). This force was considered in M2 and M3, resulting in a mean correlation of 0.90 with a mean error of less than 10%. There were fewer errors for B17_E than the other three box configurations (mean: 0.09 N/kg vs. approximately 0.25 N/kg for the other axes). Moreover, a significant interaction (p ≤ 0.001) was found between both factors for the antero-posterior axis (Fig 6). With M1, a lower error level was obtained with B17_E than the other box types, while no significant difference between box configurations was observed for the other two methods.

Fig 5

Representative example of predicted and measured GRF.

The three graphs on the left represent team member ♯1’s GRF and the three graphs on the right represent team member ♯2’s GRF. The orange curves represent the gold standard, the yellow curves M1, the green curves M2 and the blue curves M3. The different areas of the trial are highlighted: the handling phases are the ones with no background shading (areas on the left and on the right of each graph); these areas were taken into account in the statistics. The central area represents the break between two handling tasks. The phases shaded in gray are manipulation phases (pre-grip and post-deposit); these phases were not taken into account by the prediction methods.

Fig 6

RMSE for sagittal and asymmetrical L5/S1 joint moments according to prediction method and box configuration.

* significance level ≤ 0.05; ** significance level ≤ 0.01; *** significance level ≤ 0.001.

For the medio-lateral axis, few significant differences were observed between Methods or Box configurations (Tables 2 and 3); the mean correlation was about 0.84 with mean errors less than or equal to 10%. The medio-lateral axis was the axis with the lowest force amplitudes.

L5/S1 joint moment estimation

Significant differences in RMSE and rRMSE were observed for the sagittal axis. RMSE were the highest using M1 and the lowest with M2, where the mean error was 11 Nm. rRMSE were lower using M2 than with the other two methods. In addition, the RMSE for B17_E was lower than for the other three box configurations. Moreover, as with GRF in the antero-posterior axis, a significant interaction between Method and Box configuration (p ≤ 0.001) was found for the sagittal axis. With M1, a lower error level was obtained with B17_E than with the other box types, while the other two methods produced no significant difference between box configurations.

For the asymmetrical component, few significant differences were observed between methods (Tables 2 and 3). Since RMSE values were lower and rRMSE values higher, amplitudes were lower for the asymmetrical component than for the sagittal axis.

For the sagittal moment, the mean error was -27 Nm, -3.7 Nm, and 0.82 Nm for methods M1, M2 and M3, respectively. For the asymmetrical moment, the mean error was -3.9 Nm, -2.0 Nm, and -10 Nm for methods M1, M2 and M3, respectively.

Errors concerning sagittal and asymmetrical back loadings for each team member are reported in Table 4 for the sagittal axis. RMSE ranged between 16 and 52 Nm for M1, 4 and 19 Nm for M2 and 11 and 27 Nm for M3.

Table 4

RMSE for sagittal (S) and asymmetrical (A) L5/S1 joint moments for each team member and each prediction method.

RMSE [Nm]		Team member ♯1			Team member ♯2			Team (mean)
		M1	M2	M3	M1	M2	M3	M1	M2	M3
Team 1	S	29	11	11	39	19	11	34	15	11
	A	8.4	8.1	19	12	13	14	10	11	17
Team 2	S	16	10	19	52	16	27	34	13	23
	A	7.1	6.7	13	11	9.3	11	9.0	8.0	12
Team 3	S	22	4.2	15	39	12	25	30	8.2	20
	A	7.3	6.1	14	12	12	14	9.9	9.1	14
Team 4	S	29	11	20	26	11	13	27	11	16
	A	10	9.5	20	15	13	15	12	11	17
Team 5	S	19	9.5	14	16	10	11	18	10	12
	A	8.1	7.0	10	8.4	7.6	11	8.2	7.3	11

GRF estimation

The greatest difference between the prediction methods appeared on the antero-posterior axis, where M1 produced a much larger error than the other methods. This method, which considers only motion data, predicts only forces required to generate the box movement and not the traction/compression force applied by team members. This force can be applied to increase stability; it may reflect a lack of synchrony between the two handlers or may be applied voluntarily to decrease low back moments [20]. Because this force does not directly influence the motion of the box, it is necessary to acquire additional information to improve the accuracy of the estimate. The traction/compression force was used directly as an input in M2 and indirectly via the GRF of one of the team members for M3 and thus greatly reduced prediction errors in the antero-posterior direction. A recent study reported a similar problem when estimating medio-lateral force underneath each foot based solely on motion data [18].

For similar handling tasks performed individually, Larsen et al. [19] reported RMSE of 0.44 N/kg, 0.55 N/kg and 0.21 N/kg and Muller et al. [18] reported RMSE of 0.24 N/kg, 0.08 N/kg and 0.40 N/kg for the vertical, antero-posterior and medio-lateral GRF, respectively. These reported errors were the mean errors for each foot and therefore included the errors made in the distribution of forces between both feet. For this reason, these values are not directly comparable to those reported in this paper. However, when the traction/compression force is considered (M2 and M3), the errors are of the same order of magnitude as those reported in the literature for individual handling. The distribution of forces between team members is therefore well accommodated by the prediction methods, regardless of the box’s mass and the mass distribution.

L5/S1 joint moment estimation

The largest errors obtained with M1 for the sagittal L5/S1 joint moments came directly from the unconsidered antero-posterior forces. This error is particularly large for box configurations B51_E, B51_♯1 and B51_♯2, corresponding to the heaviest masses. The team members applied a higher traction force for these boxes, probably to ensure a certain stability. The greater the traction force applied, the greater the errors obtained with M1 for the sagittal L5/S1 joint moments. This phenomenon was not replicated with M2 and M3. Moreover, M3 seemed less accurate at predicting the sagittal L5/S1 joint moments than M2, even though the input data used were equivalent to the gold standard for one of the two team members. This is because, since force platform data are used for one team member, errors introduced by both team members’ kinematics and BSIP estimation interfere with GRF prediction for the other team member.

The mean errors reported for method M1 meant that, almost throughout the transfer, the moment was overestimated using the prediction method. This result is in accordance with those obtained by Barrett and Dennis [13], who indicated that the traction/compression force applied results in decreased back moments. For methods M2 and M3, no overall tendency to underestimate or overestimate was noted, but there was probably a random error around the reference value.

For asymmetrical moments, the error level was lower than for sagittal moments. However, since the amplitudes were lower, relative errors were slightly higher. Moreover, since antero-posterior forces had little influence on asymmetrical moments (only on the axial component of the moment), errors for M1 were smaller. Similarly, M3 had larger errors due to the accumulation of errors for both team members.

For similar handling tasks performed individually, Muller et al. [18] reported RMSE of 18.7 Nm, 10.6 Nm, and 12.6 Nm for the sagittal, frontal and transverse axes, respectively. These values are like those obtained with M2 and M3. Moreover, by comparing bottom-up and top-down approaches, several studies reported levels of uncertainty when estimating L5/S1 joint moments by using intrumented handles with force sensors and force platforms [32, 36, 37]. These uncertainties ranged between 4 Nm and 24 Nm for the sagittal axis, between 4 Nm and 20 Nm for the frontal axis and between 2 Nm and 15 Nm for the transverse axis. Using M2 or M3, the errors are in the range of uncertainty, regardless of box characteristics, team member or team composition.

General discussion

Team handling is commonly recommended to reduce the load during manual handling of heavy loads. Despite its benefits, some factors represent a risk for low back injuries and new methods should be developed to assess load among team members.

Based on the results of this paper, the use of only motion data (M1) did not allow an accurate estimate of L5/S1 joint moments. The lack of information on the traction/compression force applied by each team member could generate a significant error in the anterior direction and therefore in the L5/S1 joint moment estimation.

Adding the traction/compression force as an input (M2) is an improvement since it allows L5/S1 joint moments to be estimated with a relatively small error, comparable to the uncertainty values for this type of estimation. This method provides very encouraging results if the aim is to analyze realistic team handling tasks in a laboratory. Traction/compression force can be measured, for example, with strain gauges integrated into the box and positioned on its longitudinal axis. This instrumentation does not require the use of handles and does not limit handlers’ movements to an area defined by force platforms. Applying strain gauge measuring devices in the field could make it possible to study back efforts in a more natural handling situation although it would be technically complicated to do, given the work context and the loads to be moved. One viable method would be to estimate the traction/compression force and then predict GRF using method M2. Future work will be carried out to investigate this possibility.

Considering the traction/compression force in the form of force platforms under one of the team members is another interesting approach in a laboratory setting (M3). This makes it possible to estimate L5/S1 joint moments with a gold standard method for one team member and with a relatively small error for the other. When few force platforms are available, it is possible to assemble them to obtain a large displacement area for one team member. The other is then not constrained by any platform. Moreover, method M3 does not require any instrumentation of the displaced load, which can be complex, for example when instrumenting a stretcher loaded with a patient.

The proposed methods can be used to analyze the load distribution between handlers in various work situations. The experimental conditions used in this paper show that it is possible to work with handlers of different genders, sizes and weights, and to analyze the influence of the carried load’s characteristics by considering, for example, off-center loads.

Limits and perspectives

This work has several limitations. First, the tasks performed did not include foot displacements. This was constrained by the need to use force platforms. Other complementary studies involving foot displacements will have to be done, either by using several large platforms or by adding instrumented handles on the box. Nevertheless, this type of task has been used in the literature to validate prediction methods for handling tasks [18, 19]. Secondly, the traction/compression force was estimated from GRF and motion. These data were then considered as inputs for prediction method M2. It was reported that the estimation of hand forces using this type of method induces an error of about 20 N [29]. The influence of this error on the prediction method could be the object of a complementary study. Thirdly, the use of these prediction methods requires knowledge of different characteristic instants of a handling task, such as lifting and deposit, to separate the handling phases. Thus, the experimenters knew when a force was or was not applied to the handlers’ hands. Identification using a video camera might be the easiest way of doing this but requires time-consuming post-processing. To automate processing, identification methods based on handlers’ motion could be developed and implemented [18]. Fourthly, BSIP were extracted from anthropometric tables. Errors in BSIP estimation induce errors in GRF prediction [38] and in joint torque estimation [39, 40]. Improving the estimation of BSIP using calibration methods could improve the joint torque estimation, especially for subjects distant from the 50th percentile. Finally, to increase the realism of the team handling tasks to be analyzed, the ideal is to be able to perform experiments in the field in order to observe and analyze real work tasks. In this situation, the use of an optoelectronic motion capture system is no longer feasible. Previous studies have used the Lumbar Motion Monitor [2, 11] to measure low back kinematics in the field, but the proposed prediction methods require the whole body to be measured. An interesting extension of this current study would be to adapt prediction methods to wearable technologies such as magnetic and inertial measurement units [41–43], depth cameras [44, 45] or video cameras [46]. Since the accuracy of kinematic data plays a major role in prediction methods [38], studies must be carried out to verify the accuracy levels that can be obtained with this type of instrumentation.