Common modelling assumptions affect the joint moments measured during passive joint mobilizations

Joint resistance to passive mobilization has already been estimated in-vivo in several studies by measuring the applied forces and moments while manipulating the joint. Nevertheless, in most of the studies, simplified modelling approaches are used to calculate this joint resistance. The impact of these simplifications is still unknown. We propose a protocol that enables a reference 3D inverse dynamics approach to be implemented and compared to common simplified approaches. Eight typically developed children and eight children with cerebral palsy were recruited and underwent a passive testing protocol, while applied forces and moments were measured through a 3D handheld dynamometer, simultaneously to its 3D kinematics and the 3D kinematics of the different segments. Then, passive joint resistance was estimated using the reference 3D inverse dynamics approach and according to 5 simplified approaches found in the literature, i.e. ignoring either the dynamometer kinematics, the measured moments alone or together with the measured tangential forces, the gravity and the inertia of the different segments, or the distal segments kinematics. These simplifications lead to non-negligible differences with respect to the reference 3D inverse dynamics, from 3 to 32% for the ankle, 4 to 34% for the knee and 1 to 58% for the hip depending of the different simplifications. Finally, we recommend a complete 3D kinematics and dynamics modelling to estimate the joint resistance to passive mobilization.


Experimental procedure
A passive testing protocol was designed to obtain continuous joint angle and joint resistance measurements while the subject's joint was manipulated using a 3D handheld dynamometer (Sensix, Poitiers, FR, 2000 Hz), measuring triaxial applied forces and moments, through full available sagittal range of motion (ROM) (Fig. 1).Three joints were mobilized three times in five different positions (P1 to P5) at low (LV) and high (HV) velocity.It is important to note that only the stretch of the antagonist muscles is made at high velocity, the return to the starting position is made slowly: • Ankle with the knee at 90° and 0° (positions P1, P2) • Knee with the hip at 90° and 0° (positions P3, P4) • Hip with free knee only at LV (position P5) During all stretches, 3D kinematics of the pelvis and lower limbs was determined using a marker set, placed over specific body landmarks (lower limb Plug-in-Gait marker set with two additional markers placed over the iliac bone to take into account that the posterior superior iliac markers are not visible when the patient is lying) and a motion analysis system (15 cameras VICON, Oxford, UK, 100 Hz).Moreover, six markers, rigidly attached to the dynamometer, were used to track its position and orientation.Dynamic data (dynamometer forces and moments) were synchronously recorded and then downsampled at motion frequency.www.nature.com/scientificreports/EMG activities were synchronously recorded using pre-amplified dual differential surface electrodes (DE-2.1,DelSys, Inc., Boston, USA, 2000 Hz) placed over the rectus femoris, vastus lateralis, semitendinosus, tibialis anterior, peroneus longus, soleus and gastrocnemius muscles.Electrode locations, determined by the same investigator for each subject, according to the Surface Electromyography for the Non-Invasive Assessment of Muscles (SENIAM) guidelines, were prepared by shaving the skin and cleaning with alcohol.EMG activities were measured to check muscle inactivity during low velocity stretches or the occurrence of muscle reflex during high velocity stretches.

Kinematics
A seven segment, 24 degrees-of-freedom model of the pelvis and lower extremities was used to characterize joint kinematics.Six degrees of freedom were used to define the position and orientation of the pelvis, three degrees of freedom were used to define each hip, knee or ankle joint.
Body segment coordinate systems, joint centers and segment lengths were established using the marker positions collected during a calibration trial.For each segment i , a reference frame and the homogeneous matrix, T 0,i , defining both the orientation ( R 0,i ) and the position (P 0,i ) of the segment reference frame with respect to the global reference frame, was determined.Then, joint angles were calculated using the Cardan angles.
In the same way, using dynamometer marker positions, a dynamometer reference frame and its homogeneous matrix,T 0,Dyn , was determined.

