Determination of biological joint reaction forces from in-vivo experiments using a hybrid combination of biomechanical and mechanical engineering software

In biomechanical field, several studies used OpenSim software to compute the joint reaction forces from kinematics and ground reaction forces measurements. The bio-inspired joints design and their manufacturing need the usage of mechanical modeling and simulation software tools. This paper proposes a new hybrid methodology to determine biological joint reaction forces from in vivo measurements using both biomechanical and mechanical engineering softwares. The methodology has been applied to the horse forelimb joints. The computed joint reaction forces results would be compared to the results obtained with OpenSim in a previous study. This new hybrid model used a combination of measurements (bone geometry, kinematics, ground reaction forces...) and alsoOpenSim results (muscular and ligament forces). The comparison between the two models showed values with an average difference of 8% at trotting and 16% at jumping. These differences can be associated with the differences between the modelling strategies. Despite these differences, the mechanical modeling method allows the computation of advanced simulations to handle contact conditions in joints. In future, the proposed mechanical engineering methodology could open the door to define a biological digital twin of a quadruped limb including the real geometry modelling of the joint.


Introduction
In order to take inspiration from biological articulations to design bio-inspired mechanical joints usable in complex mechanism, the joint reactions forces must be known. In medicine, some researches have used instrumented prostheses to measure forces and pressures in the replaced bones [1]. Such methods would later permit to improve the upcoming prostheses comfort and lifespan. Another way to compute joint reaction forces is the use of musculoskeletal numerical models of human [2][3][4][5] or animal [6][7][8][9] limbs to understand biological joints characteristics. These models are anatomically based 3D models in which fluids and most soft tissues are often neglected.
However, in the scope of imitating highly performing biological articulations to change the designs of joints in mechanical systems, these models are essential in order to mechanically characterize the biological joints. In Becker et al. [10], the authors presented a model of the horse forelimb computing its joint reaction forces with OpenSim.
Moreover, most studies working on biological joint reaction forces estimation used OpenSim software for their modelling and simulation. Nevertheless, for the design and manufacturing of bio-inspired joints, the tools of a mechanical modeling and simulation software are necessary.
In this work, a hybrid method was proposed using both biomechanical and mechanical engineering software to compute the joint reaction forces from in vivo measurements. These joint reaction forces computations results were then compared with OpenSim results. The data used in this study are available at https://simtk.org/docman/? group_id=1728.

Specimen and measurements
A French Saddle sport horse weighing 560 kg and judged free of obvious lameness was ridden by a professional rider and their total mass with the equipment was estimated to be 650 kg. The kinematics of the horse right forelimb and the ground reaction force were simultaneously recorded for trotting and jumping a fence of 1 m high. The recordings were made for ten trials each time. For the kinematics measurements, markers were placed along the horse forelimb, principal markers at the centers of joints and secondary markers between principal ones. A board camera followed the horse parallel to its movement at trotting and jumping. The kinematics were therefore only measured in the sagittal plane. This simplification could have been done because the horse belongs to these animals, which have limbs that mostly move in a vertical plane, and therefore the hypothesis was made that movements and forces were very small in lateral direction. A video analyzer called Kinovea (OpenSource Software, Kinovea Association, France) was used to follow the markers trajectories. The reference frame was chosen with a vertical Y-axis, a longitudinal Z-axis and a transversal X-axis. Due to this non-invasive measurement system, the trajectories needed to be corrected before they were used as entry kinematics data. One of the crucial compensations was the skin artifact that was handled with a method inspired from the solidification method [11]. The entire correction process was developed and detailed in [10]. As the mechanical engineering software does not use optimization process to best-fit the kinematics, the movement of the limb was controlled with the angular coordinates. From the markers coordinates resulting from the corrections cited before, the angular variations were mathematically recalculated at each joint to define the entry kinematics data for the mechanical software simulations. These angular variations were validated by comparison with earlier studies in the discussion part of this paper. During the kinematics measurements, the ground reaction forces in the three directions were measured with a force plate buried under 100 mm of sand. This gave ground reaction forces values and the center of pressure at each step of measurement could be defined. The ground reaction forces values across the stride were published earlier [10] and the values were validated by comparison with literature.
In literature, authors often used CT (computed tomography) scans [8], or MRI (Magnetic Resonance Imaging) on cadavers [9] to obtain the bones geometry. Here, the forelimb bones of a horse skeleton were scanned with an industrial optical 3D digitizer, the GOM ATOS 3 (GOM, Braunschweig), that enabled to access to the geometry with a Computer aided design (CAD) software like CATIA V5 (Dassault Systems, Vélizy-Villacoublay). The numerical data of bone geometry obtained from the GOM ATOS 3 were uploaded on CATIA V5 software in order to fit to the dimensions of the horse used for kinematics experiments. The distances between the principal markers are supposed to correspond to the bones lengths. By dividing the distances between the principal markers at rest by the measured bones lengths, a scale factor of 1.2 ± 0.03 was deduced. This scale factor was applied to the bones CAD in all directions. On CATIA V5 software, the surfaces were filled and assigned with a homogenous material with the same mean density as bone (r bones =1800 kg.m À3 ). All experimental methods used in this study were detailed in [10] .

