Validation of Lumbar Compressive Force Simulation in Forward Flexion Condition

Abstract

Safety standard requirements must be implemented for lumbar support robots, which are mainly used for preventing low back pain (LBP) in caregivers. Usually, simulations are used to mimic actions that are not allowed for a real person. However, a comprehensive validation of a simulator in dynamic conditions has not been conducted. In this study, an ergonomic simulator is validated through forward flexion invasive experiments. The correspondence between the simulated and experimental compressive force (CF), as well as the CF obtained using two existing models about the unified angle, is investigated. The results show that the CF error between the measurements and the simulator at a flexion angle of 30∘is 11.8% and is lower than those obtained for the other two models (16.8% and 20.6%). Linear regression shows that the invasive data and estimated CF are close (slope = 1) in Merryweather’s model and CF simulator but not for Potvin’s model. We evaluate the precision of the simulator by using intraclass correlation coefficient method. Merryweather’s model is moderately consistent with invasive measurements, with R–0.685 and 0.627 at 0 and 30∘, while the CF simulator shows good consistency with Merryweather’s model with R–0.879 and 0836 at 0 and 30∘.

1. Introduction

Low back pain (LBP) is common health concern worldwide, and its prevalence increased by 18% in 2016 [1]. A major cause of LBP is the peak compressive force (CF) [2]. LBP is common among caregivers because they often exert high CFs [3]. Lumbar-type physical assistant robots are expected to reduce or mitigate the risk of LBP when carrying a heavy load [4].

In recent years, various lumbar-type assistant robots have been developed. By mimicking the back muscles using wearable soft bands, Abdoli developed a passive lifting assistive device to reduce the back muscular force [5]. Sankai and Kawamoto developed a lumbar-type assistant robot controlled by electromyography (EMG) signals [6]. Toxiri and Koopman developed an active robot using both EMG and kinematics measurements [7]. However, the safety of these robots has not been rigorously assessed.

Though safety requirements in the design stage of physical assistant robots are stated in the associated international safety standard ISO 13,482, the requirements are limited only to conceptual design guidelines, such as “A personal care robot shall be designed to minimize or reduce physical stress or strain to its user due to continuous use” [8]. In other words, the standard does not specify any quantitative criteria on the risk of LBP, although it identifies the hazards associated with physical assistant robots, which can be found in the use cases exemplified in Annex A of the standard. In this study, we evaluate whether the burden exerted at the lumbar vertebrae is safe considering the CF, the safety limit of which was originally determined to be 3.4 kN by the National Institute for Occupational Safety and Health (NIOSH) [9].

Previous studies have considered muscular damage due to LBP. van Diee¨ n reported that the trunk muscle recruitment pattern of patients with LBP is different from that of patients without LBP owing to the difference in spine stability [10]. Moreover, a postural control strategy can be adopted by people with LBP to maintain stability under trunk muscle fatigue [11]. Chan found that back pain occurs as trunk muscles exceed the limit to protect the tissues from further damage [12]. These studies imply that the spine stability should be considered to estimate the CF.

In this study, we apply only the 3.4 kN CF criterion, beyond which damage to lumbar disks is expected to cause LBP. It is worth noting that the CF exerted at a lumbar vertebra changes depending on the posture and motions even though they require a positive amount of maximum assistive force. This force is linked with a lumbar burden reduction ratio defined in JIS B 8456-1 which was determined using the machine/dummy method [13,14]. However, the machine/dummy method cannot represent the soft tissues nor mimic different motion patterns; these characteristics can be represented in human simulators. Simulators are expected to provide useful data to determine whether a caregiving posture that is captured and modeled for the later purpose of biomechanical analyses is sufficiently safe.

To avoid using any motion capture system, Potvin et al. proposed a static model to estimate the CF that exploits mechanical stability [15]. The focus in this model is the muscular and intersegmental forces concerned with the external load and human upper body weight. Merryweather et al. developed a quasi-dynamic model in which the muscular and intersegmental forces induce a static CF, and their coefficients are determined by the regression analysis of human flexion motion in a 3D simulator called 3DSSPP [16]. Finally, the quasi-dynamic CF is obtained by multiplying the static CF by a coefficient related to the flexion speed.

The existing estimation methods are limited to static/quasi-dynamic calculations, and their applicability to dynamic scenes has not been statistically verified. It has been reported that the static estimation is less than the dynamic estimation by 11–38%; thus, the static model cannot be applied in dynamic conditions [17].

Therefore, it is important to validate CF simulators with inverse dynamics functionality to estimate the CF. Currently, CF simulators can be easily constructed and combined with the moment arm proposed by Chaffin et al. [18]. In these simulators, the human movement is first recorded by a motion capture system. Second, the lumbar intersegmental force and moment are estimated using inverse dynamics [19,20]. Third, the CF is determined as a component of the resultant lumbar intersegmental force and the erector spinae force, which is perpendicular to the sacrum. By modeling the lumbar moment with a single muscle, the erector spinae is assumed to have a constant moment arm in the motion [21,22].

Although different CF estimation methods have been compared in previous studies, their dynamic performance was not comprehensively evaluated and invasive CF data were not considered [17,23]. Without the dynamic performance of the model and the invasive CF data, it is difficult to assess the accuracy of the models and their applicability in lumbar support robots.

In this study, we validate a CF simulator based on inverse dynamics computation using an ergonomic simulator. We also compare the data estimated by the CF simulator and those obtained by the regression models proposed by Potvin et al. and Merryweather et al. with the results of previous invasive CF experiments [15,16].

2. Method 2.1. Subjects and Tasks The Institutional Review Board at Nagoya University (Department of Engineering, No. 17-4) approved this project. Sixteen healthy male subjects (mean age: 25 y, standard deviation: 2.0 y; mean height: 1.71 m, standard deviation: 0.05 m; mean mass: 64.3 kg, standard deviation: 5 kg) with no history of LBP participated in the study.

Forward flexion is a practical method for estimating the load on the lumbar region [16,24]. Under conditions of no external load, subjects were asked to flex their backs forward ten times from an upright position to the highest possible posture in 2 and 10 s. The extension process took the same time.

The effect of training the subjects was not considered in the presented study. In contrast, to obtain the CF data from them in a severer situation, we considered that the basic flexion could result in worse postures than those of the trained caregiver. In the next stage, we will train subjects in terms of lifting, transfer, and positioning with a Hoyer lift. 2.2. Instrumentation

