The back exoskeleton we designed has the same structure on the left and right sides of the body (Figure 1), with the overall structure symmetrically arranged along the sagittal plane. The exoskeleton mechanism resembles a lever placed on each side of the body. When the body bends forward, the compression spring gradually compresses, generating two equal but opposite spring forces. The upward spring force is applied to the chest pad, pulling the upper body upward; we refer to this as the chest traction force. The downward spring force is distributed at the A ₃ point, acting as the driving force at the power point of the exoskeleton lever. The resistance force is distributed at the lever’s resistance point A ₄ on the back. If the back support rod is considered a two-force bar, an equal and opposite force is generated at the back pad position, directed upward along the back support rod; we call this the back traction force. The connection joint between the lever and the thigh support plate serves as the pivot point A ₅, allowing the lever to reach a balanced state through the interactions of these elements. The reaction force at the pivot is transferred to the thigh support plate.

A simplified sketch of the human-exoskeleton mechanism (Figure 2) allows to analyze the relationship between the upper body motion and the compression of the spring, with reference to the angles of the lumbar spine and hip joints under standard flexion-extension movements measured by Pal et al. (2007) (Figure 3), synergistic movement patterns of the vertebrae by Wong et al. (2006), to calculate the traction law generated by the exoskeleton under a standard flexion-extension cycle.

An absolute coordinate system is established at the hip joint O ₁ (note: all angles are referenced to this coordinate system), then two closed-loop vector equations are created to describe the mechanism motion (Equations 1, 2): x A 2 = l 6 ⁡ cos ⁡ θ 5 + l 5 ⁡ cos ⁡ θ 4 − l 9 y A 2 = l 6 ⁡ sin ⁡ θ 5 + l 5 ⁡ sin ⁡ θ 4 (1) x A 1 = l 7 ⁡ cos ⁡ θ 6 + l 8 ⁡ cos ⁡ θ 3 − l 9 y A 1 = l 7 ⁡ sin ⁡ θ 6 + l 8 ⁡ sin ⁡ θ 3 (2) where (x _A2, y _A2) and (x _A1, y _A1) are the coordinates of point A ₂, position of the exoskeleton tied to the back, and point A ₁, position of the exoskeleton tied to the chest, respectively, in the absolute coordinate system. To simplify the model, the left and right limbs are assumed to be perfectly symmetrical during flexion and extension movements, so we model the human body as a multi-rigid linkage in the sagittal plane (Manns et al., 2017). Due to the limited range of independent motion of the head and thoracic spine relative to the lumbar spine (Consmüller et al., 2012), we treat them as a single rigid body. The lumbar region is represented as a sequentially articulated mechanism with five vertebrae and the pelvis connected by revolute joints. The elbow and wrist joint angles are fixed, and the arms are assumed to remain vertical throughout the flexion-extension process to simulate the natural hanging posture of the arms during bending and lifting. Since flexion primarily requires small joint angles at the knees and ankles (Bazrgari et al., 2007; Carnegie et al., 2022), the entire lower limb is treated as a rigid body fixed to the ground. Subsequent analyses will disregard the kinematic influence of the lower limbs on the human body and exoskeleton. The whole human body model consists of 8 rigid bodies including arms, thoracic vertebrae, lumbar vertebrae, sacrum, and thighs, with 6 degrees of freedom (lighter lines in Figure 2).

The geometric dimensions and inertia parameters of the vertebrae are based on measurements by Kuo et al. (2010), Pearsall et al. (1996). The coordinates of the points A ₂, A ₁ with the movement of the upper body were obtained using the homogeneous coordinate transformation method, after which the variation of the spring length could be calculated (Equations 3, 4): l 8 = x A 1 + l 9 − l 7 ⁡ cos ⁡ θ 6 2 + y A 1 − l 7 ⁡ sin ⁡ θ 6 2 θ 3 = y A 1 − l 7 ⁡ sin ⁡ θ 6 x A 1 + l 9 − l 7 ⁡ cos ⁡ θ 6 (3) θ 5 = arcsin h 0 ⁡ sin ⁡ φ h 1 − arccos h 1 2 + l 6 2 − l 5 2 2 h 1 l 6 θ 4 = arcsin h 0 ⁡ sin ⁡ φ − l 6 ⁡ sin ⁡ θ 5 l 5 φ = arctan y A 2 x A 2 h 0 = x A 2 2 + l 9 2 + y A 2 2 h 1 = h 0 2 + l 9 2 + 2 h 0 l 9 ⁡ cos ⁡ φ θ 6 = θ 0 + θ 5 (4)

