Model-based bone material property identification

Abstract Correct torqueing of bone screws is important for orthopaedic surgery. Surgeons mainly tighten screws ad hoc, risking inappropriate torqueing. An adaptive torque-limiting screwdriver may be able to measure the torque-rotation response and use parameter identification of key material properties to recommend optimal torques. This paper analyses the identifiability and sensitivity of a model of the bone screwing process. The accuracy with which values of the Young modulus (E) of the bone were identified depended on the value of E, with larger values being less accurately identified. The error in identified σ u t s {\sigma _{uts}} (Tensile strength) values was less than 0.5 % over all the cases tested, with no discernible dependence on the co-identified values of E. Experimental validation is still required for the model and identification process, but this approach is feasible and promising from a theoretical perspective.


Introduction
A major part of orthopaedic surgery is the fixation of implants using bone screws. Correct torqueing of the screws is critical to prevent implant failure due to thread stripping [1] or screw loosening [2]. Implant failures may require revision surgery, with all associated costs and risks, or may cause further tissue damage and other complications [3].
Surgeons currently tighten screws ad hoc, without any specific torque guidelines. The success rate of this varies between surgeons, and incorrect tightening can easily occur [4]. It has been proposed that by monitoring the screwing process, the bone material properties can be identified using a model, and then used to predict the optimal tightening torque [5]; as these bone material properties are also dependent on factors like age and disease, this method takes these into account. This first requires an identifiable [6] model of the screwing process in terms of the most important material properties, and then a method for predicting optimal torque from these properties. This method is applicable to procedures where screws are fixed in bone. In many cases, the geometric placement of screws is more important than the correct torqueing. However, correct torque is still important to prevent tissue damage or screw loosening.
Other non-model-based methods have been proposed to regulate bone screw tightening torque. Reynolds et al. measures the plateau/steady-state torque while inserting a lag screw, and stops screwing when the torque reaches a multiple of this [7]. However, Reynolds' method might not be applicable to non-lag screws. Thomas et al. suggest monitoring the derivative of torque with respect to time, and halting the process when a spike is detected, indicating the screw has started binding [8]. This approach is effective at preventing overtightening but does little to ensure a minimum torque. Torque-limiting devices exist for some cases [9], especially in orthodontics, which uses small screws with low torque values [10], but these still require the surgeon to know and select the correct torque or limiter manually, in contrast to the proposed automatic model-based method.
Previously, a very simple model was demonstrated to be identifiable and robust for an arbitrary set of parameters [5]. This model was based on the material failure energy density of the bone, viscous friction between the screw and bone, and a threshold torque to advance the screw; this model was derived using an power balance between the input torque times speed and the material failure plus viscous friction. This simple model was significantly extended by Wilkie et al. [11], but the application of parameter identification was not investigated. This paper computationally analyses the model from Wilkie et al. [11], focusing on identifiability, recovery of the input variables in noisy data sets, and sensitivity to parameter values. The proposed model has greater physical descriptiveness than the model in Wilkie et al. [5], as specific material and geometric parameters rather than arbitrary lumped parameters are identified.

Model summary
The key model equations 1-11 were derived in Wilkie et al. [11] and have some basis in the model of Seneviratne et al. [12]. This model assumes that a simple screw design is selftapping into a pre-drilled cylindrical hole of greater diameter than the minor diameter of the screw (D h > D r in fig. 2), and is primary intended for trabecular bone. Equations 1-7 define constants used in the critical torque calculation in terms of the given geometry. r f and r s are the radii of the mid-point of the screw thread friction forces, and cutting forces, respectively. K f 0 is a lumped geometric parameter for friction forces dependant on screw and hole geometry. A c is the cross sectional area of the screw thread that displaces the hole material. θ is the thread helix angle. Equation 8 models how critical torque changes as the screw advances. Equation 9 models how the effective tor-sional spring constant (k) for a rigid screw in material of Young's modulus E changes with screw depth. Equations 10 and 11 model how angular velocity and torque are related to the input effort ϵ over time, using critical torque (T z (ϕ)) and k(ϕ) as functions of rotation (ϕ).
The value of ϵ is an input basis function that varies from 0-1 and represents the effort that the surgeon is exerting; this controls the torque and speed of the advancing screw by modulating a linear torque-speed relationship defined by a and b, shown in fig. 1. This was required to keep angular velocity defined in the absence of a viscous friction assumption [11]. The outputs are ϕ [rad], and T [Nm]. The model considers a critical torque (T z ) required to overcome friction/material resistance and advance the screw. Above T z , the screw advances according to the torque-speed relationship; below T z the screw does not advance, but can move with the bone as the bone elas-