An optical motion-capture (mocap) system with three MAC3D cameras (Motion Analysis Co., Santa Rosa, CA, USA) was used to acquire the motion data as shown in Figure 1. More than 40 optical markers were taped directly to the subjects’ skin to represent the feature points of each body part. The markers were placed on the head, shoulders, hands, waists, arms, legs, feet, hips, and back, including the lumbar spine (from T12 to S1.). Two force plates (M3D-EL-FP-U type, Tec Gihan Co., Kyoto, Japan) were used to measure the ground reaction force at each foot, one at the ball of the foot and the other at the heel. In this paper, the flexion level is presented as the “flexion angle”, which is the relative angle between the axis through S1 and C7 in the standing posture and flexed posture, as shown in Figure 1.

The motion and analog force data were synchronized by the mocap system and filtered by a fourth-order low-pass filter at a cut-off frequency of 10 Hz. 2.3. CF Estimation

2.3.1. CF Simulator

The CF simulator estimates the CF in two steps: inverse dynamic computation and lumbar muscular force computation. First, a virtual human model was reconstructed in the ergonomic assessment system DhaibaWorks, which was developed by the Japanese National Institute of Advanced Industrial Science and Technology (AIST) [25]. In this model, body segments are assumed to be coupled links. The intersegmental forces and joint torque of the body segments are computed by the inverse dynamic algorithm [26]. Second, the net joint torque at L5/S1^N(L5/S1) exerted by all lumbar muscles, which can be regarded as being produced by a representative muscle [21] as shown in Figure 2a, was calculated as follows:

^N(L5/S1)=_Fm(L5/S1)×r,

where_Fm(L5/S1)represents the representative muscular force (N).r is the moment arm vector; its direction is from L5/S1 and perpendicular to the representative muscle line of action, and its length is set to 5.3 cm [18].

In forward flexion motion (2D case), the representative muscle is the erector spinae, and the muscular force exerted by the erector spinae muscle maintains the lumbar stability at L5/S1. The mechanical equilibrium at L5/S1 is shown in Figure 2a. The magnitude of the CF,C_Fd(N), and shear force,S_Fd(N), can be expressed as follows:

C_Fd≡∥_CFd∥=∥^F(L5/S1)∥cos(α)+∥_Fm(L5/S1)∥,

S_Fd≡∥_SFd∥=∥^F(L5/S1)∥sin(α),

where^F(L5/S1)(N) is the intersegmental force at joint L5/S1, andα(rad) is the angle between^F(L5/S1)and the axispS1perpendicular to the sacrum surface through S1.

2.3.2. Model of Merryweather

As shown in Figure 2b, Merryweather established a model to estimate the CF. In this model, the magnitude of the CFC_Fmis expressed as follows:

C_Fm≡∥_CFm∥=k×[0.0167(BW)(H)sin(β)+0.145(L)(HB)+0.8(BW/2+L)+23]×g,

where k,BW, H,HB, L,β , and g are the speed-related coefficient (slow = 1.15, moderate = 1.3, fast = 1.4), body weight (kg), height (cm), horizontal distance from load to L5/S1(cm), load carried by the subjects (kg), angle between the trunk in different positions (C7-S1), vertical direction (rad), and gravitational constant, respectively [16].

Merryweather et al. assumed that CF consists of three components. The first is the muscular force of the upper body weight (BW/2) exerted at the upper center of mass (CoM) with respect to L5/S1; the second represents the muscular force of the external load with respect to L5/S1; the third term represents the force exerted by the upper body mass and the external load. The moment arm of the back muscle (erector spinae) is incorporated into the coefficients 0.0167 and 0.145. The coefficient of each component is determined by the regression analysis of the simulation result by a 3D ergonomic simulator 3DSSPP [16].

2.3.3. Model of Potvin

Potvin developed a static human model by directly using mechanical equilibrium and regression analyses. In this model, the magnitude of the CF,C_Fp, can be expressed as follows:

C_Fp≡∥_CFp∥=(0.54BW×g×WB+L×g×HB)/0.06+sin(γ)×(0.54BW+L)×g,

whereBWis the body weight (kg) (the upper body weight is 54% of theBWand acts at the upper CoM);WBrepresents the moment arm (m) of the upper body’s CoM; L is the external load (kg);HBis the moment arm (m) of the external load about L5/S1, andγis the angle between the trunk (C7-S1) and the horizontal direction (rad).

WBandHBare obtained by two linear regression models:

WB=−0.032+0.104LA−0.004V+0.000356γ+0.354cos(γ),

HB=−0.144+0.991LA−0.031V+0.000193γ+0.281cos(γ),

whereLA is the horizontal distance between L and the ankles, and V is the vertical distance between the ground and L. On the right side of Equation (5), the first component represents the amount of muscular force exerted by the upper body weight and the external load, the moment arm of the back muscle (erector spinae) was determined to be 0.06 m, and the second component represents the amount of the contact force between L5 and S1 [15].

The characteristics of the three models (i.e., the models by Potvin et al., Merryweather et al., and the CF simulator proposed in this study) are compared in Table 1. From Figure 1 and Table 1, it can be seen that the angles of the three models are different, and it is not convenient for us to find the relationship between the CF and flexion posture. Thus, only the angleβandγare used for calculating the CF. To display the relationship between the flexion level and the CF, we use the flexion angleθ.

The CFs at the L5/S1 level estimated by Merryweather et al. and Potvin et al. were compared with those estimated by the CF simulator [15,16].

There are other models that can be used to estimate CF, such as the statistical method, EMG method, and optimization method [27,28,29]. However, these models may yield some errors in motions. For example, in the statistical method, an error larger than 20% can occur as the torque is smaller than 10 N·m [27]. In the EMG-assisted method, a part of muscular force can be provided by passive tissues which cannot be monitored by EMG device (e.g., flexion-relaxation phenomena). Thus, muscular force and the compressive force can be underestimated. The accuracy of the optimization criteria on presenting muscular forces has not been validated. Hence, we select the models of Potvin and Merryweather because they are potentially applied to more realistic conditions than other models.

2.3.4. CF Estimated by Invasive Data

The reported invasive measurements were used to evaluate the performance of the CF simulator. The amount of CF estimated from the invasive experimental measurement isC_Fi, which is obtained from the intradiscal pressure (IDP ) measured by the invasive experiment and is used to evaluate the model’s precision [30,31,32,33].

