Abstract

Training based on muscle-oriented repetitive movements has been shown to be beneficial for the improvement of movement abilities in human limbs in relation to fitness, athletic training, and rehabilitation training. In this paper, a muscle-specific rehabilitation training method based on the optimal load orientation concept (OLOC) was proposed for patients whose motor neurons are injured, but whose muscles and tendons are intact, to implement high-efficiency resistance training for the shoulder muscles, which is one of the most complex joints in the human body. A three-dimensional musculoskeletal model of the human shoulder was used to predict muscle forces experienced during shoulder movements, in which muscles that contributed to shoulder motion were divided into 31 muscle bundles, and the Hill model was used to characterize the force-length properties of the muscle. According to the musculoskeletal model, muscle activation was calculated to represent the muscle force. Thus, training based on OLOC was proposed by maximizing the activation of a specific muscle under each posture of the training process. The analysis indicated that the muscle-specific rehabilitation training method based on the OLOC significantly improved the training efficiency for specific muscles. The method could also be used for trajectory planning, load magnitude planning, and evaluation of training effects.

1. Introduction

Rehabilitation robots that provide rehabilitation therapy following neurological injuries, such as stroke and spinal cord injury, have received increasing interest [1]. The shoulder complex is one of the most complicated joints in the human body, as it directly affects the performance of the whole upper limb movement, including hand manipulations. After suffering from a neurological disease, the motor function of the human shoulder complex can often face deterioration, and there is an increased risk of spasticity; thus, a proper approach to motor function training therapy is necessary. Previous studies have also indicated that the use of rehabilitation robotics can provide a high-intensity, task-oriented, and highly repetitive treatment in the impaired upper limb, which have all been shown to be beneficial for the restoration of shoulder function [28]. Rehabilitation robots can train patients’ limb movements in several ways by applying turnable forces to patients. Training with active-resistance mode rehabilitation robots can actively deliver resistance against movements executed by the patient. Accurate rehabilitation for human limb motor function has been widely accepted and has become increasingly popular; accurate rehabilitation requires both a higher level of control and smarter robotic designs as well as increasingly accurate models of the human motor system. However, the rehabilitation efficiency of current methods is still controversial [911] because training methods are often based on the clinical experience of the doctor or therapist who chooses from standard menus and options rather than making precise decisions based on the actual situation present in the patients’ limb [10, 12]. By contrast, most studies of rehabilitation robotics are focused on rehabilitating a whole limb, so loads often cannot be applied to a specific muscle, which is important in the rehabilitation training process because different muscles often have different levels of damage, stiffness, or deterioration. To apply efficient and targeted training to a specific muscle, Terashima et al. and Itokazu et al. [1315] proposed the concept of specific-muscle training for the upper limbs using optimization control of patients’ neural network and via EMG feedback control methods. The results showed that this training method aimed at specific muscles maximized the training effect in target muscles and simultaneously weakened the effects on other muscles. However, their research was limited to cases of small-scale movements within a single horizontal plane, and the need for real-time EMG signals as an input value has limited this method’s practicality.

Compared with most developers of rehabilitation robotics, researchers in biomedical engineering have focused on studying the physiological musculoskeletal model of the human shoulder. For shoulder models, the typical representations of the muscles’ lines of action are the line-segment model [11, 1618] and more complex 3D finite element model [19]. These models were intended to be used for several purposes, such as surgical simulation [17], wheelchair mechanics research [20, 21], neuroprostheses control [22, 23], and other similar applications. Most current studies using a musculoskeletal model focus on joint-contact forces and muscle moment arms. The direct purpose of these studies was to reproduce and simulate the patterns of muscle force generation. However, musculoskeletal models based on the anatomical structure have not been adequately or currently applied to rehabilitation robots [10, 24]. Therefore, in this paper, a three-dimensional mathematical musculoskeletal model [16, 25, 26] of the shoulder complex was used to represent the relationship of muscles and external loads. Based on this model, the static forces that acted on all of the shoulder bones, including the humerus, scapula, and clavicle as well as the muscles connected to these bones, were analysed and calculated. Muscle activation was proposed to measure the force exhibited by a certain muscle undergoing an isometric contraction. The optimal load orientation (OLO) for muscle rehabilitation was then obtained according to the results of the calculated muscle activation. Finally, a specific muscular rehabilitation training approach based on the optimal load orientation concept (OLOC) was proposed.