Inverse dynamics
Inverse dynamic analysis, using homogeneous matrix, was used to compute the hip, knee, and ankle joint moments 17,18 .A free body diagram example for the position P3 is proposed in Fig. 2. The convention used to represent the curves of moments corresponds to that recommended by the International Society of Biomechanics 19 .The reader can invert the ankle and hip curves to find a convention more generally used by the clinicians.
First, measured forces and moments from the handheld dynamometer, expressed in its reference frame, form the matrix: where F x , F y , F z are the components of forces and M x ,M y ,M z are the components of moments.
Then, this tensor is expressed in the global reference frame thanks to the dynamometer transformation matrix: where superscript t stands for matrix transpose.
The gravity acceleration matrix is given, in the global reference frame, by: (1)  9) to the shank segment.A Shank represents the inertial action matrix and is computed from the subject's shank inertial parameters and kinematics.ϕ g/Shank is the weight action matrix.ϕ Dyn/Shank represents the measured forces and moments applied to the shank by the handheld dynamometer accordingly to their kinematics.ϕ Foot/Shank represents the forces of the foot on the shank segment computed with Eq. (10).
Finally, subject's body segment inertial parameters were estimated using the regression equations developed for children by Jensen 20 .Segment pseudoinertial matrices were determined thanks to the parallel axis theorem and expressed in the segment reference frame: where [K O ] , q and m are, respectively, the segment Poinsot inertia matrix, first moment of inertia (product between the mass and the position of center of mass with respect to the reference frame origin) and mass.
Thus, through an iterative process, as described in the following algorithm, the joint moments (joint resistance in the present case), ϕ i/i+1 Ri , are obtained, at every instant, in their local frame.

If i is the mobilized segment
Else End End This inverse dynamic process with a complete 3D modelling of the joints and the dynamometer was considered as our reference.According to the literature, several modelling assumptions and simplifications were studied.
Dynamometer kinematics.We considered that the dynamometer kinematics was the same as the segment on which it was attached (Case A).Then, T 0,Dyn was modified.Its orientation part was supposed the same as the solicited segment, i, and its position was assumed fixed in the segment frame.
with P i,Dyn the averaged position of the dynamometer in the i-th segment reference frame.
Thus, this simplification will mainly affect the expression of the dynamometer forces and moments in the global reference frame.
Measured forces and moments.We ignored some components of the measured forces and moments 9,11,14,15 .First, the measured moments were considered as null (Case B).
Second, only the main force of solicitation, as if we had a monoaxial load cell, was considered (Case C).
Distal segments movement.For the knee and hip tests, we studied the effect of considering that the distal segments (foot or foot and knee) are fixed relatively to the tested segment.To do so, distal segments kinematics was supposed the same as the tested joint (Case E).Thus, for the knee joint tests (positions P3 and P4): with P knee/ankle is the averaged position of the ankle in the knee reference frame and for the hip joint test (posi- tion P5), we have: with P hip/ankle and P hip/knee are the averaged positions of, respectively, the ankle and the knee joints in the hip reference frame.

Statistics
We evaluated the different simplifications by calculating the difference between the joint moments computed through the complete 3D inverse dynamics approach and the other methods at minimal and maximal angles (M θMin and M θMax respectively).With non-normal data (verified by Shapiro-Wilk test), M θMin and M θMax of the reference inverse dynamics method were compared at those of the different cases, pooling TD and CP children and pooling all velocities (except for Case D), with non-parametric Wilcoxon signed-rank tests for paired samples.The level of significance for the tests was p = 0.05/5 = 0.01 considering Bonferroni correction for multiple comparisons.Confidence intervals were calculated according to Campbell and Gardner 21 .
In order to estimate the amount of approximation made with the different cases over the entire tested ROM, root-mean-square differences (RMSD) were calculated in respect to the reference inverse dynamics.RMSD was also expressed as a percentage of the absolute maximum joint moment value (%Max).Correlation coefficients between reference inverse dynamics and the different simplification cases were also calculated.Inverse dynamics (Eqs.1-19) and statistics calculation were implemented in Matlab 2020a (MathWorks, Natick, USA).

