Parameter estimation of a model describing the human fingers

Abstract The goal of this paper is twofold: firstly, to provide a novel mathematical model that describes the kinematic chain of motion of the human fingers based on Lagrangian mechanics with four degrees of freedom and secondly, to estimate the model parameters using data from able‐bodied individuals. In the literature there are a variety of mathematical models that have been developed to describe the motion of the human finger. These models offer little to no information on the underlying mechanisms or corresponding equations of motion. Furthermore, these models do not provide information as to how they scale with different anthropometries. The data used here is generated using an experimental procedure that considers the free response motion of each finger segment with data captured via a motion capture system. The angular data collected are then filtered and fitted to a linear second‐order differential approximation of the equations of motion. The results of the study show that the free response motion of the segments is underdamped across flexion/extension and ad/abduction.


INTRODUCTION
Mathematical models are crucial in supporting our understanding of physical processes.To be effective, a mathematical model must be able to accurately describe real-life observations and be able to make robust and reliable predictions of yet unknown events.Furthermore, the assumptions used in the modelling must be clearly indicated and accounted for in the model and analysis.In the literature there are a wide variety of mathematical models that describe the kinematics of human fingers.The most influential models of the biomechanics of the upper extremity are presented in [1][2][3][4][5][6].These models have been extremely helpful in gaining an understanding of the motion of the finger segments and in providing results on the underlying forces of the musculature from experimental procedures.However, in these models, there is no discussion of subject-specific scaling of the parameters that influence the motion.To allow a mathematical model that describes the motion of the human fingers to scale with different anthropometries, an informed decision on the geometry of the fingers is required, which allows the determination of the mass of each segment.Moreover, the finger segments rotate about their respective joints; hence the determination of the individual masses and subsequent moments of inertia are important factors.In the literature, there are only a few papers that discuss the respective moments of inertia of the finger segments [7,8].However, in these papers, the equations used to obtain the moments of inertia are not included.Furthermore, the passive moment generated at each segment due to the passing tendons, synovial fluid, surrounding tissue etc. is assumed to be the same for everyone, and the values assigned to the constants describing this interaction are not referenced [4].The goal of this paper is therefore to use existing knowledge to further the development of mathematical models that describe the motion of the human finger segments, subject to different anthropometries, and to perform parameter estimation to determine the constants that describe the passive moment generated at each segment.A preliminary version of this work was reported at the BioMedEng 22 conference [9].

