Abstract
Rehabilitation robotics has become a widely accepted method to deal with the training of people with motor dysfunction. In robotics medium training, shoulder repeated exercise training has been proven beneficial for improving motion ability of human limbs. An important and difficult paradigm for motor function rehabilitation training is the movement rhythm on the shoulder, which is not a single joint but complex and ingenious combination of bones, muscles, ligaments, and tendons. The most robots for rehabilitation were designed previously considering simplified biomechanical models only, which led to misalignment between robots and human shoulder. Current biomechanical models were merely developed for rehabilitation robotics design. This paper proposes a new hybrid spatial model based on joint geometry constraints to describe the movement of the shoulder skeletal system and establish the position analysis equation of the model by a homogeneous coordinate transformation matrix and vector method, which can be used to calculate the kinematics of human-robot integrated system. The shoulder rhythm, the most remarkable particularity in shoulder complex kinematics and important reference for shoulder training strategy using robotics, is described and analyzed via the proposed skeleton model by three independent variables in this paper. This method greatly simplifies the complexity of the shoulder movement description and provides an important reference for the training strategy making of upper limb rehabilitation via robotics.
1. Introduction
Rehabilitation robots have received increasing interest to provide rehabilitative therapy following neurological injuries such as stroke and spinal cord injury [1]. The shoulder is one of the most complex motor function areas of the human body, which has a direct impact on the recovery of upper limb motor function. When a shoulder has motor dysfunction caused by injury or disease, it often needs reasonable rehabilitation suitable for the specific conditions of the patient. For this situation, some studies have shown that a rehabilitation strategy based on task orientation that fuses high intensity repetitive exercise training will be very conducive to the improvement of athletic ability, which is one of the most effective rehabilitation methods [2–7], while the upper limb rehabilitation robot is a universally acknowledged effective way to carry out the rehabilitation strategy.
There are lots of robotics proposed currently for lower limb and upper limb rehabilitation. Most of them got fruitful achievement even though accurate motion information of bones and muscles was not considered thoroughly. There have been many attempts in robotics which evolved from a simplified 3-degrees of freedom (DoF) ball and socket, such as Carignan and Liszka [8] who have defined the shoulder as a 3-DoF ball-and-socket joint. However, it is not effective for a shoulder rehabilitation using exoskeleton robotics which has three orthogonal axes intersecting at one point to assume shoulder as a ball-and-socket joint, because the shoulder cannot be regarded as a single joint but complex and ingenious combination of bones, muscles, ligaments, and tendons. The remarkable characteristic of the shoulder complex is the movement rhythm involving the motion of the scapula, clavicle, and humerus. The rhythm analysis is important and necessary for shoulder complex rehabilitation training.
Biomedical researchers have done a lot of analysis on musculoskeletal models for shoulder training. In 1965, Dempster [9] established a simple two-dimensional series link model to describe the relative motion characteristics among bones of the shoulder; however, the model only has existing qualitative description, lacking the description with definition and quantitative analysis of the joint type. Engin and Chen [10] established a three-dimensional 6-DoF rigid-body model of the humerus relative to the human torso in 1986. Furthermore, Engin and Tumer increased the degrees of freedom of the shoulder complex model to 9 in 1989. Although Engin and Chen established a three-dimensional model of shoulder bones to analyze their kinematic characteristics, the series connection of the shoulder, without exception, simplified the scapula into two connecting rods that are connected to the clavicle and humerus.
There are complicated models in the shoulder joint. Garner and Pandy [11] established a 13-DoF human upper limb skeleton model in 1999. The model does not define the geometric constraints of the scapulothoracic articulation joint (ST) nor analyze the motion of the model. Based on the concept of this model, Maurel and Thalmann and Maurel et al. [12, 13] treat the ST joint as a single point contacting with the surface with 5 DoFs, but the study lacks systemic analysis. Tondu [14] proposed a similar shoulder mechanism model in 2005, which defines the ST joint as a planar subordinate on the basis of the Maurel model and ignores the rotational freedom of the sternoclavicular joint (SC) along the direction of the clavicle axis, modeling the shoulder girdle as a 2-DoF spatial parallel mechanism. This model simplifies the difficulty of analyzing shoulder girdles. In order to make the shoulder get a larger workspace and more activity, Lenarcic and Stanisic [15] in 2006 proposed a shoulder girdle skeleton system using a six-legged six-joint parallel mechanism to deal with the simulation, but it is difficult to determine the position of each joint and the length of the parallel bars corresponding to the actual shoulder girdle. It is also difficult to correspond to the movement of a single individual bone.
The current rehabilitation robots were designed with a simplified skeleton model and set the degrees of freedom such as the robotic-arm exoskeleton proposed by Klein et al. [16]. Therefore, this paper provides a suitable model that describes the shoulder complex motion via a model of a spatial hybrid mechanism with four rods and four joints and verifies the accuracy of the model through real data for rehabilitation robots in the shoulder joint. Among them, SC joint, acromioclavicular (AC) joint, and glenohumeral (GH) joint are all defined as spherical hinge joints, and the thorax is approximated as an ellipsoid. The constraint of the ST joint is defined as the two fixed points on the suprascapular connected to the thorax ellipsoid in a point contact constraint, and the rotation of the SC joint around the clavicular axis is considered to be an extra DoF.
2. Modeling and Analysis
In the model of the shoulder musculoskeletal dynamics, the posture analysis of the bones constituting the shoulder is the basis for shoulder motion analysis, muscle drive characteristics, muscle strength, and joint force analysis. Due to the complexity of the shoulder skeletal system, it is difficult to obtain the people’s posture of the shoulder bone quickly by using the commonly used real-time skeleton marking method. Based on results provided by Garner and Pandy [11, 17] for the anatomical results of human shoulder geometry, the human shoulder skeleton system is simplified as a hybrid spatial mechanism model to simulate and analyze the movement of the shoulder bone. The motions of each joint and bone are described by defining the local coordinate systems fixed on the corresponding bones, and the position inverse analysis of the mechanism is completed by the homogeneous coordinate transfer matrix and vector method, which simplified the analysis and prediction of the shoulder skeleton configuration. Finally, this chapter verifies the rationality of the model through the experiment of the actual human shoulder.
2.1. A Proposed Model of Shoulder Skeleton System
The kinematics of the shoulder of the human body mainly depends on the skeleton system composed of the shoulder girdle and humerus compared with the muscle system and the tendon system, so in this paper, the skeleton system for the shoulder is mainly discussed. It can be considered a closed-loop rigid system composed of the thorax, clavicle, and scapula. Joints involved in shoulder movement involve the AC joint, the SC joint, the ST joint, and the GH joint.
In this paper, a hybrid spatial mechanism model is used to simulate the shoulder skeletal system, as shown in Figure 1. In this model, the SC joint, the AC joint, and the GH joint are denoted as , , and , respectively, and are treated as the ideal 3-DoF ball-and-socket joint. The scapula achieves two-point contact with the thoracic ellipsoid surface through its medial lateral edge and , which may be equivalent to cylindrical-planar pairs with 4 DoFs. Among them, components 1, 2, 3, and 4 represent the clavicle, the scapula, the thoracic surface, and the humerus of the upper arm, respectively. The long axis of the humerus represents component 4, that is, the direction of connecting the center of the GH joint and the center of the elbow , where the elbow joint is defined as the midpoint of the outer lateral ankle EL and the inside ankle EM.
2.2. Kinematics Analysis of the Proposed Model
In the above mechanism model, components 1, 2, 3, and 4 constitute a hybrid spatial mechanism containing a closed chain and an open chain, where components 1, 2, and 3 form a closed-loop chain.
2.2.1. To Calculate DoFs of the Proposed Mechanism
To obtain the kinematics of the mechanism, we should first analyze its DoFs. For the closed loop of the shoulder girdle, the DoF of the mechanism () is obtained by the Kutzbach-Grubler formula [18]. where is the number of mechanism components, is the number of joints, and is the relative freedom of the th motion pair. The closed-loop part has three components and three motion pairs. Both SC and AC are ball-and-socket joints with 3 DoF. The ST joint is a cylinder-plane pair with 4 DoF. So a closed chain with 4 DoF in the shoulder joint system can be obtained. Considering the humerus movement, the total degree of freedom of the shoulder is up to 7. However, in the process of activity of the shoulder, the rotation of the clavicle around its axis is very small and it has an internal DoF that does not change the whole posture of link 4 [19–21]. Therefore, for joint , the rotation in the direction of the winding of component 2 is an extra freedom. If the rotation angle is known, the motion posture of the shoulder girdle can be obtained by inputting three joint variables.
2.2.2. To Establish Shoulder Global Coordinate and Local Coordinates Fixed to a Single Bone
In order to describe and analyze the kinematics of the skeleton model of the shoulder complex, every bone is fixed to one coordinate system, as shown in Figure 2. The origin point of the global coordinate system is the SC joint at point , and axes , , and are, respectively, parallel to three intersections of the human anatomical coronal plane, the sagittal plane, and the transverse plane. For the clavicle system , the origin is located on the SC joint at point . The axis is in the direction of the clavicle axis, and the axis is defined on the horizontal plane and perpendicular to . For the scapular system , its origin point is located at the point of the AC joint at , where the axis is in the direction from to the medial margin of the scapula and the direction of is perpendicular to the plane determined by points , , and . is parallel to . For the humeral coordinate system , the origin point is located at the center of the glenoid joint, where the axis is in the direction from to ; the axis is perpendicular to the plane determined by point , EL, and EM; and the axis is determined by the right-hand rule.
2.2.3. Kinematics Analysis of the Joints in the Shoulder
In the shoulder girdle, the thorax, the clavicle, and the scapula are interconnected with each other through the joints SC, AC, and ST, thus forming a closed-loop mechanism. For the closed-chain mechanism, position analysis can be attributed to determine the remaining joint variables and skeletal posture using the known mechanism geometric parameters and at least three joint variables as inputs. In this paper, the purpose of inverse analysis is to verify the accuracy of the actual shoulder movement predicted by the model.
In order to describe the matrix conveniently, the distance between joint and joint is set as , the distance between the center of joint and point is , the distance between and point is , and the length of the humerus is .
According to the above definition for the coordinate system, transition from the system to , namely, the motion of the thorax relative to the thorax, can be expressed via the following matrix: where and denote sine and cosine, respectively.
The conversion from system to system can be assumed to have moved along the -axis firstly and rotated around the -axis, then rotated around the x-axis, and finally rotated around the -axis. The conversion matrix is
The position and scale of the thoracic ellipsoid are known. Assuming that the position of the ellipsoid center point in the fixed system is , the transfer matrix of ellipsoid relative to fixed system is
If the semilong axes of the ellipsoid are, respectively, , , and and point and point are located in the thoracic ellipsoid table surface, the position coordinates can be represented by the parametric equation of the ellipsoid: where , , , and in the formula, respectively, describes the position parameters of point and point . Point and point are also located on the shoulder blade, as shown in Figure 3. The position coordinates of and in the scapula system are where in the formula is the included angle of and , because points , , and are all on the scapula, and , , and are constant. In summary, using the location transform relationship of and , the establishment of displacement constraint equation is as follows:
Formula (8) above can be expressed as follows:
Equation (9) is the joint position solution constraint equation of the girdle mechanism. Let each component be zero; the equations can be transformed into 6 constraint equations which contain a total of 10 joint variables. Clavicle rotation around its own axis is an extra degree of freedom, and if is zero, just any three rotational variables can be used as inputs to solve the remaining variables.
The kinematic analysis of the shoulder girdle is completed, while the humerus is connected to the shoulder complex through the GH joint as the end of execution. The motion posture of the upper arm is usually studied in the global reference system, using the ISB standard recommended by Wu and Cavanagh [22, 23] to describe the posture of the humerus. The rotation matrix applied to convert the global system into the humerus system is
If the position of the GH center coordinates is given, the conversion matrix for the translation from the humeral system to the global system is
2.3. Verification for Shoulder Skeletal Motion Based on Reference Data
In this paper, the size of shoulder bones refers to the anatomical data of human upper arm bones and joints obtained by Garner and Pandy [11, 17, 24]. Shoulder Database V1.1, published by Bolsterlee et al. [21] in 2014, is used as a reference source of shoulder motion data. Then, through the joint displacement constraint equations obtained above, the movement of the mechanism is analyzed by using Matlab.
2.3.1. Shoulder Mechanism Geometry
In this paper, the size of the mechanism is based on the anatomy of human bones and joints obtained by Garner and Pandy; the size of each component shown in Figure 1 is as shown in Table 1 below. The semimajor axis of the thorax ellipsoid and the center point of its ellipsoid are shown in Table 2. The coordinates of the center point are in the fixed system .
The position of GH joint center in the scapular is
2.3.2. Description of Shoulder Joint Movement
The local coordinate system of each skeleton in Figure 2 is defined to be calculated conveniently. The ISB standards recommended by Wu and Cavanagh and Wu et al. [22, 23] are widely used in biomedical engineering to describe the human upper body bones and joints in the movement. In the ISB standard, the local coordinate system of the shoulder bones and the bony landmark are shown in Figure 4.
In this paper, the shoulder skeletal configuration analysis results and experimental measurements which are unified according to the ISB standard will be compared; therefore, it is necessary to figure out the corresponding relationship between the calculation of the coordinate system and the ISB coordinate system. Klein Breteler et al. [25] obtained a more favorable shoulder morphological datasets by dissecting the body of a 57-year-old male and establishing a universal dataset of human shoulder bones. Table 3 shows the relevant parameters of bone markers and the size of the thorax of the general model and the experimental subjects in this article. All the coordinates are converted corresponding to the thorax system; the length unit is mm.
Through the definition of each local coordinate system and the characteristic sizes of each marked point of the shoulder skeleton in Tables 2 and 3, the relationship between the coordinate systems , , , and used in the above calculation and each coordinate system , , , and recommended in ISB standard can be described by the following rotation transformation matrix:
2.3.3. Definition of Shoulder Experimental Motion
Shoulder Database (Shoulder Database V1.1) published by Bolsterlee et al. in 2014 provides data of shoulder motion which has been used in the analysis [21, 26]. The experiment collected data from five healthy adult males. This article selects one of the objects (age 29, height 186 cm, and weight 109 kg) from the database. The experiment used a set of motion detection system (Optotrak System, Northern Digital Inc., Waterloo, Ontario, Canada) to detect the spatial positions of six marker groups fixed to the thorax, scapula, humerus, forearm, and palm shown in Figure 5. The locations of the bone markers are arranged according to ISB standard, and the series of the position of the markers are collected at a frequency of 100 Hz.
In order to get accurate shoulder movement data, during the experiment, the subject was asked to complete the range of motion (ROM) tasks including shoulder abduction (ABD), flexion (FLEX), and scapular plane movement (SCAP) which contain all the basic movement of the shoulder shown in Table 4. With processing, the data of the movement posture of each bone in the shoulder can be uniformly converted into the value in .
2.3.4. The Solution and Verification for Posture and Position of Each Bone
Through the processing of the spatial coordinates of each marker measured in the experiment, the joint variables are obtained, and the three independent joint displacements are used as input, and then the remaining six joint variables are obtained, as well as the shoulder girdle movement gesture, which is finally compared with the experimental measurement results.
Figure 6 shows the thorax movement gesture of the experimental subject in the ABD exercise test, and the posture of the humerus under YXY rotation sequence is described by three corresponding Euler angles HumY1, HumX, and HumY2. Figure 7 shows the measured posture of the clavicle relative to the thorax, and according to ISB standard, the posture of the clavicle is described under YXZ rotation sequence. The three Euler angles formed that depended on them are ClavY, ClavX, and ClavZ, respectively.
Thus, at any time, the posture matrix of the clavicle relative to the thorax can be obtained:
Suppose that
Substitute (17) into (12) and (13):
Let the third column of (18) and (2) be the same, and the following three constraint equations can be obtained: then
Figure 8 is the result of the input joint variables and .
Assuming that the object is at a certain moment of ABD movement, the position vector of the sternoclavicular joint marking point SC is , and the conversion relationship between and can be expressed as the following transfer matrix:
Thus, the position vector of the lower scapular AI in the ellipsoidal system can be obtained:
Combining formulas (4) and (12) and Table 3 on the thorax data, we can have the following formula:
In the ABD exercise of the experimental subject , the position of the lower scapular angle AI in the ellipsoid is shown in Figure 9.
However, the fact is that the AI point is located on the shoulder blade but not always on the surface of the ellipsoid, which means it does not exactly coincide with point in the model. In this paper, the projection point of the subscapular AI along the centerline of the ellipsoid is regarded as the equivalent point . As shown in Figure 10, the connection between AI and ellipsoid center intersects at on the ellipsoid surface, and the position and size of the chest ellipsoid of the experimental subject refer to Table 3.
It can be obtained that . Substitute (6) and (22), then there is
Another input joint variable can be found:
For the summary, for the object in the ABD motion experiment, the input joint variable for the position analysis of the mechanism is shown in Figure 11.
Then, entering the joint variables into closed solution equations, the joint variables and are obtained as shown in Figures 12 and 13.
Thus, the shoulder girdle mechanism position analysis is finished, and the results in the calculated coordinate system of the bone or joint movement posture can be obtained and then converted corresponding to the ISB standard system, so the agency model for shoulder skeletal motion prediction results can be obtained and compared with the experimental measurement results. Figure 14 shows the scapula relative to the thorax motion posture analysis results (ST joint). Figures 15, 16, and 17 are the results of the components of the scapula posture matrix after the decomposition of YXZ rotation sequence according to the ISB standard. Thanks to the fruitful achievement in biological research, the detail and accurate experimental data on shoulder motion can be provided, and in this paper, we selected a set of representative data to verify that the kinematic description of the skeleton model and experimental measurements are consistent.
Using the same method, the exercise experiments of FLEX and SCAP can be analyzed, and the agency model shoulder motion configuration prediction and experimental measurement results are basically consistent.
3. Shoulder Rhythm Model
From the view of the mechanistic model, the three rotational degrees of freedom of the humerus are independent from the three degrees of freedom of movement of the shoulder girdle mechanism. However, a great number of studies have shown that when the shoulders do some basic movements in the natural state (such as abduction, adduction, external rotation, pronation, flexion, and anteversion), there is a significant correlation and repeatability of shoulder girdle and humeral motion, following the law of the rhythm of the shoulder joint.
This article simplifies shoulder joint rhythm analysis to the problem of establishing the relationship between the humeral posture and the three independent input joint variables of the shoulder girdle mechanism.
3.1. Rhythm Analysis of Three Kinds of Exercise Task Experiment
The shoulder rhythm of an individual is highly repeatable and less affected by the size of the load when shoulders are in the condition of lighter load [27, 28], but there exist individual differences [29]. Therefore, this paper extends the traditional shoulder joint rhythm and describes the relationship between the scapula and the clavicle posture and the posture of the humerus individually via the proposed model with three independent joint variables . Through the constraints in the proposed model considering the individual size of shoulder complex, three independent joint variables can be further used to determine the corresponding relationship between the entire shoulder girdle posture and the humeral posture. This method reduces the number of variables involved in rhythms and increases the impact of individual differences.
Studies of shoulder rhythms by the literature [17, 25–28] showed that the elevation plane (HumY1) and the elevation (HumX) are relatively independent in shoulder bone movement, while the shoulder girdle posture with two joint variables has a strong regularity. When modeling the shoulder rhythms, this research also uses planes of elevation and elevation as independent variables for rhythm functions.
3.1.1. ABD Sports Rhythm
During the ABD exercise, five actions of lifting and dropping humerus are finished, as shown in Figure 6. Although the subject was asked to do abduction in the frontal plane, the actual elevation of the plane mainly focuses in the scoped of −10° to 20°, as shown in Figure 18. Although when the shoulder moved, the joint variable of elevation plane and the joint variable of elevation are considered independent from each other, the strong regularity and reflexivity are shown in Figure 18, because it represents the physical coordination of the experimental subjects in a certain motion like ABD exercise, and the coupling phenomenon will disappear in arbitrary motion. The relationship between the elevation plane and the elevation in Figure 8 is not universal, so this article will use the plane of elevation and the elevation as two independent variables of the law function.
Figure 19 shows the variation of various dependent variables with the humerus elevation angle when the subject does ABD exercise. is used to determine the rhythm of the shoulder girdle, while Figure 19(d) is used to determine the rhythm of humerus rotation around its shaft (HumY2, internal/external rotation). Cubic polynomial fitting is adopted, and shoulder rhythm function can be obtained: where denotes the polynomial fitting coefficients to the corresponding joint variables. For simplifying the calculation, the rhythm function corresponding to the ABD motion is expressed as follows:
(a) θ1
(b) θ2
(c) ϕd
(d) HumY2
Table 5 shows the fitting results of the dependent variables, and the accuracy and goodness of the fitting curve are, respectively, measured by root mean square error (RMSE) and coefficient of determination (R2).
3.1.2. SCAP Sports Rhythm
During the entire SCAP exercise, lifting and dropping actions in the humerus were completed 4 times. To analyze the ABD motor rhythm, a rhythm function is constructed.
Figure 20 shows the fitting result curve, and Table 6 shows the polynomial rhythm function coefficients, RSME, and coefficients of determination.
(a) θ1
(b) θ2
(c) ϕd
(d) HumY2
3.1.3. FLEX Sports Rhythm
During the entire FLEX exercise, lifting and dropping actions in the humerus were completed 4 times. The same as the ABD motor rhythm analysis, a rhythm function is constructed.
Figure 21 shows the fitting result curve, and Table 7 shows the polynomial rhythm function coefficient, RSME, and coefficient of determination.
(a) θ1
(b) θ2
(c) ϕd
(d) HumY2
3.2. Shoulder Rhythm Expansion
The lifting plane is not uniquely fixed during ABD, SCAP, and FLEX movements. In the ABD, SCAP, and FLEX, the humerus lifting plane angle HumY1 changes as shown in Figures 16 and 22. It can be seen that the main movement range of the plane of elevation during these three movements is as shown in Table 8.
(a)
(b)
Theoretically, the variables of the shoulder joint rhythm function are HumY1 and HumX in the whole range of activity, and the relationship between the movement of the shoulder bones and the axial rotation angle HumY2 of the humerus is not significant. At the same time, axial rotation angles HumY2 and HumX are related. Therefore, in this research, the shoulder joint rhythms corresponding to ABD, SCAP, and FLEX are extended to the motion space and , and the rhythm can be expressed as follows:
4. Conclusion
In this paper, a new shoulder complex skeletal model based on joint geometric constraints is proposed using a spatial hybrid mechanism to describe the movement of human shoulder complex, which is developed for the rehabilitation robots. The effectiveness and viability of the proposed model have been verified via the experimental data of shoulder skeletal motion. Based on this model, through the analysis of the angle and change rule of each joint in the shoulder complex during the movement, the shoulder rhythm was obtained by means of polynomial fitting with HumY1 and HumX as independent variables to make it suitable to guide the training strategy for shoulder rehabilitation via robotics.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Natural Science Foundation of China (Project nos. 51475322, 51535008, 51775367, and 51721003), Tianjin Municipal Science and Technology Department Program (Grant no. 17JCZDJC30300), and the Programme of Introducing Talents of Discipline to Universities (“111 Program”) under Grant no. B16034.