Consider an ideal compression spring that, with damping and preload force neglected, produces a chest traction force (spring force) of (Equation 5): F f = k s l s 0 − l 8 (5) where k _s is the sum of the spring stiffnesses of the exoskeleton on both sides and l _s0 is the initial length of the spring of the human body in the standing phase. Assuming that the exoskeleton consists entirely of lightweight rods, we can derive the static equilibrium equations for the lever. This allows us to determine the mathematical relationship between the two types of traction forces (Equation 6): F b = F f l 7 ⁡ sin θ 6 − θ 3 l 6 ⁡ sin θ 5 + θ 4 (6) where F _b is the back traction force. The mathematical model predicts the trend of the two traction forces generated by the exoskeleton under standard flexion and extension movements.

The moments for the body due to chest and back traction are calculated by the cross-multiplication of the vectors from points A ₁ and A ₂ to the L5/S1 joints and the vectors of the chest and back traction forces, and the sum of the two moments is the total moment of the exoskeleton.

According to the traction force model, there is a direct relationship between exoskeleton spring stiffness and traction force. The amount of traction force affects the body’s muscle force generation. Too little traction force can lead to insufficient assistive effect, and the disc pressure will remain at a high threshold, which is not effective in avoiding spinal injuries. Too much traction force may lead to negative muscle activity, which in turn may cause injury to the human body. By establishing a musculoskeletal-based kinetic model, the effect of traction force on the musculoskeletal system under the phenomenon of flexion relaxation can be simulated, and the decompression performance of the exoskeleton with optimal traction force can be obtained.

Due to the complexity and large number of low back muscle groups, we simplify the calculations by modeling only the erector spinae muscle as the main force-generating muscle group in flexion and extension, with a total of 22 muscle fibers. Within this muscle group, muscle fiber bundles are divided into the thoracic longissimus (LT), lumbar iliocostalis (LL), and lumbar longissimus (LS) based on the different attachment points of the muscle fibers (Figure 4). In contrast, the rectus abdominis (RA) was also modeled to represent the antagonist muscle of the longest muscle of the thoracic spine. Muscle contraction dynamics are described using the hill model (Millard et al., 2013), expressed as follows (Equation 7): F m = F o m a m F a m + F p m (7) where, F _m is the total muscle force, F _am is the dominant force factor, F _pm is the passive force factor, F _om is the maximum isometric force of the muscle fiber, a _m is the muscle fiber activation.

The posterior longitudinal ligament (PLL), ligamentum flavum (LF), supraspinous ligament (SSL), and interspinous ligament (ISL), which play a primary role in the flexion-extension movement, were selected (Figure 4). The attachment points of these four lumbar spine ligaments on the vertebrae were determined according to the method described by Naserkhaki et al. (2018). The ligament stiffness data followed the force-deflection curves determined by Pintar et al. (1992) for ligaments of the lumbar T12-L2, L2-L4, and L4-Sl segments. For computational convenience, nonlinear stiffness curves for all ligaments were fitted with a sixth-degree polynomial, which when multiplied by the ligament stretch can be expressed as the ligament force F _l . The intervertebral disc joint was modeled as a torsion spring with constant stiffness in the coronal axis to calculate the disc’s resistance moment M _v , and their stiffness data were obtained from (Senteler et al., 2016).