MODELLING
The equations of motion for human finger segments are established as a serial linkage of three compound cylindrical rod pendula with a density of  = 1.16 gr∕cm 3 [10].Applying Lagrangian mechanics with four degrees of freedom-three in the flexion/extension of the metacarpophalangeal (MCP), proximal interphalangeal (PIP), and distal interphalangeal (DIP) Conceptual image of a finger showing the three degrees of freedom in flexion/extension. 1 ,  2 ,  3 correspond to the flexion/extension angle of the proximal, middle, and distal segments, respectively.m i , I i , L i , T i corresponds to the mass, the moment of inertia, the length, and the overall torque of segment i. IReproduced with permission. [11]  FIGURE 2 The spherical coordinate system used in the modelling.
Angles  1 ,  2 ,  3 corresponds to the flexion angles of the MCP, PIP and DIP joints, respectively.xy 1 , xy 2 , xy 3 are the azimuthal projections of each end point.Angle  is the abduction angle of each finger.
joints, and one in the adduction/abduction of the MCP jointin a similar manner to the approach adopted in [1].A spherical coordinate system is used, with its origin attached to the MCP joint.The elevation angle corresponds to the flexion/extension motion of the segments, with its positive direction being perpendicular to the dorsal side of the palm facing away from it and the azimuthal direction to the ad/abduction.Figure 1 shows the flexion/extension angles in their respective frame of reference alongside the mass, moment of inertia, segment length and overall torque at each joint.
The convention used in this present model is that when all segments are fully extended, the elevation angle of the respective local reference frame is zero.The set of equations describing the system as a linkage of three rigid cylindrical rod pendula for the midpoint of each cylinder is shown in (1). Figure 2 shows the spherical coordinate system used in our model.
The variables in (1) are functions of time, but the time notation is omitted for convenience.Let  i , m i , L i ,  denote the flexion angle, mass, length of segment i (i = 1 to 3) and the corresponding abduction angle, respectively.Then the equations describing the position of the finger segments in the given spherical coordinate system are given by: x y 1 = − L 1 2 cos ( 1 ) ) ) ) x 1 = xy 1 cos () x 2 = xy 2 cos () x 3 = xy 3 cos () In ( 1) z i is the Cartesian representation of the projection of each point into the elevation plane, xy i is the projection of each point in the azimuthal plane and x i , y i are the individual Cartesian coordinates of the xy i vector (all for i = 1 to 3).
The moment of inertia, I o,i , of a uniform cylinder, of mass m i , length L i and radius R i , about its midpoint and the calculation of the mass of each segment are shown in (2) [12] The total kinetic energy of the system is given by (3): The gravitational potential energy is given by (4): In the literature, there are papers that include a passive moment generated at each joint from the surrounding tissue, tendons and synovial fluid and this has been modelled as a linear torsional spring and damper [4,8,13].The torsional spring effect arises in the potential energy term, as shown in (5): where K i , K a ,  eq,i ,  eq are: the torsional spring constant for the flexion/extension movement of segment i, the torsional spring constant for the abduction movement, the equilibrium angles of the flexion/extension of segment i and ad/abduction movements that correspond to the functional position of the hand, respectively [14].
The Lagrangian is expressed as the difference between the kinetic and potential energy of the system, that is, Lastly, the torsional damper effects are introduced into the Lagrangian as a Rayleigh dissipation function [15]: where B i , θi , B a , φ are the torsional damper constants of the flexion/extension and the flexion/extension angular velocity of segment i, and the torsional damper constant and velocity for the ad/abduction movement respectively of segment i.
Let q = { 1 ,  2 ,  3 , } denote the vector of the generalized coordinates, then the equations of motion for each degree of freedom are obtained via the Euler-Lagrange equations derived using the D'Alembert principle of virtual work [16] and given by: Here T i corresponds to the sum of the torques that is present during motion.During free movement, these are the muscle torques exerted at each degree of freedom and these are of the form: where F j is the muscle force production which is assumed to be based on a Hill-muscle model formulation, and the partial derivative corresponds to the moment arm of muscle j that spans the joint of segment i, with r being the musculotendon length as a function of the generalized coordinates [4,17,18].For the ad/abduction movement the muscles that span the MCP joint are considered as actuators.From ( 8), the well-known form of the equations of motion as shown in [4] is obtained.The main difference between the two models is that in the one presented here, an informed approximation of the geometry of the finger segments is used, which allows for a consistent determination of the moment of inertia of each segment based on the anthropometry of the individual.An approximation of the equations of motion for each degree of freedom can be obtained that can be solved analytically.This approximation is called the IBK approximation, where I denotes the moment of inertia, B the torsional damper and K the torsional spring constants respectively.In the literature, this type of model has been utilised to determine the spring and damper constants for different types of movements [19][20][21].The linear approximation of ( 8) is shown in (10) below∶ where I i , B i , K i are the moment of inertia, the torsional damper and the torsional spring constants for each degree of freedom.Before performing model parameter estimation, a structural identifiability analysis is required to determine whether all unknown parameters in the model can be identified uniquely from the given model observations/measurements [22].One important remark to note that will be useful in the identifiability analysis, is that the moment of inertia is uniquely determined from the first term of (8) for each degree of freedom.By measuring the length and radius of each finger segment, one can obtain an analytical expression for estimating its theoretical moment of inertia.The equations that derive the theoretical moment of inertia for the proximal (i = 1), middle (i = 2), and distal (i = 3) segments and for the abduction movement (i = a) are given below, respectively, as follows: ) I a = I 1 cos 2 ( 1 ) + I 2 cos 2 ( 2 ) + I 3 cos 2 ( 3 )