Parameter identification overview
The overall parameter identification and sensitivity analysis process is outlined in Figure 3.

Input variable recovery
As it is not possible to measure the input effort of the surgeon directly, the effort must be determined from the measured angular speed and torque. It is not possible to do this exactly, as the relationship between T, ϕ, and ϵ depends on the unknown parameters, E and σ uts , but a heuristic method can be used. In this case it is assumed that a and b are known. The basic heuristic was to calculate the value of k at each time point using Hooke's law, k = ΔT/Δϕ. If the absolute value of the estimated k is relatively low, the position is changing a lot without much change in torque, indicating that the screw is advancing into the material, so the standard torque-speed model is used: ϵ = (φ + bT)/a [11]. Otherwise if the k value is higher than a threshold (3 N/rad was arbitrarily defined), that indicates the screw is only moving as torque changes, and therefore elastic behaviour is occurring, so the zero-speed assumption was used with the torque-speed model, giving ϵ = bT/a [11]. Due to high noise in T, it must be filtered before estimating k values; this was achieved with a 20 sample median filter, followed by a 20 Hz filter using the MATLAB "lowpass" function with a steepness of 0.8 and a stopband attenuation of 80 dB.

Time variant objective-function weighting
Changes in E have a very small effect on the T and ϕ response [11]. Hence, to ensure that E was not overfit to random noise, it was desirable to place more weighting on time periods that specifically exhibit elastic behaviour. This was done by changing the weight onφ for 0.1 s as the screw is loaded and de-loaded (where the elastic behaviour occurs), as per eqn. 12. The angular speed threshold of 0.07 [rad/s] was manually determined. The weights used were manually adjusted for the given sensor noise levels to give reasonable identification over a range of parameters (eqns. 12, 13, and 14). The objective function used for parameter identification used weighted 2-norms of the   [18] differences between the simulation (using current parameters) and data (original simulation with noise), shown in eqn. 15.

Identification methodology
Initially data is generated using the model in eqns. 1-11.
To simulate the effects of sensor noise and sampling, the model was first forward-simulated with a 20-cycle trapezoidal input signal ϵ(t) (ranging from 0 to 0.6, dwell of 0.5 s and rise/fall time of 0.25 s) at 10000 Hz, using the parameters in Tab. 1. Then the data was down sampled to 100 Hz to simulate sampling, and zero-mean Gaussian noise specified in Tab. 1 was added to simulate sensor noise. To simulate a gyroscope sensor, which would be a likely sensor in practice, the forward simulation angular position was numerically differentiated to get angular velocity, and the noise was added to this angular velocity signal. The parameter identification used only this noisy, down-sampled, angular velocity data. The identification process first takes the measured variables and performs the input recovery (Section 2.3) to get the approximated input effort over time. The input variables were then analysed to determine the time periods for the time-variant weighting (Section 2.4). Then simulated annealing [13] parameter identification was used to identify E and σ uts . The simulated annealing ran for 1000 iterations from a constant initial guess of (E, σ uts ) = (100[MPa], 1.5[MPa]), and allowed the objective function to move uphill up to 10 times in a row before being reset to the previous best parameter guess. Zero-mean Gaussian samples were used to perturb the parameters at each iteration. The standard deviations were a proportion of the parameters' current values. This proportion decreased exponentially over the 1000 iterations. Different rates were used for these exponential decays as σ uts converged much faster than E. The standard deviation proportion for σ uts went from 0.5 to 0.005, and the proportion for E went from 0.5 to 0.05.
The parameters/constants from Wilkie et al. [11] were used. The constants a and b were selected based on expected values. Variations in a and b are very unlikely to have a significant effect if they are consistent between the initial forward simulation, input variable recovery, and forward simulation during the identification; this is because a is mainly for scaling, and, with rearrangement, can always be combined with ϵ to write the equations in terms of aϵ, and b is usually multiplying T except for one case in eqn. 10.