C_Fican be obtained as follows:

C_Fi=IDP×_Adisc×c,

where_Adiscis the cross-sectional area of the intervertebral disc, and c is the correction factor to convertIDPto CF. CF is calculated from theIDPusing a conversion coefficient because the property of the intervertebral discs is not uniform. The conversion coefficient, or correction factor, is used to determine the meanIDP per unit area. The proposed range of this correction factor is 0.55–0.77 [30], and a mean value of 0.66 was used in recent studies [34,35]. Dreischarf et al. established a finite element model and pre-determined several values (0.77, 0.66, and 1) for the correction factor in the L5/S1 disc [36]. The CF had a strong correspondence with the pre-load simulated value at the L5/S1 level when the correction factor was 0.66. Therefore, in this study, a correction factor of 0.66 was adopted to convert theIDPto CF.

The estimated and experimental CFs can be compared with respect to the unified flexion angle (θ ). In Figure 3, the precision of the CF simulator is compared with data from previous invasive experiments [30,31,32,33]. Detailed information of these experiments is presented in Table 2.

2.3.5. Validation

The estimated CFs from each model for ten subjects and invasive data at two flexion angles (0^∘, 30^∘) were compared. We used the Kruskal–Wallis test/Mann–Whitney test with Boferroni procedure to analyze whether estimated and invasive data were significantly different from each other.

A scatter plot regression method was employed to represent the relationship between simulation and experimental results [37]. The CF determined from invasive data were considered to be the real CF. The CFs obtained from the models with ten subjects were compared with invasive result under the same flexion angle. The CF was estimated at flexion angles of 0^∘, 10^∘, 20^∘, and 30^∘ , which were selected from the literature [30,31,32,33]. Subsequently, we built a relationship between the invasive experimental and the modeled CFs, representing the horizontal and vertical axes, respectively.

The intraclass correlation coefficient (ICC) is a reliability evaluation method that is used for quantitatively evaluating the precision of the CF simulator in this study [38]. However, the body parameters and kinematics data reported from previous invasive experiments are not sufficient to reconstruct the motions in the CF simulator. To solve this problem, we adopted an indirect method: first, we attempted to find a statistical model that is at least moderately consistent with the invasive data by the ICC method; and second, we evaluated the precision of the CF simulator with respect to this statistical model by the ICC method again. We labeled the statistical model in this procedure as “precision evaluating model”. In the second use of the ICC method, sixteen participants were randomly divided into two sub-groups: the sub-group for the CF simulator (mean age: 25.5 y, standard deviation: 2.0 y; mean height: 1.72 m, standard deviation: 0.045 m; mean mass: 64.5 kg, standard deviation: 5.5 kg), and the sub-group for Merryweather’s model (mean age: 24.5 y, standard deviation: 2.1 y; mean height: 1.70 m, standard deviation: 0.05 m; mean mass: 64.2 kg, standard deviation: 5.5 kg).

3. Result 3.1. CF Estimation

Figure 4 shows the comparison between the estimated CF obtained with the proposed simulator, the model by Potvin [15], and the model by Merryweather [16], and the modified invasive measurements for forward flexion angles of 0^∘, 10^∘, 20^∘, and 30^∘ . The mean values of L5/S1 at the different forward flexion angles for all estimations (models and experiment) are shown in Figure 4. The CF estimated by the CF simulator at 0^∘, 10^∘, 20^∘, and 30^∘ is 416, 841, 1244, and 1546 N, respectively. Table 3 presents the difference in the CF between adjacent flexion angle ranges of 0^∘–10^∘, 10^∘–20^∘, and 20^∘–30^∘, hereinafter referred to as AR1, AR2, and AR3, respectively. These differences are different among the different models: the increase in CF obtained by the CF simulator between AR1 and AR2 and between AR2 and AR3 is −2.2 and −10.1 N/deg, respectively.

3.2. Influence of Dynamic Motion

Figure 5 shows the comparison between two- and ten-second flexions in terms of the change in the mean angle, angular velocity, and mean angular acceleration in all the trials from one subject. In Figure 5, the whole motion, including flexion and recovery within 2 and 10 s, is normalized between 0 and 100%. The largest difference in the angular velocity was observed at point A, and the corresponding flexion angle was approximately 48^∘. The largest difference in angular acceleration was observed at point B for a flexion angle (vector S1-C7 relative angle before and after being rotated) of 90^∘ . Figure 6 shows the estimation with and without the 20 kg load at different flexion speeds. The CFs at 48^∘and 90^∘for the two- and ten-second flexion motions were similar.

3.3. Regression Analysis

Figure 7 shows the linear regression fit of the models obtained by comparing the CF estimated using the models with that obtained from invasive measurements under the same flexion angle. It can be seen that the intercept becomes 0, and only the slope reveals the relation of data close to the vertical axis (estimated CF from models) and the horizontal axis data (CF obtained by invasive measurements). As the slope between the CF and flexion angle (Figure 7a–c) approaches 1, the correlation between the invasive measurements and the data estimated using the corresponding model increases. As shown in Table 4, the slopes obtained using the dynamic CF simulator, Merryweather’s model, and Potvin’s model are 1.012, 1.070, and 0.641 with^R2of 0.683, 0.658, and 0.675, respectively; their root mean square errors (RMSEs) are 269, 259, and 205 respectively. Thus, we selected Merryweather’s model as the precision evaluating model.

3.4. Precision Analysis

In the first step of precision evaluation, the reliability of Merryweather’s model when comparing with invasive measurement using the ICC method was 0.685 and 0.627 at 0 and 30^∘, respectively. Therefore, the consistency of Merryweather’s model with invasive measurement is moderate (0.5–0.75), which can be accepted as the precision evaluating model. In the second step, the R between Merryweather’s model and the CF simulator is 0.879 and 0.836 at 0 and 30^∘ , respectively. This shows a good consistency (0.75–0.9) between the two models [38].

4. Discussion

The accuracy of simulators that provide validation data of the safety standards of physical assistant robots must be high. We compared the performance of three types of biomechanical models with the data obtained from an invasive experiment involving ten healthy subjects. We hypothesized that a CF simulator can improve the precision of CF estimation over the models proposed by Merryweather and Potvin under dynamic conditions [15,16]. Under forward flexion conditions, the results of each model were compared with those of invasive measurements [30,31,32,33].