Results
Differences in M θMin and M θMax and RMSD for all cases and positions are presented in Figs. 3 and 4. Figure 4 represents the relationships between the joint angle and the joint moment.Before pooling the results, although most CP patients demonstrated muscle reflex at HV, no significant differences, in the amount of joint moment approximation caused by the simplifications, according to pathology or velocity were found.Conversely, significant differences between the reference inverse dynamics and the different cases (pooling all results) were found.Results are presented below for each case considered.Detailed results are available in supplementary material.
RMSD values lead to a median approximation up to 7.6% of the absolute maximum joint resistance value measured.Consistency of the method relative to reference inverse dynamics was verified with high coefficient of correlation values (R > 0.98).( 14)   When we consider only a monoaxial load (Case C), significant differences with respect to the reference inverse dynamics method were found on M θMax for positions P1 and P2 (median differences When we consider only a monoaxial load (Case C), RMSD values increase and remain important for P1 and P2 and become non-negligible for P3 and P4 (%Max > 5%).Consistency of the method relatively to reference inverse dynamics remains verified with high coefficient of correlation values (R > 0.92).However, for some trials in positions P1 and P2, correlation was lower (0.4 < R < 0.9).
For both cases B and C, M θMax at the ankle joint is likely to be overestimated and M θMin slightly underestimated at positions P1 and P2 (Figs. 3, 4).

Inertia and gravity (Case D)
Significant differences with respect to the reference inverse dynamics method were found on M θMax for positions P1, P3, P4 and P5 (median differences [IQR] of 0. M θMax and M θMin at the knee joint are likely to be underestimated at position P3 (Figs. 3, 4).M θMax is likely to be also underestimated at position 4 whereas results are more variable concerning M θMin .M θMax and M θMin at the hip joint are likely to be underestimated at position P5.Ignoring inertia and gravity leads to important RMSD for positions P3, P4 and P5 with median approximations between 30 and 57% of the absolute maximum joint resistance value measured.For this case, the impact of this modelling assumption seems affected by the movement velocity.At low velocity, consistency of the method relatively to reference inverse dynamics remains verified (R > 0.98) but at high velocity the correlation decreases for positions P3 and P4 (R of 0.76 [0.54-0.89]and 0.83 [0.64-0.87],respectively).
In this case, low RMSD values were found for positions P3 and P4 with mean %Max inferior to 2.0%.For position P5, RMSD values were more important and lead to a mean approximation up to 6% of the maximum joint resistance value measured.Consistency of the method relative to reference inverse dynamics was verified with high coefficient of correlation values (R > 0.99).

Discussion
This study aimed to evaluate how different modelling assumptions that can be found in the literature affect the computation of the joint resistance to mobilization.To do so, we proposed a protocol that enabled a reference inverse dynamics approach to be implemented.Applied forces and moments were measured through a 3D handheld dynamometer, simultaneously to its 3D kinematics and the 3D kinematics of the different segments.
Significant differences have been found between the reference inverse dynamics and the different simplifications studied for the joint moment at minimal and maximal angles (M θMin and M θMax ).However, the discussion will be focused on simplifications that led to the approximations that would impact the evaluation of joint resistance.

Oversimplified modelling assumptions
Measuring the dynamometer kinematics is rarely done in the literature.Instead, its kinematics is supposed to be the same as the segment on which it is attached (Case A), assuming that the axes of the dynamometer are aligned with those of the manipulated segment 8,9,12,14 .However, this simplification leads to approximations on the orientation and position of the measured forces and moments, which cause differences with respect to the reference inverse dynamics method.Particularly, it was the case during the knee stretches (positions P3 and P4), where RMSD values or median differences on M θMax and M θMin are the higher.Indeed, in these test conditions, applying a movement profile to the dynamometer that enables to maintain a constant position and orientation of the dynamometer in the shank reference frame while mobilizing the knee joint seems more complex.Simplifying the dynamometer kinematics could lead to a slight underestimation or overestimation of the maximal knee extension moment at positions P3 and P4 respectively, and a slight overestimation of the maximal knee flexion moment at position P4 only [up to a median difference of 6% of the reference absolute maximal joint moment in position P3].Some authors do not consider measured moments or use monoaxial load cells 9,11,14,15 , whereas mobilizations are not purely in the sagittal plane (Supplementary Material-Fig.S1).Interestingly, ignoring the measured moments alone or also the tangential measured forces components leads to non-negligible differences for the ankle stretches (positions P1 and P2).For other positions, results seems to be more variable with some trials showing important differences with the reference inverse dynamics method.These results show that when applying a stretch, the evaluator does not apply a purely monoaxial force.This is especially true for ankle stretches, where both the measured forces and moments have a great impact.Thus, Case B and Case C enhance the importance of using a 3D dynamometer.Not measuring the full tensor underestimate the maximum plantarflexion and dorsiflexion moments [up to a median difference of 63% of the reference absolute maximal joint moment in position P2 for Case C].
Ignoring the inertia and the gravity of the segments (Case D) is a well-known effect in the literature [22][23][24] .Thus, in the literature, inertia and gravity are more frequently ignored for small segments (i.e. the foot) than for larger segments.This study confirms that this simplification is acceptable for the ankle stretches, where small differences with the reference are found, but not for the knee and hip stretches.This is in agreement with what is done in the literature, where this simplification is made only for the ankle stretches 8 .Moreover, it is interesting to see the effect of the velocity on the resulting approximation for the knee stretches (positions P3 and P4).Indeed, differences on M θMax and M θMin were greater at high velocity, and similarly, consistency with the reference inverse dynamics method, was impacted, with mean R values decreasing importantly.The inertia and the gravity of the segments must be considered when computing knee or hip resistance to mobilizations [up to a median difference of 55% of the reference absolute maximal joint moment in position P5].
The last assumption tested consists of ignoring the kinematics of the distal segments (Case E).Authors assume that during the mobilization, distal segments are fixed with respect to the manipulated one, neglecting the inertial and gravitational efforts due to the distal segments, as they have their own kinematics.This simplification appears to be acceptable for the knee stretches (positions P3 and P4), which is in agreement with Bar-On et al. 8 but starts to be non-negligible for hip stretches (P5), as shown by RMSD and M θMax and M θMin values.Measuring distal segment kinematics is recommended especially for hip assessments [up to a median difference of 7% of the reference absolute maximal joint moment in position P5].
Our results question the use of oversimplified approaches to compute joint resistance to mobilization.Certainly, some of the tested simplifications may remain arguable, justified by cost or time limitations.First, more www.nature.com/scientificreports/trained evaluators might manage to apply a movement profile to the dynamometer that strictly follows the joint plane (Case A) or limits the non-axial forces and moments (Case B and C), reducing the differences with the reference inverse dynamics method.Thus, it would require a training to compensate the differences of the cases A, B or C. Nevertheless, reference data (based on dynamometer kinematics or 3D load repartition during the mobilization) would be necessary for such a training.Second, the 3D inverse dynamics approach could be considered more complex to implement.In a clinical context, time is a valuable resource.However, we did not evaluate if 3D inverse dynamics approach takes more time than other literature methods.It should also be noted that some of the proposed methods in the literature require the measurement of specific distances, i.e. distance between the dynamometer or lad cell and the joint center that could also be time-consuming 8,9 .Moreover, the biomechanical models used in the literature methods differ from the ones used in gait analysis, which can lead to approximations when comparing the joint moments measured during passive stretches and the ones during gait 25 .With the reference inverse dynamics method, the biomechanical model can be the one used in gait analysis (use of the same marker set).Finally, it must be considered that differences in the joint resistance measure may lead to an altered estimation of mechanical properties of interest such as joint passive stiffness or musculo-tendon parameter such as muscle slack length 1,2 .For all these reasons, we recommend a complete kinematics and dynamics modelling appears before considering analyses for clinical purposes.

Limitations
Various limitations are present in this study.First, more than 16 subjects could have been tested; this could have led to a greater homogeneity of the sample.However, both box-plots and non-parametric tests were able to identify confidently the effects of the tested simplifications.Second, simplifications have been studied separately, whereas in the literature, there are several examples where these ones have been combined.For example, in most of the studies dynamometer kinematics (Case A) and distal segments movement (Case E) are simultaneously ignored.Nevertheless, since inverse dynamics is a linear problem, we expect that the resulting approximations add cumulatively.

Conclusion
Joint resistance to mobilization has been estimated in-vivo in several studies, through inverse dynamics computations, by measuring the applied forces and moments while manipulating the joint.Nevertheless, in most of the studies, simplified approaches are used to calculate joint moments.This study assessed the impact of these modelling assumptions.
Most of these assumptions lead to non-negligible differences with respect to a reference 3D inverse dynamics approach.Thus, literature results should be considered with caution.Finally, we recommend a complete kinematics and dynamics modelling before considering analyses for clinical purposes.

Figure 1 .
Figure 1.Test positions.Joints were mobilized through sagittal range of motion using a 3D handheld dynamometer.Dynamometer and body segment kinematics was measured with reflective markers and motion analysis system.Muscle activity was controlled with surface electromyography.

Figure 2 .
Figure 2. Free body diagram of the position P3 test.The joint resistance of the knee ϕ Shank/Thigh is calculated by applying Eq. (9) to the shank segment.A Shank represents the inertial action matrix and is computed from the subject's shank inertial parameters and kinematics.ϕ g/Shank is the weight action matrix.ϕ Dyn/Shank represents the measured forces and moments applied to the shank by the handheld dynamometer accordingly to their kinematics.ϕ Foot/Shank represents the forces of the foot on the shank segment computed with Eq. (10).

Figure 3 .
Figure 3. M θMax , M θMin and RMSD differences (10 −2 Nm/kg) of healthy and pathological subjects obtained for the different studied cases in all the positions.Asterisks indicate statically significant differences with the reference 3D inverse dynamics approach.

Figure 4 .
Figure 4. Relationships between the joint angle and moment of the different studied cases at low (LV) and high (HV) velocity for all the positions for one typically developed child.