STRUCTURAL IDENTIFIABILITY
Structural identifiability analysis is concerned with the question of whether the parameters of a model can be uniquely determined in the theoretical situation of noise-free, continuous measurement/observation data [22].Regardless of the complexity of a model, if it is shown that the model is unidentifiable then it is not possible to determine the model parameters uniquely or locally, that is, there is an (uncountably) infinite set of model parameter values that can give rise to the same input-output problem.It is therefore important to undertake an identifiability analysis on the model of ( 10) as a prerequisite to estimating the parameters of the model and to support the design of an appropriate experimental procedure to ensure parameter identifiability for the given model outputs.
The free response of ( 10) is chosen as the model for our parameter estimation that corresponds to the loading/unloading of the torsional spring component.In real-life this is equivalent to the following type of motion: 1. Participants keep their hands in a relaxed position with no active movement or external force applied.2. The researcher moves each segment to an arbitrary position within its range of motion.3. The segment is then released.4. The response of the segment to its return to its natural steady state is recorded.
The structural identifiability analysis of the free response of (10) is determined using the Laplace transform approach [23].Dividing by the moment of inertia in (10), its free response is then given by: The initial conditions for each generalised coordinate will be the random initial angle within the range of motion of each segment q 0,i and a zero initial angular velocity.Taking the Laplace transform of ( 14) with these initial conditions yields the following Laplace transform of the response: From [23], each coefficient of the individual powers of  in the numerator and denominator of ( 16) is uniquely determined if and only if, the coefficient of the highest power of  in the denominator is one.From (16), the given response function coefficients (moment invariants) that can be uniquely determined are as follows: q 0,i , q eq,i , The model is therefore structurally unidentifiable since only the parameter combinations K i ∕I i and B i ∕I i are identifiable (can be uniquely determined) and not the individual parameters within these combinations.A-priori knowledge of one parameter is therefore required for all of the remaining parameters to be uniquely determined.Recall that the parameter I i for the moment of inertia can be estimated from (11)(12)(13)(14) for each segment.Thus, if the moment of inertia I i is known, then the remaining two parameters K i and B i can be determined uniquely for each segment from the given input-output relationship.

METHODS
The free response of each finger segment (index through little fingers) in the four degrees of freedom model was studied.
As mentioned previously, the free response movement corresponds to the involuntary motion of each segment, which is attributed to the loading/unloading of the torsional spring component of the passive moment.The free response movements of the thumb are not discussed in this paper.Participants were instructed to hold their forearm at a 90 • angle with respect to their arm.For the flexion movement, the participant's wrist was in its neutral position.For the abduction movement, their wrist was rotated to 90 • of pronation.For that wrist orientation, gravitational contributions are neglected since gravity is assumed to be perpendicular to the plane of motion.For the abduction movement participants were asked to keep their individual segments at full extension (all flexion angles set to zero), because the moment of inertia of the abduction movement is dependent on the individual segment flexion angles, as shown in (14).Throughout the trials, participants were instructed to have their fingers relaxed and not to apply any voluntary resistance.Great care was taken to reduce the movement of the proximal segments, minimizing the contributions of the Coriolis forces during the free response movement.A novel marker set was developed, as shown in Figure 3 that allows for the determination of the segment angles from motion capture data.Segment angles were determined using vector calculus.Each segment is represented by a vector between the proximal and distal markers for the joint.A vector between the mid-distance of the radial and ulnar styloid processes (the marker shown with a red arrow in Figure 3) and the distal marker to the MCP joint of each finger (shown with yellow and white arrows in Figure 3 is also used to determine the flexion/extension angle of each proximal segment.)Ad/abduction angles are determined between the vector defining the middle of the palm (mid-distance between ulnar and radial styloid marker (red arrow) to the distal marker of the middle segment's MCP joint (white arrow)) and the vector defined by the mid distance between the ulnar and the radial styloid marker (red arrow) and the proximal marker on the proximal segment (green arrows) of each finger.The L-shaped marker is used to define the anatomical planes of the dorsal side of the palm.Each vector is projected onto the sagittal (flexion/extension) or radioulnar (ad/abduction) planes and the dot product between the tail of the distal vector translated to the tail of the proximal vector, with respect to the joint of interest, is used to determine angle [24].Angle extraction from the motion capture data was performed in Vicon Nexus Procalc [24].
At least 15 trials were recorded per degree of freedom per finger.In total, each participant had at least 240 trials for the index to the little finger segments.Before marking the participant, the length and diameter of all their segments were measured and these values were used to determine the mass and moments of inertia of each segment.Motion capture data were recorded at 150 Hz using 12 infrared cameras in the Motion Capture Laboratory within the School of Engineering at the University of Warwick.A total of 23 able-bodied participants were recruited for the study.The exclusion criteria for the study were participants that had been diagnosed with arthritis, or had a digit amputation, or who were less than 18 years old.Written informed consent was taken before the experimental procedure.This study was granted full approval from the In Figure 4, the characteristic response of the free response movement where the finger segment returns to its steady state can be seen.In Figure 5, the characteristic response of the free response movement where the finger segment returns to its steady state can be seen.In Figure 6, the characteristic response of the free response movement where the finger segment returns to its steady state can be seen.In Figure 7, the characteristic response of the free response movement where the finger segment returns to its steady state can be seen.In Figure 8, the characteristic response of the free response movement where the finger segment returns to its steady state can be seen.The raw data of the free response of each degree of freedom can be seen.
For all fingers, except the middle finger, the abduction-free response movement of the whole digit was performed by moving each finger towards the middle one.This is evident from Figure 7, where the free response movement data are inverted  compared to the data shown in Figure 8, where the middle finger was moved away from the midline of the palm.
The angular data extracted from Procalc, were low-pass filtered using a fourth order Butterworth filter with a cut-off frequency of 15 Hz for the MCP and PIP joint movements and 12 Hz for the DIP joint and abduction movements.The cutoff frequencies were determined by plotting the Hilbert-Huang spectrum of randomly selected trials of the signal and averaging the highest frequencies that had the highest instantaneous energies, as shown in Figure 9. Instantaneous energy is proportional to the amplitude of the transformed signal and has units of (rad 2 ).Consecutively, the peaks of the data were identified, FIGURE 6 Index finger DIP joint raw data.The characteristic response of the free response movement where the finger segment returns to its steady state can be seen.