The CF was estimated under the flexion angleθof 0^∘, 10^∘, 20^∘, and 30^∘. A Wilcoxon test was performed to evaluate the difference between the CFs obtained with the three models and the invasive measurements at L5/S1 and 0^∘and 30^∘. The angles of 10^∘and 20^∘ were excluded because of the lack of invasive data at those angles. Among all body parameters, the body height and weight have the greatest influence on the CF estimation [9]. The results of the Mann–Whitney test imply that there is no significant difference in the height and weight between the subjects in the present study and those in the invasive experiments. Thus, we consider the two groups of subjects to be in the same population and the group difference in height and weight to be negligible small. As a result, the CF of each model significantly differed from the CF obtained from the invasive measurements. Furthermore, as shown in Figure 4, the CF discrepancies imply that the model by Merryweather and the proposed CF simulator yield results that are more similar to those of the invasive measurements than those of the model by Potvin. The ratio of the increment of CF increasing with respect to flexion angle observed in ratio of the invasive experiment decreases from 49 to 11.5 N/deg asθ increases from AR1 to AR3. Only the CF simulator and the model by Potvin et al. show a trend similar to that of the invasive measurements, while the ratio of the increment of compressive force CF with respect to flexion angle obtained by Merryweather’s model increases from 35 to 37 N/deg from AR1 to AR2 and decreases from AR2 to AR3. Therefore, the CF simulator has a good correspondence with the invasive experiment data, in terms of both the trend and the absolute value. In Figure 4, the error of the CF simulator, compared with the experimental results, at 30^∘is 11.8%, which is much lower than the 16.8% and 20.6% errors of the models by Merryweather and Potvin, respectively.

The differences between the CF simulator and the model by Potvin is probably due to the dynamic/static conditions: at 0^∘, 10^∘, 20^∘, and 30^∘, the CF estimated with Potvin’s model is 241 N (52%), 326 N (38%), 426 N (34%), and 450 N (29%) smaller than that estimated with the CF simulator, respectively. This result is in good agreement with the results of van Diee¨ n, who reported an error range of 11–38% of dynamic estimation [17].

We can conclude that Merryweather’s model performs better than Potvin’s one when comparing with invasive data. The possible reasons are as follows. First, Potvin’s model is originally based on static condition where the contribution of dynamic motion is not taken into account, while Merryweather’s model corrects the estimated static CF by using a dynamic coefficient k > 1. In the stoop case, which is the worst condition in forward flexion, k approaches 1.15. Second, the model of Potvin is based on the mechanical equilibrium at lumbar, while Merryweather modified the constant term of the mechanical equation for compensating the error of the estimated result with respect to the actual human body condition. Third, in Potvin’s model, muscular force could be smaller than zero. This situation happens because the negative value of WB and HB can be obtained in Potvin’s model under the condition that LA approaches zero or when the load is horizontally close to the ankle. However, muscular force should always be positive. The muscular force with a negative sign implies that the CF estimated in Potvin’s model is smaller than the actual value; however, Merryweather’s equation does not contain the negative moment arm; thus, the muscular force is positive.

IDPs are usually measured at the L4/L5 level rather than the L5/S1 level. However, the difference between theIDPs measured at these positions is small. The in vitro experiment showed that theIDP obtained at L4/L5 was in good agreement with that obtained at L5/S1. As shown in Table 5, the errors measured from the upright posture and flexion motion are as small as 0.02 and 0.03 MPa, respectively [39,40,41]. Therefore, the difference in the CF acting on L4/L5 and L5/S1 is considered negligible, and the invasive measurements measured at the L4/L5 level are used for comparison. The difference of 0.02 and 0.03 MPa is small and corresponds to only 3% of the total L5/S1 stress in both upright and flexion postures. Therefore, we used the invasive data at the L4/L5 level. Furthermore, based on a Mann–Whitney test, the Cabello’s measurement with six specimens shows no significant difference (p> 0.05) between the stress at L4/L5 and L5/S1 at the significance level of 0.05 regardless of the posture. Thus, it is reasonable to use theIDPmeasured at the L4/L5 level instead of that at the L5/S1 level, which is difficult to measure in the invasive experiment.

Muscle co-activation also affects the forward flexion. We focused on the extent of the contribution of the co-activation on the CF. van Diee¨ n performed an EMG-based co-contraction analysis using lumbar muscles’ EMG [42]. He measured the preliminary, maximum isometric contractions of those muscles and estimated the muscular force of each muscle with the EMG-based model. Their result suggested a linear relationship between the peak net moment and peak force at L5/S1 during both 2D flexion and 3D flexion movements. This implies that the single back muscle model can also be used for both 2D and 3D models only considering CF [42]. Thus, we consider that the single muscle model is suitable for forward flexion motion.

In the single-back muscle model, the representative muscular moment arm is important for muscular force estimation. The erector spinae is considered the major muscle in forward flexion [43]. The muscular moment arm is the distance between the vertebral body center and the erector spinae center in the sagittal plane. To obtain the muscular moment arm, the distance between the erector spinae geometrical center and the neutral line, which passes through the vertebral body center in the transverse plane, is calculated by magnetic resonance imaging (MRI) or computed tomography (CT) [44]. Table 6 shows the muscular moment arm estimated in previous studies. Most estimations are in the 5.3–7 cm range [18,44,45,46,47]. In this study, the lower limit (5.3 cm) was used to estimate the CF and evaluate extreme lumbar conditions. This value is also in the range of representative moment arm length estimated by van Diee¨ n, which is suggested to be between 5–6 cm [42].

The geometry of trunk muscles including erector spinae (ES), quadratus lumborum (QL), pasos (PA), internal oblique (IO), external oblique (EO), rectus abdominis (RA), and latissimus dorsi (LA) has been reported by McGill, Cholewicki, and Stokes [48,49,50]. McGill estimated the relative displacement between L4/L5 level and the muscles’ endpoints [48]. Cholewicki presented the L4/L5 but each level from the ribcage to the sacrum and described a more detailed model for multifidus and erector spinae [49]. Stokes presented a more detailed model of EO and IO than Cholewicki [50]. However, the pattern of the muscle recruitment varies among people in lateral flexion and twisting. It is necessary to obtain how the muscles quantitatively affect in lateral flexion and twisting, as it has been found in forward flexion where erector spinae is the major muscle and the moment arm can be considered as a constant. Even though the representative moment arm is difficult to use aside from forward flexion, it is possible to establish the relationship between the CF and the dominant factors, such as flexion angle, twisting angle, and the external load based on the regression analysis. Therefore, the performance of this model in lateral flexion and twisting motion remains to be investigated in the future.