The Newton-Euler method was employed to establish the multi-body dynamics equations for each rigid body (Equation 8), which can visualize the contact forces between different vertebrae, as well as the contributions made by active and passive muscle forces, ligaments, and intervertebral discs to the generation of extension moments. This, in turn, facilitates the assessment of intervertebral disc pressures. The L5/S1 intervertebral disc pressure, which is a combination of disc moments, muscle forces, ligament forces, and vertebral joint forces F _x , F _y , can be written as: F c 5 = F F m , F l , F x , F y , M v (8)

Since the unknown amount of muscle fiber activation is much larger than the number of equations in the kinetic model, optimization algorithms are required to obtain solutions, while different objective functions lead to variability in the optimization results (Ostraich and Riemer, 2024). We formulate the muscle as well as the stiffness solution problem as a multi-objective optimization problem (Figure 5) with constraints based on the dynamical equations of Newton’s Euler method and a range of muscle activation a _m . Considering that with a spring stiffness of 2000 N/m the exoskeleton traction force will be close to the effective traction force threshold of 500 N (Zaïri et al., 2021), the stiffness does not exceed this value. It is constructed as follows (Equations 9, 10): min   f a m , k s = f 1 , f 2 , f 3 (9) s . t .   m a = ∑ F x + ∑ F y + ∑ F m + ∑ F l J θ ¨ = ∑ F x × l x + ∑ F y × l y + ∑ F m × l m + ∑ F l × l l + ∑ M v 0 ≤ a m ≤ 1 (10) where m represents the mass, J represents the rotational inertia, a represents the acceleration, and θ ¨ represents the angular acceleration of the joint. l _x , l _y , l _m , l _l are the force arms of each component of the force relative to the disc rotating the joint, respectively.

The first objective function is the minimization of the sum of squares of muscle activation (Equation 11), which means that the sum of work done by human muscles is minimized. It is widely applied in various muscle optimization solution problems (Rahmati et al., 2017; De Groote et al., 2016; Manns et al., 2017). Unlike other methods, the optimization values of muscle activation in this function are related to the corresponding maximal isometric force, indicating that muscle fibers with greater force capacity will tend to exert more force. f 1 = ∑ i = 1 22 200 a m , i 2 F o m , i (11)

The second objective function employs the total sum of squared differences to minimize the variation in activation levels among all muscle fibers (Equation 12). This reflects the body’s ability to neurally marshal muscle, with similarity in activation changes of muscle fibers that play the same role when muscle groups are stimulated through the same or approximate reflex pathway. f 2 = ∑ i = 1 22 ∑ j = i + 1 22 a m , i − a m , j 2 (12)

The third objective minimizes the cumulative reaction force and moment at the L5/S1 intervertebral disc over a flexion-extension cycle (Equation 13), with w ₁, w ₂ and w ₃ representing weight coefficients, F _L5/S1, M _L5/S1, and M _hip representing the L5/S1 force, the L5/S1 moment, and the hip moment, respectively. This reflects the body’s tendency to generate active muscle forces in a way that minimizes skeletal load (Brown and Potvin, 2005). Although the optimization principles of these three objective functions may conflict, they together capture the true activation patterns of the erector spinae muscles during bending and lifting movements. Including exoskeleton spring stiffness in the optimization framework does not shift the optimization direction; on the contrary, it represents an optimal stiffness that reduces spinal load and muscle activation without compromising muscle coordination. f 3 = w 1 F L 5 / S 1 + w 2 M L 5 / S 1 + w 3 M h i p (13)

The simulation framework was constructed and executed in MATLAB 2022b, with optimization conducted using a multi-objective genetic algorithm (MOGA) (Deb, 1999). MOGA balances multiple objective functions by finding a set of non-dominated solutions (Pareto optimal solutions) to achieve the best compromise between different objectives. The L5/S1 intervertebral disc pressure before wearing the exoskeleton was solved by assigning the degree of muscle activation through a multi-objective genetic algorithm, so as to assign a reasonable muscle force to each joint to satisfy the equilibrium of kinetic states. The spring stiffness was optimized according to the peak pressure of L5/S1 joints, and the L5/S1 disc pressure after wearing the exoskeleton was optimized by incorporating the traction force. The simulation timestep was set to 0.1 s, with the model starting in an initial upright posture, bending to the lowest point, then returning to an upright stance, covering a total duration of 4 s. This setup allows the simulation to capture real-time loading on intervertebral discs and changes in activation levels across individual muscle fibers of the erector spinae. For ease of analysis, the muscle fibers were grouped into three major muscle groups—thoracic longissimus, lumbar iliocostalis and longissimus, and rectus abdominis. The activation level of each muscle group was represented by the sum of activations of the fibers within that group.