Identifiability/sensitivity analysis
The effects of varying the two parameters (E and σ uts ) were investigated. For each step of each value, the parameter identification was run 100 times with different sets of sensor noise. The identified values for all runs of each parameter/noise level step were recorded, and the distributions of both parameters were plotted in constellations to show how the distribution of identified values changes as the parameters change. Principle component analysis [19] was used to find the principle components of the constellations, to analyse any correlations between E and σ uts identification. To test sensitivity, the two parameters were varied around their nominal values from Tab. 1. To quantify sensitivity, the 2-norm of the difference in data between the nominal and varied simulations was used. These methods allowed the model to be tested for identifiability and sensitivity, but not validity and accuracy, which requires experimental data.

Results
The result of the input effort recovery with default parameters is shown in Figure 4. Figure 5 shows how the distributions of the identified E and σ uts respond to changes in the true E and σ uts . Each distribution was normalised with the known correct value, and a logarithm was taken before plotting to make large and small variations legible side-by-side. The principle component analysis results for the constellations from     Figure 4 shows that the input variable recovery works as expected. The recovered results with the given settings closely match the known data. Figure 5 shows that the spread of the identified σ uts was minimal in all cases. This is further shown in the principle component analysis, where the ∼90°components relating to σ uts spread all have standard deviations of around 0.002 log(MPa), or about 0.5 %. Figure 5 also shows that it is possible to identify E to  a reasonable accuracy when it is a low value (3 MPa). However using the model in eqns. 1-11, E becomes hard to identify at higher, more physiological values such as 300 MPa (typical femoral head trabecular bone [14]), as there is very ambiguous estimation variability over many orders of magnitude. This identification inaccuracy is due to a higher Young's modulus leading to less significant elastic behaviour in the data, reducing the signal-to-noise ratio to almost zero. However, the orthogonality of the principle components in Tab. 2 shows that there was no discernible correlation between the identified σ uts and E values. Therefore, the poor identification of E did not affect the good identification of σ uts . Figure 7 and Tab. 3 also show that large changes in the value of E around an expected value for bone (300 MPa) had almost no effect on the output of Table 3: Quantitative differences for sensitivity analyses in Fig. 6 and Fig. 7. the simulation, suggesting that E is of minor significance. Whereas Fig. 6 and Tab. 3 show that even small changes in σ uts have significant changes on the output, suggesting σ uts is a much more clinically significant parameter than E. It should be noted that the default values for σ uts and E are plausible for the application of this method, but there is a large range of biologically significant values [20]. The model was demonstrated to be capable of identifying σ uts and E from noisy data in simulation, with plausible material properties. The model is especially promising for identifying σ uts , but quite limited for identifying E. Fortunately, the low sensitivity of the model output to changes in E suggests that accurate identification of E may not be very important to the bone screw-torqueing application, and the lack of correlation between E and σ uts identification shows that the poor identifiability of E is not detrimental to the identifiability of σ uts . Furthermore, it may also be possible to use a relation between σ uts and E [21] to infer E from the accurately identified σ uts and a priori information such as age, and apparent bone mineral density (e. g., from quantitative computed tomography imaging [22]).