Our experiments were conducted in a 2D plane, which cannot represent realistic movement performed during caregiving tasks. We consider, however, that the forward flexion motions can represent actual caregiving transfer task motions. The critical issue is whether the CF estimated in 3D caregiving motions is similar to that obtained from a 2D model. Kingma reported that, in manual material transfer tasks, there is a minor L5/S1 moment difference between 2D and 3D tasks for the asymmetric handling task with the angle ranging from 0^∘to 90^∘ [20]. Considering Kingma’s data, if the erector spinae’s moment arm is taken as 5.3 cm, as proposed by Chaffin et al. (2006), the error in the CF at 10^∘, 30^∘, 60^∘, and 90^∘of CF at 0^∘ is 1.7%, 1.7%, 0%, and 1.3%, respectively [20]. These errors imply that transfer motions can be simplified by 2D flexion when taking CF as the safety assessment standard. By studying 2D flexion, we can estimate the peak CF achieved during caregiving tasks.

In this study, we assume that CF is one of the major factors contributing to back pain and the CF should be similar at the same flexion angle in symmetric forward flexion and extension process. The second assumption is supported by Wilke’s invasive data [51]. However, extension could result in back pain not because of the lumbar but instead the hip: in the flexion-extension process, LBP patients are more likely to have an increase in gluteal fatigability related to LBP in extension than healthy people [52]. In our study, we aim to validate the CF simulator for standardization of lumbar safety. Although the hip muscle is related to LBP, this factor cannot be quantitatively evaluated for LBP currently. We emphasize that compressive force is suitable for standardization.

For each sample, the values along the horizontal axis in Figure 7a–c were obtained from different flexion postures in the invasive experiment, and the values in the vertical axis were obtained by each model under the same conditions. The linear relationship helps to understand whether the estimated data are in good agreement with the experimental data. In Figure 7d, the difference between the models is presented by the distance (D) between the marker and the line y = x shown in Figure 7a–c. It can be seen that the CFs estimated using the Kruskal–Wallis test on the D of the Merryweather, Potvin, and CF simulator are significantly different (p≤ 0.001). In Table 4, the slope implies that both the CF simulator and the model by Merryweather perform better than the model by Potvin. Thus, we conclude that dynamic models can obtain more accurate CF estimates at L5/S1 than static ones. The discrepancy between the model of Potvin and the invasive measurements is probably due to the difference in the CF during static and dynamic motion.

The Kruskal–Wallis test results show a significant difference (p≤ 0.01) among the estimates of the models and the invasive data. The Mann–Whitney test is performed between the estimated CF by each model and the invasive data, respectively. Considering the small amount of invasive data available, the significance level of invasive data and estimated CF by each model is set asα= 0.01. Potvin’s data show a highly significant difference with the invasive data in both 0^∘and 30^∘flexion conditions. A small difference is observed between the CF simulator and the Merryweather’s model. Even though the anthropometric data vary, the CF estimated with the simulator and Merryweather model have a smaller difference from the invasive measurement than the estimates of Potvin’s model. However, in the 3D condition, Merryweather’s model cannot estimate the external forces other than those exerted by the hands. During caregiving tasks performed with arms, which shorten the moment arm of external loads, Merryweather’s model can overestimate the moment at L5/S1. Under this condition, the CF simulator can estimate the external force acting on any point of the body. Because the back muscle model can also be used in 3D cases, the simulator can be applied in various caregiving conditions.

None of the previous studies quantitatively evaluated the precision of a model by using invasive data [17,23]. An important reason is that invasive data are reported in only a few conditions and the subjects are limited. Thus, we compared only 0 and 30^∘conditions between invasive data and estimated data, where most of the invasive data are reported. Although the consistency between Merryweather’s model and invasive measurement is moderate, it is promised expected to promote the reliability if more invasive data are reported in the future. Even so, not all simulators can reconstruct the motion and simulate the conditions while relying only on the reported data of invasive experimental conditions. A simpler “precision evaluating model” is taken as the bridge between invasive measurement and some complex model, such as the CF simulator; if the precision of the precision evaluating model is accepted, a complex model in good or excellent consistency with the precision evaluating model is also accepted.

5. Limitation

Validation was performed with a relatively small group of ten subjects. They did not represent a wide range of ages, body weights, and body heights; hence, these parameters are not investigated in this study. These healthy subjects did not include the elderly. In different age groups, the compressive strength at the lumbar disc of the elderly is less, approximately 1.7 kN, than that of young people under the 10% risk level [53]. Under the same CF, the elderly may have a higher risk of low back pain.

With the exception of lumbar criterion, arm strength is an important factor of flexion/lifting motion. The criterion of arm strength is not included in this study because the strengths of nonpreferred and preferred arms show a considerable difference. This difference may cause errors in standardization for the safety requirement. We did not simulate the restricted hip motion of caregivers with hip pathology, which is another limitation of this study. One reason is that the different severity of hip pathology may lead to different restricted range of hip motion; thus, it is difficult to select an appropriate restriction level. Another reason is that the restriction degree of one joint motion may affect other joint motions. Caregivers with hip pathology have to adopt inappropriate posture in caregiving tasks, thus experiencing more CF at the lumbar. Therefore, in this study, we asked the subject to maintain a stoop posture to present the worst case in forward flexion. Considering that the peak CF usually occurs at the flexion/lifting period in transfer tasks, we validated the CF simulator in the dynamic flexion condition. The twisting motion is not presented in this model because the erector spinae muscular force cannot provide most of the required twisting torque; hence, the evaluation of twisting motion requires a muscle model that can accurately present all the activated muscles at the lumbar region. Because the recruitment pattern of the lumbar muscles in twisting cannot be accurately obtained, the way in which twisting with external load in the musculoskeletal model is presented remains to be investigated.

The lumbar correction factor c range from 0.55 to 0.77 with the mean value of 0.66 [30,36]. Even though the correction factor varies for each individual, we cannot clearly determine the variations. Thus, we directly took the mean value 0.66 for all individuals as it was used by Dreischarf et al. and Nachemson et al. [30,36].