A total of eight healthy male subjects were invited to participate voluntarily in the trial, with a mean age of 24 ± 1.5 years, mean height of 176 ± 8.3 cm, and mean weight of 73.3 ± 5.8 kg. All subjects had no report of low back pain in the last 3 months and had not been involved in large lifting tasks. Before the experiment, all subjects read the instructions and precautions related to the experiment and signed an informed consent form for the experiment. The experiment was approved by the Ethical Review Committee of Hebei University of Technology.

Raw EMG data were recorded using a Noraxon Ultium Wireless Surface EMG Tester (Noraxon Inc., United States) with a sampling frequency of 2000 Hz. The skin of the subjects was shaved before the experiment, and then cleaned of grease and other contaminants by wiping with an alcohol cotton ball. The surface EMG sensors were attached to the thoracic erector spinae (TES), lumbar erector spinae (LES), and rectus abdominis (RA), of which the thoracic erector spinae and lumbar erector spinae were mainly responsible for trunk flexion and extension, and the rectus abdominis muscle was responsible for trunk flexion. Before the experiment, the sensor light was covered with black tape to prevent the motion capture system from misjudging it as a Marker point. The kinematic data of the joints were recorded using a Vicon T40-S Motion Capture System (Oxford Metrics Ltd., United Kingdom) with a sampling frequency of 200 Hz. Reflex point markers were attached to T7, T12, S1, PSIS, ASIS, hip, mid-thigh, knee, mid-calf, ankle, and tip of the foot. The angle of the vector formed between the labeled points was used as the angle of joint movement during the experiment.

The experiment was divided into four groups: No weight in hands without exoskeleton, 8 kg weight in hands without exoskeleton, no weight in hands with exoskeleton, and 8 kg weight in hands with exoskeleton. During the complete experiment, subjects were asked to keep their arms parallel to the direction of gravity at all times. Without the exoskeleton, the subjects were bending from the initial standing position in a natural descending motion to the subjective conscious limit, and then returned to the standing position in a natural rising motion without artificial intervention of synergistic movements between lumbar vertebrae and hip joints. There is no artificial interference with the synergistic movement between the lumbar spine and the hip joints. The phases from standing to bending, and from rising to standing, were all approximately 2 s. The choice of 0 kg was made in view of the fact that even without a load, the disc pressure during bending is not optimistic, and back pain in farmers has been reported in the field of agricultural harvesting or cultivation (Fathallah, 2010). 8 kg was chosen because it is the load level at which the NOISH lifting criteria for L5/S1 disc compression can reach the injury risk threshold (McGill and Kippers, 1994).The two loads represent most of the typical light-duty load-bearing work situations. Subjects wore the exoskeleton and started the experiment with a bending and rising motion without the exoskeleton as the norm, causing a supervisor to observe the subject’s posture in real time and to remind and correct his movements, using a metronome to control the tempo of the task. Subjects were made to not bend their knees and ankles as much as possible. If the knee and ankle flexion angles exceeded a certain threshold, the data from this set of experiments were nullified. This is because if the knee and ankle joints were overflexed it would affect the phenomenon of flexion relaxation (Shin et al., 2004),and would not be conducive to the observation of the connection between the exoskeleton and the EMG signals. All experimental situations were repeated six times, and after completing one set of experiments, the subjects were allowed to take a full rest for 10 min to prevent muscle fatigue.