Rehabilitation robots were controlled to apply forces to each joint axis separately via intelligent application of the exoskeleton mechanism. For a given rehabilitation movement trajectory of the shoulder, the optimal load orientation of a specific muscle of the whole range of motion (optimal load orientation cluster) was determined by calculating an inverse dynamic problem. Then, the shoulder movement trajectory of rehabilitation was designed or evaluated according to the optimal load orientation cluster above. The results showed that training based on the OLOC enhanced the activation of specific muscles and reduced the activation of other muscles and thus enabled efficient training of specific muscles.

This paper also assessed the influence of the magnitude of the external load on the training effects of OLOC training through simulation testing.

2. Materials and Methods

When implementing active-resisted rehabilitation training for an impaired shoulder using a robotic device, robots delivered resistance against active movements executed by the shoulder.

During shoulder movement, for a specific muscle, a certain external load with different orientations leads to different muscle forces. Therefore, for a specific muscle, muscle activation can be designed by controlling the orientation of the external load.

2.1. Musculoskeletal Model of the Shoulder Complex

A three-dimensional mathematical musculoskeletal shoulder model [16, 25, 26] was used, as shown in Figure 1 [27], to represent the geometric architectural properties of the skeleton and muscles of the shoulder. The geometric parameters of the model were developed using CT images of bones and muscles collected from the Visible Human Project (VHP) database [25].

2.1.1. Skeleton Model

The shoulder skeletal structure consists of the following bones: the thorax, clavicle, scapula, and humerus, and the following joints: the sternoclavicular joint (SC), acromioclavicular joint (AC), glenohumeral joint (GH), and scapulothoracic joint (ST).

In this paper, a hybrid mechanism model was used to simulate the structure of the skeletal system of the shoulder, as depicted in Figure 2, where “1” represents the clavicle, “2” represents the scapula, “3” represents the thorax, and “4” represents the humerus; point a represents the SC articulation; b represents the AC articulation; c and d represent the upper feature point and lower feature point of the scapula, respectively; and e represents the midpoint of the EL (lateral epicondyle) and EM (medial epicondyle) in the humerus.

Considering that translations are negligible compared with rotations, the SC, AC, and GH joints were assumed to be ball-and-socket joints. As an exception, the ST joint was considered to be a joint that allowed the scapula translation and rotator movement with respect to the thorax due to the compliance of the surrounding muscles. The thorax was represented as an ellipsoid, as shown in Figure 2.

To describe and analyse the skeleton model above, every bone is fixed to a coordinate system, as shown in Table 1. The ISB standard recommended by Wu et al. [28] is widely used in the field of biomedical engineering to describe the movements of bones and joints in the human upper limb. The relationship between each coordinate system described in Table 1 and its corresponding coordinate system recommended in the ISB standard can be described by a rotation transformation matrix.

2.1.2. Muscle-Driven Model and Muscle Activation

The eighteen muscle groups associated with shoulder motion were divided into 31 muscle bundles according to the results of anatomical measurements of shoulder muscles performed by Garner and Pandy [26], and these muscle bundles were numbered from M1 to M31, as shown in Table 2.

To calculate the muscle force, each muscle bundle was modelled as a 3-element Hill-type model that was widely applied in muscle-driven simulations [2932]. Four parameters were used to represent each muscle’s force-generating properties, including the tendon slack length (), pennation angle (), optimal muscle-fibre length (), and peak isometric muscle force (). The values of the parameters of each muscle bundle are shown in Table 2 and were determined from the reports by Garner and Pandy [33] and Yanagawa et al. [34].

According to the Hill model, the actual muscle-fibre length () and isometric muscle force () can be calculated using the data presented in Table 1, and the total muscle length of the muscle bundle was determined in the musculoskeletal model by the obstacle-set method proposed by Garner and Pandy [26].

During shoulder movements, the muscle length changes with bone movements. When the shoulder is in a certain position and posture, the muscles can exhibit contractions that are considered isometric, so the actual muscle force can change while the muscle length remains unchanged. Muscle activation is defined by the equation below to describe the force state of a muscle during an isometric contraction:

When the actual force is maximum, and ; when the actual force is minimum, and .

2.1.3. Musculoskeletal Model