The individual differences affect the estimation results; the participants in this study were different from those in the invasive experiments. Because the variance of the invasive measurements affects the reliability of the data, it is necessary to investigate the factors causing the variance. The coefficient of variation of the erector spinae’s EMG is reported in resisted back extension with a wide range of 0.16–0.38 [54]. However, only a slight difference was observed in the body parameters of the subjects: ages of 18 ± 1.2 yr , heights of 176.5 ± 3.2 cm, and body masses of 71 ± 4.5 kg. This implies that the error on invasive measurements could not be reduced by only normalizing the body parameters or fixing the motion pattern. Even though this study includes almost all previous invasive data, the difference between the invasive data of our participants and those of previous studies is difficult to estimate and introduces inaccuracy in the quantitative evaluation of the three models.

6. Conclusions In this study, we validated a CF simulator by applying it to lumbar-type assistant robots with standard ISO 13,482 safety requirements. The validation was performed in forward flexion using previous invasive measurements for comparison. In addition to the CF simulator, the performance of two existing models, namely those proposed by Merryweather and Potvin, were evaluated and compared. Using the ICC method, Merryweather’s model was selected as the “precision evaluating model” to indirectly evaluate the precision of the CF simulator which was established based on the inverse dynamics method and single muscle model.

When considering the increment in CF with respect to the flexion angle, only the CF simulator was in agreement with the decreasing trend of the increment in CF observed in the invasive measurements (as shown in Table 3). But none of the models can perfectly show the same tendency with the invasive measurement. During the forward flexion, a significant difference between the estimated CF from Merryweather’s model and the CF simulator in the regression analysis was found; the model by Merryweather et al. and the proposed CF simulator showed similar performance and yielded results that were more similar to those of the invasive measurements than those of the model by Potvin et al. Potvin’s model is regarded as different from other two models. In the ICC precision analysis, Merryweather’s model showed a moderate consistency with invasive measurement at 0 and 30^∘; in the second step, the CF simulator had a good correspondence with Merryweather’s model at both angles. From these results, we conclude that the CF simulator yields results that are similar to those of invasive measurements.

Figure 1. Schematic of experimental setup for obtaining the ground reaction force and motions. Four force plates were used to measure the ground reaction force on both feet, and optical markers were placed all over the body and monitored by the motion-capture system. The external box with a mass of 20 kg was simulated in the compressive force (CF) simulator. The flexion angleθwas considered as an indicator of the flexion level.

Figure 2. Computation of the compressive force (CF) at L5/S1 using (a) the CF dynamic simulator (N(L5/S1)andF(L5/S1) are first obtained by inverse dynamic calculations), (b) Merryweather's model [16], and (c) Potvin's model [15]; the calculation is based on the global coordinate system oxy.

Figure 3. Flowchart of compressive force (CF) estimation procedure and validation of the CF simulator. Each motion was recorded by a motion-capture system, and kinematics data were used for estimating the CF in the CF simulator and in the models by Potvin and Merryweather [15,16]. The data were compared with invasive experimental data to evaluate the precision of each model.

Figure 4. Comparison between different models and invasive measurement under free flexion. The error bar signifies the standard deviation (STD) (N) (*p≤ 0.05; **p≤ 0.01).

Figure 5.(a) Flexion angle versus flexion process for 2 and 10 s movements; 'A' represents the greatest difference in angular velocity between 2 and 10 s, which is approximately 48∘, and 'B' represents the greatest difference in angular acceleration between 2 and 10 s, which is approximately 90∘; the flexion cycle is period that starts in an upright posture via flexion and ends in an upright posture again; (b) Trunk angular velocity versus flexion process for 2 and 10 s flexions; (c) trunk angular acceleration versus flexion process for 2 and 10 s flexions.

Figure 6. Compressive force (CF) estimated by CF simulator at 0∘, 48∘, and 90∘with and without a 20 kg external mass.

Figure 7. Regression analysis between the estimated compressive force (CF) obtained from (a) the model by Merryweather et al. (2009), (b) the CF simulator, and (c) the model by Potvin et al. (1997) with respect to the CF determined from invasive measurements; (d) Kruskal-Wallis test results: significance level (p≤ 0.001) among the distances between markers and y = x line estimated by the three models.

Model	CF Simulator	Merryweather [16]	Potvin [15]
Type	Inverse dynamic model	Regression model	Regression model
Conditions	Dynamic	Quasi-dynamic	Static
Torso angles base line	C7/S1 in standing/flexion	Vertical and C7/S1	C7/S1 and horizontal

Researcher	Information	Disc Level	Adisc(cm2)	Postures	IDP (MPa )	CF (N)
Nachemson [30]	Two subjects 49–52 y 67–75 kg 163–175 cm	L4/L5	19.7	Standing	0.60 *	780
Wilke [31]	One subject 45 y 70 kg 168 cm	L4/L5	18	Standing Forward flexion at 36∘ Standing (19.8 kg) lifting (19.8 kg)	0.5 1.08 1 2.3	594 1283 1188 2732
Sato [32]	Eight subjects 22–29 y 60–90 kg 166–181cm	L4/L5	15.9	Standing Forward flexion at 30∘	0.53 1.32	556 1385
Takahashi [33]	Three subjects 24–26 y 70–77 kg 170–180 cm	L4/L5	19.1	Standing Forward flexion at 30∘Standing (10 kg) Forward flexion (10 kg)	0.34 1.21 0.416 1.46	556 1385 428 1525

* Nachemson’s measurement is corrected by Dreischarf et al. [36].

N = 10	Increment of the CF at L5/S1 (N/deg)
Angle range (∘)	AR1 (from 0∘to 10∘)	AR2 (from 10∘to 20∘)	AR3 (from 20∘to 30∘ )
Invasive measurement	49.0	26.6	11.5
CF simulator	42.5	40.3	30.2
Merryweather [16]	35.0	37.2	31.5
Potvin [15]	34.0	30.3	27.9

Model	Slope	R2	RMSE (N)
Merryweather [16]	1.070	0.658	259
Potvin [15]	0.641	0.675	205
CF simulator	1.012	0.683	269