The original EMG signals were first imported into MATLAB, and the raw data were filtered using a fourth order Butterworth bandpass filter. Afterwards, the data were subjected to full-wave rectification and the RMS values of the EMG data were plotted. To ensure the comparability of electromyography (EMG) data among different subjects, the EMG signal data were normalized. Prior to the experimental tasks, the maximum voluntary contraction (MVC) of each subject was measured using the manual resistance method on a Roman chair (Biviá-Roig et al., 2019). The experimental signal values of each group were divided by the MVC signal values to obtain the normalized data, expressed as MVC%. This normalization allows signals from different individuals or different experimental conditions to be adjusted to the same scale, enabling direct comparisons between signals from different datasets or conditions and facilitating subsequent processing and analysis steps. The Kolmogorov - Smirnov test was employed to examine the normality of the electromyography (EMG) signal parameters. Subsequently, a paired t - test (α = 0.05) was utilized to evaluate whether there were significant differences in the peak values of the MVC% signals before and after wearing the exoskeleton.

Simulation results showed that the total stiffness of the optimized spring under 8 kg load was 1233 N/m. The chest and back traction provided by the exoskeleton increased with increasing flexion (Figure 6). The exoskeleton provided almost no traction and moment at the 0%–15% and 85%–100% stages. This is because at this stage the lumbar joints are less than 22.1°, the hip joints are less than 14.2°, and the angle of inclination of the human trunk is small enough that no more traction or extension moments are required. The exoskeleton provided maximum traction when the body reached 50% of the final stage of flexion, with thoracic traction of 173.4 N and dorsal traction of 125.5 N. The flexion moment generated by the dorsal traction partially offset the extension moment generated by the thoracic traction, making the overall extension moment generated by the exoskeleton appear not to be excessive (Figure 6). At the point where the body reaches 50% of the final stage of flexion, the maximum extension moment is just 13.1 N m. Lower support torque and higher traction performance may be better adapted in workplaces oriented to bending up under light loads or 0 kg loads (Zhang et al., 2023). Multiple traction forces demonstrated an excellent mechanical effect of reducing disc pressure (Figure 7), with the main reduction reflected in the peak phase of disc pressure, where peak compression forces were reduced from 3,326.4 N to 2,873.2 N.

The exoskeleton assist mechanism met the expectations of the objective function planning. At 8 kg weightbearing, muscle activation before and after wearing the exoskeleton exhibited the following changes (Figure 8): during the 0%–25%, 50% and 75%–100% phases of the flexion-extension exercise cycle, which are phases in which the lumbar spine is in mild flexion or deep flexion, the muscle activation had a similar trend with or without wearing the exoskeleton, consistent with the lack of active muscle force during the flexion-relaxation phenomenon. Since the traction force of the exoskeleton can also generate a partial extension moment, a trend of reduced muscle activation after wearing the exoskeleton was observed during the partial flexion-extension movement interval. Simulation results showed that the activation of the rectus abdominis muscle did not change much before and after the subjects wore the exoskeleton, and never exceeded the upper limit of 0.2. The experimental results (Figure 9) showed that, in the 0 kg experimental group, there was no significant difference in the level of EMG signals in the thoracic erector spinae muscle after wearing the exoskeleton (p = 0.62). The peak EMG signals of both the lumbar erector spinae and rectus abdominis muscles were slightly reduced, by 6.41% in the lumbar erector spinae (p = 0.036) and by 18.08% in the rectus abdominis (p = 0.02). In the 8 kg experimental group, the level of electromyographic signals in the thoracic erector spinae muscle after wearing the exoskeleton was also not significantly different (p = 0.226). Both the lumbar erector spinae and rectus abdominis muscles also had a slight decrease in peak EMG signals, with the lumbar erector spinae decreasing by 7.32% (p = 0.015) and the rectus abdominis decreasing by 16.03% (p = 0.023). The peak EMG signal of the rectus abdominis was consistently much lower than that of the erector spinae in both the 0 kg and 8 kg experiments and no significant elevation of the peak EMG signal of the rectus abdominis was observed with the exoskeleton. Because of the low extension moments generated, the results of both simulations and experiments observed a low impact of the exoskeleton on the muscles, which is consistent with the original design intent that the traction generated by the exoskeleton is more focused on reducing disc pressure than on reducing muscle activation.