Mechanical multibody modelling
Simcenter mechanical engineering software (Siemens, Munich), combines system simulation, 3D Computer Aided Engineering (CAE) and is widely used in industry to simulate mechanical systems and predict their behavior. The joint reaction forces are the combination of kinematics, ground reaction forces and muscular and ligament forces, so all these elements would be needed to compute the joint reaction forces. First, the limb geometry measured earlier with the GOM ATOS 3 was imported to the mechanical engineering software. A six segments (humerus, radius, wrist, metacarpal and carpal bones, first phalanx, second phalanx and foot) model was defined, leading to five articular joints: shoulder, elbow, wrist, metacarpophalangeal (MCP) joint, inter-phalangeal joint (P2-P1). The kinematics frames defining the axis of rotation of each joint were defined with the same methodology than the one used in the previous study [10]. These five joints were modeled with revolute joints leading to a five degree-of-freedom model. Even if the kinematics was recorded in 2D, the revolution axes of each joint are not perpendicular to the sagittal plane and not aligned, so the simulated limb would not exactly move in a plane. To determine the joints axes of rotation, the surfaces of each articulation were best fitted with simple surfaces with CATIA V5 software (spheres, cylinders) to define the characteristic axes and centers of the joints.
In the mechanical engineering software, it was first necessary to define a FEM super-element (containing masses, inertial properties) of the model called a skeleton [12], which enabled to define the joints kinematics. Two nodes were therefore defined at the proximal and distal center of joint of each bone, and one third node was defined at the middle of these two. The FEM super-element of the model is represented in Figure 1. On this representation, only the proximal and distal nodes of each bone are represented. These three nodes were linked by rigid body element connections because bones are considered as rigid. Hinge joints were defined between the two coincident nodes at each joint with the corresponding frame (for example, a hinge joint was defined between Scapula [3] and Humerus [1]. The bones geometries were meshed with 3D tetrahedral elements and the resulting mesh nodes were connected to the middle node of each segment with rigid body elements. The 3D meshing would be necessary for the insertion of muscles. To control the kinematics of the limb, a constraint of displacement containing the top marker corrected trajectory was assigned to the proximal node of the scapula (node named Horse [1] on Fig. 1) in order to control the overall displacement of the limb. The angular displacements previously computed were used to command the hinges rotations at each joint.
Concerning the ground reaction forces, the values resulting from the force plate measurements were directly applied to the distal node of the foot segment. Indeed, at each instant of measurement, the moment of the force was recalculated from the center of pressure at the time to the coordinates of the distal node of the foot segment. The corresponding kinematics and ground reaction forces are given in Figure 2 for trotting and in Figure 3 for jumping.
The modelling of muscles and ligaments was a crucial part of the estimation of joint reaction forces. In the previous work, the muscles were first all modeled in 3D on CATIA V5 by using a software of anatomical observation (Biosphera). These muscular elements were then assembled on the bones geometry on the CAD software. This 3D modelling of the muscles on CATIA V5 enabled to deduce the precise coordinates of muscles insertion points and also the inertial properties of the limb. In the previous model built with OpenSim [10], NMS builder was used to build a model actuated with the 23 muscle-tendons and 5 ligamentous passive elastic structures. They were modeled with their insertion points, their inertial properties, their muscular parameters (maximal force in the muscle, fiber length, and angle of penation) and wrapping surfaces to prevent them to cross bones. The Static Optimization tool of OpenSim enabled to obtain the muscular activations and forces over time during the complete simulation. More precisely, the Static Optimization tool computed the muscular forces necessary to reach the kinematics inputs and the ground reaction forces applied to the limb. Some of these muscular forces are represented for the trot in Figure 4b and for the jump in Figure 4c. In the muscular forces computation with OpenSim, it was necessary to reduce the reserves (virtual torques used by the software if  the model cannot reach the entry data) to have joint reaction forces values as realistic as possible. The methodology to reduce these reserves was detailed in Becker et al. [10], and it consisted particularly in blocking the MCP joint. There are very few muscles at the MCP level, the muscles are mostly localized in the proximal part of the limb. The long muscle-tendon structures going down to the distal part, are attached at different insertion points and are guided by wrapping surfaces. This allowed to assume that the blocking of MCP would more impact the damping of the limb than the muscular forces values. These muscular forces could therefore be used for the Simcenter modelling, where the MCP is moving. In the current model using Simcenter, the muscles and ligaments were represented by FEM super-elements containing their inertial properties. The OpenSim resulting muscular forces were directly applied at the insertion points of muscles located on the limb bones. In practice, forces were applied directly on the mesh nodes corresponding to the insertion points of muscles (Fig. 4a).

Kinematics results
It is important to remind that the 2D measurements of the kinematics lead to a 2D model even if some data are in 3D. The mean trotting speed was 4 m/s and the mean speed of the body of the horse at jumping was about 5 m/s. As explained in materials and methods section, the trajectories of each markers were corrected and then the angular openings of each joint were automatically recalculated. For the shoulder, the rest angle was 125°and this angle varied from À4.5 to +14% at trotting. The initial position of the elbow was 150°and across trotting it went from À39 to +5.8%. The wrist was initially opened of 175°and it closed up to À38% and opened up to +9.2%. The MCP joint has an initial opening of 154°, and closed up to À1.3% but opened up to +56.2%. Finally, the P2-P1 angle joint had an initial position of 172°and at trotting it closes only up to 0.7% but it opens up to 29%.
At jumping, the angular variations were globally wider than for trotting. At the shoulder joint, the angle varies from À12 to +8.2%, at the elbow joint, it varies from À52 to +7.2% and at wrist joint it goes from À58 to +9.4%. The MCP joint closes up to À13.4% and opens up to +49.8%, and finally the P2-P1 angle varies around its rest position from À0.5 to +29.1%.

Joint reaction force comparison at trotting
For the comparison of joint reaction forces results between the two models, the norm in the YZ place was considered. The comparison of joint reaction forces at the five forelimb articulations for the trotting stride obtained with both models is given in Figure 5. In this figure, joint reaction force obtained with OpenSim is represented with a red dashed line while the Simcenter modelling joint reaction force is represented with a green solid line.
It is important to keep in mind that the joint reaction forces are a combination of kinematics, muscular forces and ground reactions forces. This is confirmed by comparing the ground reaction forces (Fig. 2) and the muscular forces ( Fig. 4b) with the joint reaction forces. It appears that the joint reactions forces at swing phase are rather low (when ground reaction forces and muscular forces are low) and that they really increase at stance phase (when ground reaction forces and muscular forces are high).
The percentages of difference for the extremum values at trotting are given in Table 1.
The joint reaction forces computed with the multibody model were quite similar to those computed with the OpenSim model. For the wrist joint, the MCP joint and the P2P1 joint, the results present at most 2.6% of difference. For the proximal joints, the results are quite more different with 17% of difference at the shoulder joint and 18% of difference at the elbow joint. This results in a mean difference of 8% between the two models.

Joint reaction force comparison at jumping
The jumping joint reaction forces norm comparison between the two models at the five forelimb articulations is given in Figure 6. In this figure, joint reaction force obtained with OpenSim is represented with a red dashed line while the Simcenter modelling joint reaction force is represented with a green solid line. Here, one can notice again the synchronization between the ground reaction forces (Fig. 3), the muscular forces (Fig. 4c) and the joint reactions forces (Fig. 6). The percentages of difference between the two modelling methods are given in Table 2.
At jumping, the joint reaction forces computed with the multibody model are broadly lower than the results from OpenSim model. These differences are mainly noticed for the three proximal joints (shoulder, elbow and wrist). Indeed, the joint reactions forces computed with Simcenter

Validation of kinematics results
The recalculations of joint angles are very close to the previous studies at trotting [10,13] and at jumping [10]. In Becker et al. [10], the correction methods applied to the kinematics were the same but the angular variations were obtained from OpenSim. These comparisons would therefore validate the correction methods applied to the kinematics data because of the resemblance with Dutto et al. [13], but also to validate the inverse kinematics tool of OpenSim which earlier enabled to obtain almost the same values as here. Peak shoulder and elbow joint angles for trotting were 142.5°and 158.7°respectively, which were very similar to the values reported by Dutto et al. [13] (138°± 5°and 150°± 10°) and Becker et al. [10] (144.5°and 158.7°). For the wrist and MCP joints, the peak recalculated values are 192.2°and 240.6°respectively and those are also very close to the values given by Dutto et al. [13] (186°± 9°and 231°± 4°), by Harrison et al. [2] (182°± 2°and 241°± 4°) and also by Becker et al. [10] (190.9°and 238.5°). Finally, for the phalangeal joint, the peak recalculated value was 221.9°which shows high resemblance with Dutto et al. [13] who reported a value of 220°± 3°and with Becker et al. [10] which gave a value of 220.5°.
At jumping, the joint angular variations were larger than at trotting, due to wider movements to reach the jumping kinematics. The recalculated values were very close to those obtained with the inverse kinematics tool of Opensim [10].

Discussion about joint reaction force results
As explained earlier, the joint reaction forces are a combination of kinematics, ground reaction forces and muscular forces. It can be observed that the joint reactions forces norm values obtained with both methods (Figs. 5 and 6) are linked with the evolution of the ground reactions forces (Figs. 2 and 3) and also with muscular forces (Fig. 4). During the swing phase, there is still some co-contraction of muscles to keep it upward resulting in some joint reaction force. This explains why there is still some low reaction force in joints during the swing phase.

Comparison between Simcenter mechanical engineering software and OpenSim joint reaction force results
There were some observable differences of joint reaction force values between the two models. In OpenSim, when the kinematics and the muscles parameters did not enable the model to reach the entry data, it used virtual torques at the center of joints named the reserves. These reserves should be as low as possible. In the OpenSim model [10], the kinematics used for the estimation of joint reaction force were simplified to reduce these reserves: the MCP joint rotation was blocked. Here, in the multibody model, the kinematics was not simplified for the joint reaction force computation: the five hinge joints rotations were allowed. Moreover, it seems that the MCP joint plays an important role in the landing phase. Indeed, when looking at this joint in the jumping stance phase (Fig. 7), it appears that it is the damper for the landing phase and the propeller for the next stride. Over the complete strides, its angle has an opening interval of 85°for trotting and 94°for jumping [10]. Therefore, it means that in the OpenSim model, the MCP joint cannot complete its damper-propeller mission and the other joints need to compensate this lack by handling more joint reaction force than needed.
This would cause differences in the observed joint reaction results between the two models. It can be assumed that taking into account a kinematics closer to reality leads to a more realistic estimation of the joint reaction force. To validate this hypothesis about MCP joint, the mechanical engineering software model was simulated with blocking the rotation of this joint. This led to joint reaction forces closer to the first model as seen in Table 3.
With the blocking of the MCP joint, the difference is reduced from a mean value of 8% to a mean value of 3.7% at trotting and from a mean value of 16% to a mean value of 8.5% at jumping. The influence of the blocking of MCP joint rotation on the joint reaction forces is therefore demonstrated. Nevertheless, only some in vivo measurements could be used to verify the joint internal loads without all modelling hypotheses, but this implies surgical procedures.

Limitations of the study
Unfortunately, this study is suffering from the difficulties of multi-scale modelling [14][15][16]. Indeed, biomechanical systems are difficult to model and simulate because it requires the establishment of many hypotheses regarding the biological system (kinematic uncertainties associated with markers position, muscle parameters, difficulty to define the rotation axes of biological joints...). In our study, the limitations of the mechanical modelling of a biological system were handled by simplifications and hypothesis that modified the simulation results somehow.
The first limitation is that our model mostly uses 3D data (geometry of bones, inertial parameters, ground reaction forces, muscle parameters…), but the kinematics measurements were planar. The trajectories recorded were only the projections of the markers on the camera image. This consequently led to 2D computation results. Nevertheless, 2D results are acceptable because the horse movement is mostly in a sagittal plane and that lateral movements and therefore forces are very low. In consequence, only the norms in the sagittal plane (Y axis and Z axis) were studied and compared.
It could also be considered as a limitation to use muscular forces resulting from the OpenSim model to compute joint reaction forces in the mechanical model. This limitation first comes from the difference in the entry kinematics used in both models. First, OpenSim optimizes the kinematics to recalculate angles of rotations. Nevertheless, when considering all the processing and corrections made on the kinematics data before importing them in the OpenSim software (corrections of speed, corrections of hoof marker trajectory, skin artifact correction), the compensation of kinematics resulting from OpenSim optimization process led to very low modification on the angles of rotations. This could explain partly the difference between the magnitudes of loads obtained in both procedures used in this study (Open Sim and Mechanical software). Concerning the MCP joint, the kinematics of the joint has been treated differently between the two models and this could concur to a variation in the joint reaction force calculation.
Finally, the torque reserves used in OpenSim model to accommodate the associated uncertainties of the experimental data are needed to have perfectly mechanically equilibrated systems. This contributes to the variations in the joint reaction values between OpenSim and Simcenter models. Nevertheless, this hybrid method constitutes a new strategy for the mechanical modelling of a horse forelimb and a solid basis for future works.

Conclusion
The overall goal of this work was to build a multibody model of a horse forelimb with mechanical engineering software to compute not only the joint reactions forces, but also to be able to perform usual stress and contact pressure analysis when required. A model developed with OpenSim software in a previous work [10] was used to define the applied muscular and ligament forces and compare the joint reaction forces. This new hybrid model used a combination of measurements (bone geometry, kinematics, ground reaction forces…) and also OpenSim results (muscular and ligament forces). For the muscular and ligament forces, only biomechanical software could manage muscles modelling and the OpenSim Static Optimization tool enabled to compute the muscular forces in the horse forelimb over time. These results were directly introduced in the mechanical engineering software simulations by applying the forces at the corresponding muscles insertions points on bones.
This hybrid model led to joint reaction forces values that were compared to OpenSim published values. Some differences were observed between the two models. These differences were mostly explained by the difference in the kinematics: in the OpenSim simulation the MCP joint was blocked while it was not in the mechanical engineering software simulations. In reality, this joint seems to have a damper-propeller mission and should therefore be considered. Other limitations are described in a dedicated Section 4.5.
The comparison between the two models showed values with an average difference of 8% at trotting and 16% at jumping. These differences can be associated with the differences between the modelling strategies. Despite these differences, the mechanical modeling method allows the computation of advanced simulations to handle contact conditions in joints. In future, the proposed mechanical engineering methodology could open the door to define a biological digital twin of a quadruped limb including the real geometry modelling of the joint.