Inside the human body, muscles are often observed to be bar-like fibres and are connected to the bones via tendons. Muscles will always bypass some bones, joints, and surrounding tissue that forms the muscle path that passes the origin, via point, obstacle, and point of insertion [26]. Based on the skeletal and muscle models depicted above, the muscle path was determined using the obstacle-set method proposed by Garner and Pandy [26]. Here, the data used in the obstacle-set method, such as the position of the feature points (origin point or insertion point) and type and size of the obstacles, were determined from the results reported by Garner and Pandy [16, 26], and all data were transformed into the coordinate system that is described in Table 1. Once the path of a muscle bundle was determined, its total muscle length was determined. The actual muscle-fibre length () and isometric muscle force () were then calculated.

2.2. Static Analysis and Prediction of Muscle Activation

Muscle activation was used to describe the force condition of the muscle bundles. In this paper, inertial forces were ignored because the motion of the shoulder is slow; thus, the shoulder could be considered to be in a situation of static balance. Some multisolution static equilibrium problems were solved to calculate the muscle forces for each prescribed posture of the shoulder under a certain external load. Any shoulder bone is simultaneously affected by gravity, joint force, joint torque, and muscle torque. The exerted torque of a muscle on any shoulder bone could be determined by static analysis of the bone. Once the muscle torques were obtained, the static balance equations describing the muscles were established and the activation of each muscle bundle was calculated by solving the equations.

2.2.1. Static Analysis of Bones

The gravity of the bones and muscles of the upper limbs was considered to be one force, while the mass and centroid position of the upper limbs were determined based on the results of anthropometrical data from the study by Shan and Bohn [35]. When the shoulder is in a certain posture, the static balance equations describing the humerus, scapula, and clavicle can be established, and then, the muscle input torques that are needed to drive each bone can be calculated by solving the equations.

The humerus was analysed as an example. To indicate the force between two bones, was used to represent the force that is applied to bone by bone (Figure 2). The free body diagram of the upper limb bones is shown in Figure 3, where is the gravitational force on the upper limb, is the external load, is the force applied to the humerus by the scapula, is the position vector of the upper limb centroid, is the positioning vector for the point where the external force is applied, and is the input torque applied to the upper limb produced by the muscles. Since the GH joint was assumed to be a ball-and-socket joint, the joint torque applied to the humerus by the scapula through the GH joint was described as zero.

The value of each vector was determined in the global coordinate . The joint force and the input torque were calculated by solving the static equilibrium equations

Similarly, static analyses of the scapula and the clavicle were completed using the same methods, which are shown in Figures 4 and 5.

The static analysis described above gave the input torques , , , , , and that were needed to drive the bones when the shoulder was in a specific position and posture was under a certain external load, and all of the input torques were generated by shoulder muscles that are connected to the bones.

2.2.2. Static Analysis of the Muscles

The input torque was produced entirely by the muscle bundles attached to the bones, so nine equations were established based on the torque balance conditions of all 31 muscle bundles.

The literature shows that only when the real muscle-fibre length is more than 1.5 times the optimal muscle-fibre length will the muscle lose its active function of contraction and begin to generate passive force [33]. In this paper, the muscles were in the normal stretching range and in a state of active contraction, so the muscle force was always a tensile force. Since the muscle paths were determined already in the musculoskeletal model using the obstacle-set method, the muscle line of action could be determined as well. Therefore, the direction of the muscle force acting on a bone is always from the feature point (origin point or insert point) to the nearest via point.

Figure 6 shows the condition of the muscle forces on the clavicle. The muscle bundles connected to the clavicle are M1, M5, M14, and M20; among these, M14 and M20 are connected to the clavicle at the origin point, while M1 and M5 are connected to the clavicle at the insertion point. In Figure 6, is the activation of the muscle bundle Mi; is the isometric muscle force of Mi; is the vector from point to the origin point of Mi, is the vector from point to the insert point of Mi; is the muscle force vector when Mi is connected to the bone by the origin point, is the muscle force vector when Mi is connected to the bone by the insert point, and and can be determined using the musculoskeletal model. Therefore, the actual muscle force of Mi can be represented as when the feature point of Mi was the origin point (M14, M20) or when the feature point of Mi was the insertion point (M1, M5). The input torque of the clavicle, which had already been calculated by the static analysis of the bones, was entirely produced by the muscle bundles above, so the torque balance equation could be given as follows:

The muscle bundles connected to the scapula could be analysed through the same process, as shown in Figure 7. The muscle bundles connected to the origin point were M21–M31, while the muscle bundles connected to the insertion point were M2, M3, M4, and M6–M13. The input torque was entirely produced by these muscle bundles, so the torque balance equation could be given as follows:

The muscle bundles M14–M31, which are connected to the humerus, were analysed, as shown in Figure 8, where was the adjacent feature point connected to the humerus. All of the feature points of the humerus were at the insertion point because the humerus is located at the end of the shoulder bones. The input torque was entirely produced by these muscle bundles so the torque balance equation of the muscles could be given as follows:

By solving (4)–(6) simultaneously, nine equations were established to calculate muscle activation () in all 31 muscle bundles, so these problems had multiple solutions. According to the studies of Crowninshield and Brand [36], minimizing the sum of the squares of all muscle stresses was chosen as the objective function, so the multisolution problem above was transformed into an optimization problem with certain boundary constraints. The procedure of the static optimization was implemented using the FMINCON function in Matlab, and the optimal objective function was given as follows, where is the physiological cross-sectional area of the muscle bundle Mi given in Table 1:

The real-time physiological cross-sectional area of a muscle bundle always changes during its contraction process, while the muscle belly volume () can be considered to be constant. A more appropriate optimal objective function, therefore, could be given as follows, where represents the actual muscle-fibre length:

2.3. Algorithm of Optimal Load Orientation (OLO)
2.3.1. Definition of the Optimal External-Load Orientation

The diagram of the upper limb under an external load is shown in Figure 9. The force coordinate system , the origin of which was located at the midpoint of the EL and EM, was parallel to the humerus’ frame of reference , which was defined in [25]. The upper arm and forearm were assumed to be relatively static, and the external load that was considered to be a pure force was assumed to be applied on the origin of the coordinate system . During shoulder movement, only when the direction of the external load was perpendicular to the long axis of the humerus would the shoulder muscle suffer the greatest load effect. Therefore, the external load, the orientation of which was described by the angle , was assumed to be in the normal plane and was considered to be perpendicular to the long axis of the humerus.

If a shoulder was moving under an external load with a constant magnitude, the muscle forces would depend on the posture of the shoulder and the orientation of the external load. Therefore, for a specific muscle to be active in a certain posture of the shoulder under a constant-sized external load, there would always be a specific external-load orientation that would lead to maximal muscle activation. This orientation of the external load was defined as the optimal external-load orientation ().

The curve shows the relationship between the level of muscle activation and the orientation of an external load. For example, the curve of a deltoid acromial activation (DLTa) with an external-load magnitude of 2.0 kg (approximately 20 N) is depicted in Figure 10 (where an abduction of 70 deg is performed). The activation plotted with a black line was the result of low-pass filtering (using a moving average filter with the span of 15). [] was the load orientation interval that located the activation level in the interval range from 90% to 100% of its maximum. Thus, the midpoint of the interval [] was defined to be the optimal load orientation of DLTa.

2.3.2. Algorithm Expressing the Optimal External-Load Orientation

For a certain posture of the shoulder, this calculation results in that varies with the magnitudes of in several cases. For some shoulder postures, when the assumed value of was small, the curves of some muscles would appear to be impulse noise around its maximum.Due to the small value of activation and narrow bandwidth of , calculated activations less than 0.02 were considered unreliable and discarded. For instance, the activation level of the triceps brachii long (TRCl) head tended to be very low (<0.02), with an abduction from 20° to 40°, with an external load that was less than 2.5 kg (approximately 25 N). With the external load value increased to be greater than 2.5 kg, the activation of TRCl became significant and an optimal orientation cluster could be determined. Similar patterns happened with many other muscles, such as the supraspinatus, teres major, and deltoid scapular. Additionally, when the load value was assumed to be large, a well-distributed relatively high activation level might appear in the curves, thus resulting in uncertainties relating to . Different load magnitudes could result in different distributions of . The analysis results showed that changed with external-load magnitude but did not change much in certain intervals of external-load magnitudes (0.5~6.0 kg in this paper), which means that the distribution of has significant consistency over a range of loads. To obtain more complete results, was calculated under different external loads ( = 0.5 kg, 1.0 kg, … , 6.0 kg), and the least squares-fitting result of values were calculated as the mean optimal external-load orientation ().

The flow chart shown in Figure 11 describes the algorithm used to determine the optimal external-load orientation. The curve of a specific muscle under a certain magnitude external load for a certain shoulder posture was obtained through iteration, where the value of could be simultaneously determined. Another iteration was then performed by changing the load magnitude from 0.5 kg to 6.0 kg to calculate . Finally, another iteration was completed to obtain all of the values of present in the whole trajectory of a rehabilitation movement ( cluster) by changing the shoulder posture. The clusters of all of the muscle could be similarly determined.