Researchers	Characteristic	Lumbar Disc Level	IDP(MPa )	Load (N)
Dolan [39]	22 subjects 19–96 years old 174 cm	L4/L5 L5/S1	0.74 (standing upright) 0.73 (standing upright)	1000
Cabello [40]	6 subjects 49.5 (33–66) years old	L4/L5 L5/S1	0.84 (standing upright) 0.78 (standing upright) 0.99 (20∘flexion) 0.97 (20∘flexion)	100–750

Researcher	Subjects	Length (cm)
McGill [45]	13	8.5(0.68)
Kumar [44]	8 males, 5 females	6.00(0.37), 5.85(0.49)
Nèmeth [46]	11 males, 10 females	7.1 ± 2, 6.5 ± 2
Daggfeldt [47]	4 males	5.3

Author Contributions

Conceptualization, Y.Y.; Formal analysis, X.X.; Investigation, X.X.; Methodology, X.X.; Project administration, Y.A.; Resources, Z.T.; Supervision, Y.Y.; Visualization, N.K.; Writing-original draft, X.X.; Writing-review & editing, Y.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Japan Agency for Medical Research and Development (AMED) ⟨Japan⟩ (⟨JP16he1202004⟩).

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the colleagues of the National Institute of Advanced Industrial Science and Technology (AIST), ⟨Japan⟩.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Vos, T.; Abajobir, A.A.; Abate, K.H. Global, regional, and national incidence, prevalence, and years lived with disability for 328 diseases and injuries for 195 countries, 1990-2016: A systematic analysis for the Global Burden of Disease Study 2016. Lancet 2017, 390, 1211-1259.

2. Marras, W.S.; Jorgensen, M.J.; Davis, K.G. Effect of foot movement and an elastic lumbar back support on spinal loading during free-dynamic symmetric and asymmetric lifting exertions. Ergonomics 2000, 43, 653-668.

3. Burdorf, A.; Koppelaar, E.; Evanoff, B. Assessment of the impact of lifting device use on low back pain and musculoskeletal injury claims among nurses. Occup. Environ. Med. 2013, 70, 491-497.

4. De Looze, M.P.; Bosch, T.; Krause, F. Exoskeletons for industrial application and their potential effects on physical work load. Ergonomics 2016, 59, 671-681.

5. Abdoli-E, M.; Agnew, M.J.; Stevenson, J.M. An on-body personal lift augmentation device (PLAD) reduces EMG amplitude of erector spinae during lifting tasks. Clin. Biomech. 2006, 21, 456-465.

6. Kawamoto, H.; Sankai, Y. Power assist method based on phase sequence and muscle force condition for HAL. Adv. Robot. 2005, 19, 717-734.

7. Toxiri, S.; Koopman, A.S.; Lazzaroni, M. Rationale, implementation and evaluation of assistive strategies for an active back-support exoskeleton. Front. Robot. AI 2018, 5, 5.

8. ISO 13482:2014. Robotics. Robots and Robotic Devices. Safety Requirements for Personal care Robots; International Organization for Standardization: Geneva, Switzerland, 2014.

9. Waters, T.R.; Putz-Anderson, V.; Garg, A.; Fine, L.J. Revised NIOSH equation for the design and evaluation of manual lifting tasks. Ergonomics 1993, 36, 749-776.

10. van Dieën, J.H.; Selen, L.P.J.; Cholewicki, J. Trunk muscle activation in low-back pain patients, an analysis of the literature. J. Electromyogr. Kinesiol. 2003, 13, 333-351.

11. Johanson, E.; Brumagne, S.; Janssens, L. The effect of acute back muscle fatigue on postural control strategy in people with and without recurrent low back pain. Eur. Spine J. 2011, 20, 2152-2159.

12. Chan, S.T.; Fung, P.K.; Ng, N.Y. Dynamic changes of elasticity, cross-sectional area, and fat infiltration of multifidus at different postures in men with chronic low back pain. Spine J. 2012, 12, 381-388.

13. Nabeshima, C.; Ayusawa, K.; Hochberg, C.; Yoshida, E. Standard performance test of wearable robots for lumbar support. IEEE Robot. Autom. Lett. 2018, 3, 2182-2189.

14. JIS B 8456-1:2017. Personal Care Robots, Part 1: Physical Assistant Robots for Lumbar Support; Japanese Standards Association: Tokyo, Japan, 2017. (In Japanese)

15. Potvin, J.R. Use of NIOSH equation inputs to calculate lumbosacral compression forces. Ergonomics 1997, 40, 691-707.

16. Merryweather, A.S.; Loertscher, M.C.; Bloswick, D.S. A revised back compressive force estimation model for ergonomic evaluation of lifting tasks. Work 2009, 34, 263-272.

17. van Dieën, J.H.; Faber, G.S.; Loos, R.C.; Kuijer, P.P.F.; Kingma, I.; van der Molen, H.F.; Frings-Dresen, M.H. Validity of estimates of spinal compression forces obtained from worksite measurements. Ergonomics 2010, 53, 792-800.

18. Chaffin, D.B.; Andersson, G.; Martin, B.J. Occupational Biomechanical Models. In Occupational Biomechanics; John Wiley & Sons: Chicago, IL, USA, 2006; pp. 133-134.

19. Hof, A.L. An explicit expression for the moment in multibody systems. J. Biomech. 1992, 25, 1209-1211.

20. Kingma, I.; De Looze, M.P.; Toussaint, H.M. Validation of a full body 3-D dynamic linked segment model. Hum. Mov. Sci. 1996, 15, 833-860.

21. Chaffin, D.B.; Baker, W.H. A biomechanical model for analysis of symmetric sagittal plane lifting. AIIE Trans. 1970, 2, 16-27.

22. McGill, S.M.; Norman, R.W. Effects of an anatomically detailed erector spinae model on L4L5 disc compression and shear. J. Biomech. 1987, 20, 591-600.

23. Rajaee, M.A.; Arjmand, N.; Shirazi-Adl, A.; Plamondon, A.; Schmidt, H. Comparative evaluation of six quantitative lifting tools to estimate spine loads during static activities. Appl. Ergon. 2015, 48, 22-32.

24. Marras, W.S.; Davis, K.G.; Kirking, B.C.; Bertsche, P.K. A comprehensive analysis of low-back disorder risk and spinal loading during the transferring and repositioning of patients using different techniques. Ergonomics 1999, 42, 904-926.