FIGURE 7
Index finger abduction movement raw data.The characteristic response of the free response movement where the finger segment returns to its steady state can be seen.and the model was fitted from the onset of the peak up to its steady state, where the natural response motion of the finger segments occurs.
Figure 10 shows the raw and filtered data alongside the fit with the lowest root mean square error for the same trial.

UNDERDAMPED MOTION
The underdamped motion of the free response is achieved when the discriminant of the characteristic polynomial of the denominator of ( 16) is negative.The resulting analytic expression is shown below: By fitting the analytic expression in (18) to the raw data, the unknown model parameters can then be determined from the following expressions:

CRITICALLY DAMPED MOTION
The critically damped motion of the free response is achieved when the discriminant of the characteristic polynomial of the denominator of ( 16) is zero.The resulting analytic expression is shown below.
By fitting the analytic expression of (20) to the raw data the unknown model parameters can then be determined from the following expressions:

OVERDAMPED MOTION
The overdamped motion of the free response is achieved when the discriminant of the characteristic polynomial of the denominator of ( 16) is positive.Let a, b be the roots of the corresponding characteristic polynomial with a > b.The resulting analytic expression is shown below: Since a, b are the roots of a second order polynomial then the parameters can be estimated from (22) using: A custom MATLAB code was written that performed the above-mentioned steps.The parameters were extracted for the fit that had the lowest root mean square error.For each segment, a mean and standard deviation for each parameter were obtained.Any missing parameters from the tables below were due to participants having a non-disqualifying injury at the specific segment or the number of available trials where the markers are shown being fewer than five.

RESULTS
Tables 1-4 show the mean values of the torsional spring and damper parameter estimates as determined from the analysis above for each segment for all participants, alongside the respective moments of inertia and damping ratios.
One particularly interesting outcome of the study, as can be seen from the raw data in Figures 4-8, is that the free response of the segments, for all of the degrees of freedom, oscillates before reaching its steady state.This behaviour is characteristic of underdamped motion.Thus, the average damping ratio was calculated for each segment for each degree of freedom using the following expression: In Table 5 the mean damping ratio values are provided for all of the trials.