2.3.3. OLOC Training Simulation and Evaluation of the Training Trajectory

The weight-lifting abduction and abduction with OLOC under the same magnitude of external load were simulated to evaluate the effect of OLOC presented above to allow the promotion of the training efficiency. The mean activation over different abduction angles under a certain magnitude of external load was used to represent the training efficiency in a specific muscle.

For a given rehabilitation movement trajectory, the training efficiency differences among different muscles were evaluated, and the applicability of the rehabilitation movement for different muscles was also evaluated.

2.3.4. Influence of the Load Magnitude

Different muscles’ training effects are different in response to different load magnitudes. To determine the response characteristics of the external-load magnitude of the shoulder muscles, the mean activations of each main muscle bundle under different load magnitudes were calculated and compared.

To obtain the optimal load orientation cluster and study the influence of the load magnitude, shoulder abduction was simulated in the coronal plane from 20 to 80° with load magnitudes ranging from 0.5 kg to 6.0 kg, and the results were analysed.

3. Results

3.1. Optimal External-Load Orientation

Despite different distributions of due to different load magnitudes for a specific muscle, the results indicated a significant consistency in the activations calculated under different magnitudes of external load. Figure 12(a) shows the of TRPc, DLTc, and DLTa under a series of load magnitudes. The grey circle shows the distributions of , and the solid line shows the curve obtained by regression fitting. The analysis suggested that within a certain range of loads, the distribution of in these muscle bundles showed good consistency and the values of were unique. The 3D representation of the results is shown in Figure 12(b), where the coordinate system used in this figure coincided with the ground reference described in [25]. The GH centre trajectory (blue line), elbow movement trajectory (pink line), and virtual long axis of the humerus (black dotted line) are all vividly shown. The red line and arrow represent the distribution of . The distribution of has good continuity with the path of the elbow in space.

Like the three muscle bundles described above, many other muscles’ are distributed in a single region, such as PMJs, DLTs, and SUPR. However, except for scenarios involving a single region, there were also multiregion distributions for several muscles, including SRAm, TRPt1, LTDt, TMJ, and TRClg, among others. Figure 13 shows the optimal load orientations of TMJ, where L1 and L2 were the two optimal orientation paths available for selection obtained by regression fitting. In this range of abduction, the muscle activation of TMJ caused by an external load with the same magnitude along L1 or L2 had no large difference (the difference is less than 10%). Moreover, the activations of TMJ performed uniformly during an abduction from 20 to 32 degrees (region L3).

3.2. Active-Resisted Abduction with Loads Based on OLOC

Figure 14 shows a comparison of the mean activations for the main muscles of the shoulder obtained from a simulation of abduction with weight lifting of 3.0 kg that would always apply a vertical force on the humerus and abduction with OLOC under the same magnitude of load (). The results suggested that the muscle-specific rehabilitation training method based on the OLOC significantly promoted the training efficiency of specific muscles (the average of the increased proportion of the mean activation of all muscles was 537%, and the mean activation of all muscles was promoted by 165%).

There were several muscles whose mean activations were promoted more significantly, which included the SRA (from 0.00 to 0.21), LTDc (from 0.02 to 0.26, promote 1640% relatively), TRP (0.04 to 0.37, 905%), TMJ (0.11 to 0.67, 532%), LTDi (0.17 to 0.78, 361%), DLTa (0.13 to 0.57, 325%), and SUPR (0.25 to 0.97, 280%). There were some muscles whose promotion was only moderately increased but was kept at a relatively high level, including the RMN (75%), RMJ (99%), LTDt (95%), and LTDl (68%). However, there were also instances of an insignificant promotion of the mean activation, including the INFR (1.07%) and BICl (5.89%). The results of this comparison indicated that there were differences in the training efficiency among many muscles for the same rehabilitation movement using the OLOC.

3.3. Evaluation of Rehabilitation Training Movements

In an OLOC training movement of a certain trajectory, the mean activations of different muscles were different. The differences in activation were governed by the movement trajectories used in the shoulder’s rehabilitation. Figure 15 shows the average activation of different muscles in OLOC training (abduction from 20 to 80 degree, ). The average activation is calculated by the mean activation and standard deviation (), which is shown as error bars.