25. Endo, Y.; Tada, M.; Mochimaru, M. Dhaiba: Development of virtual ergonomic assessment system with human models. In Proceedings of the 3rd International Digital Human Symposium, Tokyo, Japan, 20-22 May 2014.

26. Imamura, Y.; Ayusawa, K.; Endo, Y.; Yoshida, E. Simulation-based design for robotic care device: Optimizing trajectory of transfer support robot. In Proceedings of the ICORR 2017-15th IEEE International Conference on Rehabilitation Robotics, London, UK, 17-20 July 2017.

27. McGill, S.M.; Norman, R.W.; Cholewicki, J. A simple polynomial that predicts low-back compression during complex 3-D tasks. Ergonomics 1996, 39, 1107-1118.

28. Granata, K.P.; Marras, W.S. An EMG-assisted model of trunk loading during free-dynamic lifting. J. Biomech. 1995, 28, 1309-1317.

29. Jäger, M.; Luttmann, A.; Göllner, R.; Laurig, W. "The Dortmunder"-Biomechanical Model for Quantification and Assessment of the Load on the Lumbar Spine. SAE Trans. 2001, 2163-2171.

30. Nachemson, A.L.F.; Morris, J.M. In vivo measurements of intradiscal pressure: Discometry, a method for the determination of pressure in the lower lumbar discs. JBJS 1964, 46, 1077-1092.

31. Wilke, H.J.; Neef, P.; Caimi, M.; Hoogland, T.; Claes, L.E. New in vivo measurements of pressures in the intervertebral disc in daily life. Spine 1999, 24, 755-762.

32. Sato, K.; Kikuchi, S.; Yonezawa, T. In vivo intradiscal pressure measurement in healthy individuals and in patients with ongoing back problems. Spine 1999, 24, 2468-2474.

33. Takahashi, I.; Kikuchi, S.I.; Sato, K.; Sato, N. Mechanical load of the lumbar spine during forward bending motion of the trunk-a biomechanical study. Spine 2006, 31, 18-23.

34. Putzer, M.; Ehrlich, I.; Rasmussen, J.; Gebbeken, N.; Dendorfer, S. Sensitivity of lumbar spine loading to anatomical parameters. J. Biomech. 2016, 49, 953-958.

35. Bruno, A.G.; Mary, L.B.; Anderso, D.E. Development and validation of a musculoskeletal model of the fully articulated thoracolumbar spine and rib cage. J. Biomech. Eng. 2015, 137, 1-10.

36. Dreischarf, M.; Rohlmann, A.; Zhu, R.; Schmidt, H.; Zander, T. Is it possible to estimate the compressive force in the lumbar spine from intradiscal pressure measurements? A finite element evaluation. Med. Eng. Phys. 2013, 35, 1385-1390.

37. Henninger, H.B.; Reese, S.P.; Anderson, A.E.; Weiss, J.A. Validation of computational models in biomechanics. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 2010, 224, 801-812.

38. Koo, T.K.; Li, M.Y. A guideline of selecting and reporting intraclass correlation coefficients for reliability research. J. Chiropr. Med. 2016, 15, 155-163.

39. Dolan, P.; Luo, J.; Pollintine, P.; Landham, P.R. Intervertebral disc decompression following endplate damage: Implications for disc degeneration depend on spinal level and age. Spine 2013, 38, 1473-1481.

40. Cabello, J.; Cavanilles-Walker, J.M.; Iborra, M.; Ubierna, M.T.; Covaro, A.; Roca, J. The protective role of dynamic stabilization on the adjacent disc to a rigid instrumented level. An in vitro biomechanical analysis. Arch. Orthop. Trauma Surg. 2013, 133, 443-448.

41. Yonezawa, T. Development and clinical application of an intradiscal pressure sensor. Biomed. Eng. 1997, 35, 249-253. (In Japanese)

42. van Dieën, J.H.; Kingma, I. Effects of antagonistic co-contraction on differences between electromyography based and optimization based estimates of spinal forces. Ergonomics 2005, 48, 411-426.

43. Bogduk, N.; Macintosh, J.E.; Pearcy, M.J. A universal model of the lumbar back muscles in the upright position. Ergonomics 1992, 17, 897-913.

44. Kumar, S. Moment arms of spinal musculature determined from CT scans. Clin. Biomech. 1988, 3, 137-144.

45. McGill, S.M.; Patt, N.; Norman, R.W. Measurement of the trunk musculature of active males using ct scan radiography: Implications for force and moment generating capacity about the l4l5 joint. J. Biomech. 1988, 21, 329-341.

46. Nèmeth, G.; Ohlsèn, H. Moment arm lengths of trunk muscles to the lumbosacral joint obtained in vivo with computed tomography. Spine 1986, 11, 158-160.

47. Daggfeldt, K.; Thorstensson, A. The mechanics of back-extensor torque production about the lumbar spine. J. Biomech. 2003, 36, 815-825.

48. McGill, S.M. A myoelectrically based dynamic three-dimensional model to predict loads on lumbar spine tissues during lateral bending. J. Biomech. 1992, 25, 395-414.

49. Cholewicki, J.; McGill, S.M. Mechanical stability of the in vivo lumbar spine: Implications for injury and chronic low back pain. Clin. Biomech. 1996, 11, 1-15.

50. Stokes, I.A.; Gardner-Morse, M. Quantitative anatomy of the lumbar musculature. J. Biomech. 1999, 32, 311-316.

51. Wilke, H.J.; Neef, P.; Hinz, B.; Seidel, H.; Claes, L. Intradiscal pressure together with anthropometric data-a data set for the validation of models. Clin. Biomech. 2001, 16, S111-S126.

52. Kankaanpää, M.; Taimela, S.; Laaksonen, D.; Hënninen, O.; Airaksinen, O. Back and hip extensor fatigability in chronic low back pain patients and controls. Arch. Phys. Med. Rehabil. 1998, 79, 412-417.

53. Kudo, N.; Yamada, Y.; Ito, D. Age-related injury risk curves for the lumbar spine for use in low-back-pain prevention in manual handling tasks. Robomech. J. 2016, 6, 12.

54. Hibbs, A.E.; Thompson, K.G.; French, D.N.; Hodgson, D.; Spears, I.R. Peak and average rectified EMG measures: Which method of data reduction should be used for assessing core training exercises? J. Electromyogr. Kinesiol. 2011, 21, 102-111.