DISCUSSION
The torsional spring and damper constants for all the degrees of freedom for the index through to the little finger segments were estimated.Comparison between the values determined here and those found in the literature can only be done for the       torsional spring constant of the DIP and MCP joints for the flexion/extension movements of the index finger.In the work presented in [21] a mean and standard deviation for the torsional spring constant for the middle segment were obtained using a robot that flexed the DIP joint of the index finger while the MCP joint of the same finger was held at either 0 or 60 degrees.In the experiment presented in [21] the authors used the same wrist orientation for the flexion/extension movement as presented here.From the data presented in [21] the extension stiffness value from their trial-to-trial reproducibility study is compared to the mean value of all the torsional spring constants of the PIP flexion of the index finger from Table 1.The extension value given in [21] was used because the data that were fitted to the free response equations correspond to the extension of the segments returning to their steady state.

Ring
In Table 6, the parameter values from this study and those found in [21] can be seen.The percentage difference between the two values is 29.08%.Even though the percentage difference is relatively high, this can be attributed to the difference in sample size and anthropometrical variation between the two studies.Another possible explanation for the percentage difference can be attributed to the Coriolis forces that come from the relative movement of the proximal and distal phalanges during the experiment.However, given the high standard deviation from our study, the mean value from [21] is within the boundary set by the standard deviation.
In the work presented in [25], the authors reported the passive MCP joint stiffness that was determined with two different methods using a standard stiffness measurement device and a soft robotic actuator.A similar wrist orientation was assumed in their study as that used in this paper.In the following table, the mean values together with the standard deviation (SD) for the results from [25] and the mean of the results presented in this paper are provided for comparison.  1 and the extension stiffness determined from the data available at [25].

Method K (N m/rad) SD (N m/rad)
Soft robotic actuator [25] 0.0324 0.0058 Stiffness measurement device [25] 0.0391 0.0057 This study 0.0366 0.0127 In Table 7, the parameter values from this study and those found in [21] can be seen.The percentage differences between the two different methods in [25] and this study are 12.18% for the soft robotic actuator and 6.61% for the stiffness measurement device.The percentage errors for the MCP joint are smaller compared to the PIP and this can be attributed to the fact that the contribution of the Coriolis forces on the proximal phalanx from the relative motion of the second metacarpal is diminished.However, from the previous two tables, it is apparent that the method presented in this paper provides results that are in good agreement with other methods found in the literature.This study has been shown to agree with different experiments for estimating the torsional spring constant of the passive moment of the fingers.The values found in Tables 1-4 have the potential to be used as a reference for the values of the passive moment contributions.The main advantage of this study is the reproducibility of the experimental procedure as it does not rely on any robotic actuators compared to those found in [21,25].Potential improvements to the experimental protocol might be to include a medical clamp to fix the elbow and wrist in place, which reduces fatigue for the participant, and potentially the use of an anesthetic agent that can numb the sensation on the fingers.The latter is believed to allow the involuntary movement of the free response to be exactly that, involuntary, without the participants having any reflex control over the moving segment.
Another interesting outcome of the study is the observed variation in the parameters across the recruited population.This is in contrast to the results presented in [4] where the authors assumed the same spring and damper constant values for all fingers for all participants.Given the extent to which these parameters vary, it is therefore important that appropriate scaling functions are derived in order to allow subject specific estimation of the passive moment parameters to be determined based on anthropometric data.On this basis, scaling functions between the torsional spring and damper parameters and the anthropometric values obtained were also determined.Linear fits between the estimated parameters and different anthropometrical measurements were examined.These measurements were the segment length (sl), segment radius (sr), the sum of the palm length (pl), palm breadth (pb), palm width (pw), the product between pl and pb as suggested in [26], the product between the hand length (hl), which is the sum of the palm length and the sum of the lengths of all the segments of the middle finger and pb.Lastly, the product of segment length and segment radius, the sum of the individual segment lengths (Ssl), the sum of the individual radii (Ssr), and the sum of the individual products of segment lengths and radii (SLR) were calculated (the last three were considered for the abduction movement only).For example, the SLR, Ssl, and Ssr for the index finger are given by, respectively: The normalized inverse of the standard deviation of the parameters was used as a weighting factor in the fitting process.Parameters with a low standard deviation had a higher weighting value than those with a high standard deviation.Fits that had the highest R 2 value were chosen.The scaling functions can be seen alongside the R 2 value below.The anthropometric parameters are given in, mm and the torsional spring and damper constants are in N m∕rad and N ms∕rad, respectively.In Table 8, the fitting coefficients with their associated R 2 value, alongside the equations used for scaling the parameters, can be found.
With the equations shown in Table 8 the torsional spring and damper parameters of the passive moment can be determined for any individual.
As mentioned, a particularly interesting finding of this study was that the free response of the finger segments is underdamped for all degrees of freedom.This is evident from the mean damping ratio values as shown in Table 5.This characteristic movement of the segments has not been reported previously in the literature.From the type of experiment performed, it is evident that these values are characteristic of the underlying structure of the human finger segments.Since this is an involuntary movement, it provides a glimpse of the passive structure and a quantification of the passive moment generated at each joint.An underdamped motion is characterized by lower force/torque production on a system compared to other types of free response, hence the characteristic oscillation about its steady state.However, since this has seemingly not been reported elsewhere in the literature, the conclusion is drawn here that, during voluntary motion of the finger segments, the passive moment increases and becomes either critically damped or overdamped.Potentially, this means that the passive moment components are dynamic in nature rather than static scalar parameters as considered previously in the literature.