On average, there were several muscles whose activations were relatively high (), which included the RMN (), PMN (), LTDt (), LTDl (), SUPR (), and TMN (). These muscles’ levels of training were much higher than those of the others; thus, the results indicate that abduction as a rehabilitation movement was beneficial for the training of these muscles. However, the effects of an abduction in the coronal plane for the training of some muscles whose activations were lower than 0.22 were ineffective, including the SRA, INFR, and BICI, and of these, mostly the INFR was most noteworthy (). Furthermore, considering that a small external load is generally used in rehabilitation training, the SRA, INFR, and BICI may not be efficiently trained.

3.4. Influence of the Load Magnitude on the Effect of OLOC Training

For the aim of practical applications, an appropriate magnitude of external load must be determined before planning the rehabilitation. Figure 16 shows the effect of varying the magnitude of the external load of the mean activations of the main shoulder muscles involved in abduction with OLOC, and Figure 17 depicts the impact of the load magnitudes on the mean activation for some muscles more intuitively.

There were some muscles whose activations showed a low level of activation when using a small load but experienced a significant promotion with an increase in load size, such as the SRA, DLTs, TRP, BICs, and other similar muscles. For these muscles, a small increase in the load could lead to a significant promotion in muscle activation. Some muscles’ activations barely changed with the changing load magnitudes, such as the PMN, LTD, DLTc, and DLTa, among others. For these muscles, the training levels slowly promoted activity; thus, a change in the rehabilitation movement should be considered for these muscles. Therefore, when determining an external-load plan for training a specific muscle, an appropriate choice should be made between increasing the level of training and reducing the value of the load.

4. Discussion

The overall aim of this study was to propose a muscle-specific rehabilitation training method for the shoulder based on the optimal load orientation concept (OLOC). Biomedical engineering researchers developed human musculoskeletal movement models to reveal and imitate the structural characteristics and mechanisms that generate force in the musculoskeletal system. However, these results have not currently been fully utilized in the rehabilitation robots’ field. Similarly, training strategies and training trajectories are often based on the clinical experience of a doctor or therapist rather than making precise decisions based on existing research results from a musculoskeletal model, so the training effects cannot be guaranteed and can be difficult to quantify and evaluate. Likewise, most of the rehabilitation robots designed to rehabilitate human limbs are designed to train an entire malfunctional limb rather than training specific muscles in it. Different muscles have different situations related to their damage, stiffness, or degeneration, so uniform training cannot produce the best training effects. There have been some studies on specific-muscle training [1315], but the current research has great limitations because of the lack of use of the musculoskeletal movement mechanism. Therefore, a rehabilitation training method that is muscle-specific and based on the musculoskeletal model is proposed in this paper to provide a basis for the design of rehabilitation robots and further development of training strategies. Muscle activation was used in the musculoskeletal model to describe the force state on a specific muscle, and the optimal load orientation was proposed and calculated to maximize the training effect in some specific muscles.

The simulation results showed that training based on the OLOC could significantly improve the training effects in specific muscles more than simple weightlifting training. This method could be applied in rehabilitation robots designed to achieve a specific-muscle training function, which may significantly improve the practicality of using robots. This method also provided a method for quantifying the training effects on specific muscles during a given training process, which could be used to evaluate training effects and trajectory planning. For a specific muscle, an optimal rehabilitation movement may exist and can theoretically be designed to maximize the mean activation of the muscle.

Different muscles’ training effects were different in response to different load magnitudes. The training effect on some muscles increased significantly with an increased magnitude of the external load, while the training effect in other muscles was much smaller. Therefore, when designing a rehabilitation training program for a specific muscle, a reasonable choice should be made between raising the training effect and reducing the external-load magnitude according to its characteristics of how the training effect responds to the external-load magnitude.

Bifurcation may occur to the spatial distribution of some muscles, which may provide more possibilities for designing rehabilitation robots and developing training strategies.

This paper focused on the analysis of shoulder muscles, but the same methods could be used in other joints of the human body as long as there is enough anatomical data.

Data Availability

Most of the parameters of the shoulder model could be found in the website https://simtk.org/projects/dsem.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Project nos. 51475322, 51535008, 51775367, 51721003), the Tianjin Science and Technology Committee Program (Grant no. 17JCZDJC30300), and the Programme of Introducing Talents of Discipline to Universities (“111 Program”) under Grant no. B16034.