Nominal Value Comparison Value
As the model was derived from physical principles, one may expect it to be a reasonable mathematical description of the physical phenomena. However, as mentioned in Wilkie et al. [11], there are a few assumptions used to simplify the model. Bone is a complex anisotropic, non-homogeneous, and discontinuous material, and this analysis modelled bone as an isotropic, homogeneous, and continuous material. This will result in some error in the identified parameter values. However even imperfect values may be useful if experimental testing shows that the error is sufficiently small. A noteworthy discrepancy is in the non-homogeneous structure of bone, often there is a harder outer layer and softer inner layer. This paper focused on the simple case of homogeneous bone, but the model, and specifically eqn. 8, could be modified to account for a varying σ uts with depth; as it is important to minimise variables for identification, this could be done by modelling σ uts in several discrete steps as a function of screw depth, or as two values (i. e., surface and bulk σ uts ), and a transition depth.
Due to noisy torque signals, filtering was required to remove noise before input variable recovery, however it also removed some desired higher frequency components of the signal. This can potentially introduce artefacts in the processed signal precisely in the regions governed by elasticity, and the parameter identification may over fit these artefacts. There are some notable fluctuations in the recovered data, some of this will be random, and due to noise, but there will be some systematic overshoots in the trapezoidal transitions due to the bandwidth limiting nature of the filtering approach. This was mitigated using median filtering, which can smooth out noise without smoothing out sharp transitions in the signal. However, the effectiveness will be dependent on the shape of the underlying signal. A possible solution might be an adaptive filtering span that reduces when sensor noise is low (by increasing cutoff frequency). Additionally, the threshold of k used for detecting elastic and non-elastic sections during input variable recovery was a constant in this paper, but there may be better heuristics for this. If the threshold is too high, it will treat non-elastic sections assuming elasticity, and if too low it will treat elastic sections with the non-elastic assumption.
The sampling rate and noise levels used were based on relatively low-cost general-purpose sensors (0.5 % of fullrange torque sensor such as NCTE-2300 Series [23], and 0.002 rads −1 noise at 92 Hz (∼100 Hz) from an MPU-9250 gyroscope [18]), using more accurate sensors with less noise and/or higher sampling rates may make it possible to identify higher Young's moduli with reasonable accuracy. Conversely, this paper only looks at simulated results, and experimental noise may be non-Gaussian, and may contain other corruption, such as gyroscope drift, torque hysteresis, sensor non-linearities, and temporal shifts across sensors.
The weights used for the different variables in the objective function were determined subjectively. More work is required to determine what the optimal weights are, and it may be possible to heuristically choose better values by analysing the measured data, and the tendencies of clinicians during the screwing process. These optimal values will also depend on the levels of sensor noise, as too much weight on a noisy sensor may result in overfitting to that sensors' noise and overwhelming the small weight of the other data. The threshold used for identifying timeperiods for time-variant weighting was chosen for simplicity in this analysis, and may not work well when the input is not a perfect repeating trapezoidal signal; however this is simply a form of edge detection, and there are potentially many other algorithms that could be used to make the process more robust.
Experimental verification is required to determine whether the modelling assumptions lead to unacceptable variation in the identified values; some variation is expected, but further work will determine how much variation there is, and how much is acceptable. Experimental verification is also important for the rest of the identification process, as the models of sensor noise may differ from reality, and other real-world limitations will apply, for example, measuring rotation at the handle end of a screwdriver will also measure any elastic behaviour due to torsion in the screwdriver shaft. If a significant relationship can be found between screwing speed and torque, all of the complexity introduced with the torque-speed model of the operator, and input variable recovery could potentially be eliminated. Chumakov [24] suggests that there is no significant relationship between screw speed and torque, although they used aluminium, which is not very strain rate dependant [25], while bone is [26], so more research is required.

Conclusion
This paper analysed the use of a bone screwing model for identifying material properties in simulation. Input variable recovery was successfully demonstrated, as was timevariant objective-function weighting. It was demonstrated that σ uts and E could be identified. σ uts was identified with great accuracy and consistency. E could be identified effectively at low values such as 3 MPa, however at more realistic values such as 300 MPa, the identified parameters had a prohibitively large variance; this may improve if more accurate sensors are assumed, or if other aspects of the modelling approach are improved/tuned, particularly in the heuristic sections of input recovery and time-variant weighting. Further study may also reveal E to be less significant than σ uts for optimal torque prediction, or it may be possible to determine E from σ uts using a priori relationships. Experimental testing is also required to validate all parts of the identification and model and ensure that they operate correctly under real-world conditions. Funding: Partial support by grants "CiD" and "Digitalisation in the OR" from BMBF (Project numbers 13FH5E02IA and 13FH5I05IA).