CONCLUSION
In this paper, a novel mathematical description of the equations of motion of human finger segments based on Lagrangian mechanics and cylindrical approximations of the segments of human digits has been derived.A linear second-order differential equation approximation to the underlying equations of motion for each degree of freedom is established, and parameter estimation has been performed using this model for the free response movement.An experiment has been designed that corresponds to the free response movement, and a novel marker set has also been applied to obtain corresponding angular data.It has been shown that the free response of the finger segments is underdamped for all the degrees of freedom and scaling functions between the torsional spring and damper parameters and subject anthropometry have been derived.The intention is that this model will be used to support parameter estimation during voluntary movement which can be controlled by a musculoskeletal model with surface electromyography (sEMG) as its actuator.This has the potential to further our understanding of how passive moments change, if at all, during voluntary movement controlled by the central nervous system.Moreover, the model has the potential to be used in the design of a neuroprosthesis as a means to incorporate feedback control for joint state estimation from sEMG data and haptic feedback mechanisms for both partial and complete hand amputees.

FIGURE 3
FIGURE 3 Marker set developed for determining flexion/extension and ad/abduction angles from able-bodied participants.

FIGURE 4
FIGURE 4Index finger MCP joint movement raw data.The characteristic response of the free response movement where the finger segment returns to its steady state can be seen.

FIGURE 5
FIGURE 5Index finger PIP joint movement raw data.The characteristic response of the free response movement where the finger segment returns to its steady state can be seen.

FIGURE 8 FIGURE 9
FIGURE 8Middle finger abduction movement raw data.The characteristic response of the free response movement where the finger segment returns to its steady state can be seen.

FIGURE 10
FIGURE 10Raw and filtered data of the MCP joint movement of the ring finger of participant, alongside the best fit to the filtered data.The fit to the filtered data has a RMSE of 1.2 • and a R 2 value of 0.9927 suggesting an excellent fit.

TABLE 1
Index finger mean parameter estimate values for all degrees of freedom.

TABLE 2
Middle finger mean parameter estimate values for all degrees of freedom.

TABLE 3
Ring finger mean parameter estimate values for all degrees of freedom.

TABLE 4
Little finger mean parameter estimate values for all degrees of freedom.

TABLE 5
Mean values of the damping ratio for each segment in each degree of freedom alongside their standard deviation (SD).

TABLE 6
[21] value with standard deviation (SD) of the torsional spring constant from all trials of the PIP flexion as shown in Table1and the extension stiffness determined from the trial-to-trial study in[21].

TABLE 7
Mean value with standard deviation (SD) of the torsional spring constant from all trials of the MCP flexion as shown in Table