I recently conducted a study in which I developed a detailed spinal model to evaluate the dynamic characetrisitcs of a spinal disc implant design as it would perform in vivo. This model was validated through carefull comparison with publisehed data derived from cadaveric studies as well as other simulation models and in vivo intradiscal pressure measurements.
In addition to this a detailed FEM model of a cervical disc implant was created to which the inputs derived from the detailed spine simulaiton model could be imported and applied as boundary conditions. The slip velocities and pressures obtained from this analysis was then imported into a triangulation model which could calculate the estiamted wear rate over a typical lifetime of 10 000 000 cycles (as prescribed in the ISO standard for testing wear on spinal disc arthroplasty).
The aim of the study was thus to develop a chain of simulations in order to evaluate the wear imposed on a cervical disc implant during a lifetime of cycles and to evaluate how the obtained values compare to the wear obtained through the prescribed ISO test method. The following paragraphs summizes the work conducted.
The design of spinal disc replacements is challenging due to the complexity of the loads and motions imposed on the implants. Wear is a critical implant design factor, and wear testing evaluates the future efficacy of the device (Goel et al., 2006). Two wear test protocols are in general use today, i.e. ISO 18192-1 (2005) and ASTM F2423-05 (2006). They prescribe simultaneous cyclical rotations (axial, flexion/extension and lateral) and axial load profiles (Figure 1) that must be applied to the implant for ten million cycles at a rate of 1 Hz. The centre of rotation of the motions must correspond to the implant's centre of rotation. The implant must be immersed in a diluted bovine serum within a given temperature range for the duration of the testing.
Wear testing is expensive and time consuming (more than three months per test). Although some deductions can be made from the wear patterns, they provide limited insight for the implant designer into the mechanisms involved. The extent to which wear testing reproduces in vivo loads and motions, can also be questioned. Simulations can help to bridge this gap, but the in vivo wear of spinal disc implants has to date not been determined through a process of simulations, although it has been investigated in total knee arthroplasty after explantation (Fregly et al., 2005).
The chain of simulation models presented in this article is aimed at simulating the in vivo conditions as well as the wear tests, and can assist designers in understanding the loads, motions and wear mechanisms that spinal disc implants must contend with. They also provide a quicker route (compared to testing) for comparing alternative implant designs, even though numerical simulations involve various approximations.
The article's main focus is on describing how the individual simulations in the chain contribute to these objectives and how to transfer the data between the links of the chains. The main links of the chain are geometric modelling using reverse engineering (Materialise MIMICS) and CAD software (Autodesk Inventor), implant motion and load estimation using kinematic simulation (BRG LifeMOD), contact stress and velocity modelling using finite element modelling, FEM (Marc), and wear prediction using a linear wear model and reverse engineering triangulation software. The main elements in the chain of simulation models are described in the following sections and a critical evaluation of the chain as a whole in the discussion section.
Kinematic Model Inputs
The kinematic model's role in the simulation chain is to determine the loads and motions that the implant is exposed to. These parameters are strongly influenced by the motion of the head. Externally applied loads can also contribute, but are not considered in this paper.
Although studies have been conducted on the range of motion (ROM) of specific levels of the cervical spine under specific load conditions (Chin et al., 2003), most of these studies were conducted in vitro using muscle-less and/or ligament-less cadaveric cervical spines, in which spinal loads and kinematics are quite different from in vivo conditions (Miura et al., 2002). There remains a lack of published data regarding the actual in vivo range of motion, and the muscle and ligament loads responsible for the motion in the cervical spine (Goel et al., 2006). Even the wear test protocols use the findings of cadaveric studies (specifically those of White and Panjabi, 1978).
A kinematics simulation can potentially give more representative loads and motions at a cervical implant if a suitable range of motion for the body can be selected. This selection is, however, complex since it has to represent the implant's whole life (which will exhibit large variations amongst test subjects), and in some sense the "worst case" that the implant must withstand. The selection of a representative range of motion is not part of the simulation chain and therefore is beyond the scope of this paper, but De Jongh (2007), McGuan and Friedrichs (1996) and Jasiewicza, et al. (2007) describe procedures for measuring these motions.
Cervical Spine Kinematics Simulation
The kinematic model of the cervical spine with single-level disc replacement is used to determine the relative motions and loads that the implant is exposed to. LifeMOD, an extension of the dynamics simulation software program MSC.ADAMS, was used by the authors.
The kinematic model was created using LifeMOD’s standard GeBod anthropometric database. A median model was created (height = 1.778 m, mass = 77 kg), but modified to retain only the upper torso, central torso, lower torso, arms, neck and head (Figure 2). A further reduction of the model could be considered to further reduce the computation time, but the authors chose to rather include too much than too little, to ensure that constraints artificially imposed where the model is truncated do not influence the kinematics in the cervical spine.
Vertebrae, Muscles and Ligaments
Since LifeMOD assigns the neck as one segment, the various vertebrae segments had to be created manually. Computed Tomography (CT) scan data, with 2 mm spacing, of a complete cervical spinal unit (24 year old volunteer, 1.80 m, 90 kg) was imported into MIMICS (segmentation software) and converted to STL 3D models. A typical two part cervical disc implant design (Figure 3) was exported from Autodesk Inventor into MIMICS as STL files and the superior and inferior endplates were manually placed at the recommended positions (Wasserman, 2006) on the C5 and C6 vertebral bodies respectively, the most common implanted level (Wasserman, 2006). Each implant disc endplate was then merged with its vertebra, becoming a single geometric entity. Since the simulation is aimed at investigating the wear in the implant, the details of the interaction of the implant with the vertebra were not modelled. The merged parts were imported in STL format into LifeMOD as geometric entities, where they replaced LifeMOD's C5 and C6 shell element vertebrae. LifeMOD's shell elements of the other vertebrae were retained without changes. The 3 DOF joints (described below) between the other vertebrae and the imported entities were not affected by the small size differences between the LifeMOD's C5 and C6 shell elements and the imported ones.
LifeMOD generates the main muscle groups that is responsible for the control of large motions of each body segment, in this case flexion, extension and lateral motion. LifeMOD assigns muscle tissue properties, which include the physiological Cross Sectional Area (pCSA) and the maximum allowable stress in each muscle, from its database. Each muscle contains a contractile element in series with a spring-damper element, storing the input motion and effectively ‘training’ the muscles to reproduce the necessary force to recreate the desired motion. The pCSA values given by LifeMOD correspond to those given by Rezasoltani (1999) and Rezasoltani (2002), i.e. 155 to 244 mm2 for non-active subjects. The model with the muscle groups in place is shown in Figure 4. The generation of the cervical ligaments were performed manually since LifeMOD does not generate these soft tissues automatically. Biomechanical properties of the relevant cervical ligament tissues, as described by Yoganandan et al. (2000) were entered into LifeMOD. Ligament attachment points from cadaveric analysis were used (Smith et al., 1996). This anatomical data was verified by Van De Graaff (2002).
Intervertebral Joints
The intact intervertebral discs were represented by passive 6-degree-of-freedom bushing elements. The joints were placed on a line estimated to follow the Instantaneous Axes of Rotation (IAR) (Jansen and DeAngelo, undated; Bogduk and Mercer, 2000) of the cervical Functional Spinal Unit (FSU) as indicated in Figure 5. This aspect of the model neglects the movement of the IAR relative to the FSU when the FSU changes shape, since the exact location of the IAR and its relative motion remains unclear (Jansen and DeAngelo, undated; Bogduk and Mercer, 2000). The effect of the IAR location on the implant loads and motions can be investigated in future using the kinematic model, thus providing further insight.
The 6-degree-of freedom joint elements allow translations in the transverse plane (horizontal plane), which allows positional adaptation of the IAR to unbalanced forces. The stiffness properties of the intervertebral discs under compression were specified for the joint stiffnesses in al three directions, since no data was available for the other two directions (Yoganandan et al., 2000).
The interface definition between the implant's endplates, including contact stiffness (N/m), damping (N/m/s) and static friction coefficient, was based on the properties of cobalt chromium molybdenum, CoCrMo (Liu et al., 2005; Udofia et al., 2005; Soncini et al., 2003, Wei et al., 2005). An axi-symmetrical static FEM estimated the contact stiffness to be 106 N/mm (De Jongh, 2007). Damping is not expected to play a significant role in the kinematic simulation since the movements are relatively slow, but is required by LifeMOD for numerical stability.
It should be noted that in the kinematic model the whole of the C5 and C6 vertebrae, merged with the respective disc endplates, were defined as contact bodies, and thus includes the facet joints at the posterior end of the segments. The geometry as obtained from the CT scan was used for the facet joints, which may lead to errors when contact occurs at the facet joints, due to inaccuracies in the scanned shapes. Further, LifeMOD's shell elements of the remaining vertebrae do not make provision for the facet joints. The errors introduced by the geometric errors in the facet joints (that are modelled), and neglecting the contribution of facet joints of the remaining vertebrae to the kinematics, requires further investigation in future since some expect the errors to be small in the movement ranges considered here (Wasserman, 2007), while in vitro studies have indicated that the facet joints carry substantial loads in certain movement ranges (Goel et al., 1998).
Calibration and Kinematic Results
The muscle forces, ligament stiffnesses and their attachment points in the kinematic model were calibrated by imposing the measured head motion (as described in the preceding section) as input. This process trains the muscles and ligaments to generate the forces that replicate the input motion. To validate the simulation, its values for some key parameters computed for a typical range of motion (De Jongh, 2007) were compared to published data: The simulation predicted a peak intradiscal pressure of 3.03 MPa, which compares well with the 3.05 MPa reported by McGuan and Friedrichs (1998), and 3.10 MPa reported by Kumaresan et al. (1999). Further, Ha (2005) gives values of -1.5° to +1.5° for lateral rotation and -3° to +4° for extension/flexion, which is well matched by the simulation's prediction of 3° to +1.5° for lateral and -2.5° tot +4.7° for extension/flexion vertebral rotations. The simulation found the peak applied bending moments for the prescribed motion to be 3 Nm during flexion/extension and 4.5 Nm for lateral movement, which similar although somewhat larger than the 1.5 to 3 Nm reported by Goel et at. (2006), 1.8 Nm reported by Ha (2006) and 1.5 Nm reported by Yoganandan et al (2000). It should be noted that some of the differences in the reported quantitative values are due to differences in the range of motion used for the different studies.
The following parameters were extracted from the kinematic simulation results: intradiscal force/pressure at the C5/C6 level and the rotations of the vertebral bodies. Figure 1 shows a comparison of these parameters, for a typical ROM, with those prescribed in the ISO/DIS 18192-1 wear test protocols. Although the results shown in Figure 1 obtained from the kinematics simulation should not be considered to be definitive (since they were determined from a typical, but arbitrarily chosen, ROM), the large differences between the kinematics results and the wear test protocol's values indicates that there could be a large difference between the in vivo conditions and the wear test protocols' conditions.
The high rotational values of the test protocols may be due to their reliance on biomechanical experiments using cadavers, which exhibit no muscular or ligamentous reaction forces (Gilmour, 2006). In contrast to the range of motion, the maximum load prescribed by the test protocols is significantly lower than that found by the kinematic simulation (the latter agrees with Kumaresan et al., 1999). It should however be noted that this maximum force is for extreme flexion in the kinematic model, while the average force over the simulated time correlates better with the test protocols’ force input. These results emphasise that the selection of a representative range of motion, as input to the kinematic simulation, is not a simple matter.
Although LifeMOD "trains" the muscles to reproduce the ROM, it should be noted that it is not currently possible to validate the detailed kinematics since the in vivo tissue properties (i.e. ligaments and muscles) after surgical intervention may differ from the properties used for the simulation. Further, the specific design of the implant may introduce abnormal loading and motion conditions in the cervical spine because of possible restrictions in natural intervertebral motion (Ahn, 2005).
FEM Contact Analysis
The FEM contact analysis is the link in the simulation chain that determines the contact pressures and slip velocities required by the wear model. The analysis can use the loads and motions determined by kinematic modelling as described in the previous section, or those prescribed in a wear test protocol (ISO/DIS 18192-1 and ASTM F 2423 – 05).
For the results presented here, the central part of the CAD model of the implant (Figure 3) was discretized mostly using tetrahedral elements. A “multi-point constraint” (MPC) was applied to the superior endplate from where the motions could be applied. For kinematic model inputs, the MPC was situated where the motions were determined in LifeMOD, while for standard wear test protocol inputs the MPC was situated at the centre of rotation of the implant in the neutral position as described in the test protocol (Figure 6). The external forces were applied to the inferior endplate, and the endplate was constrained with an inferior/superior and translational directional stiffness of 829.7 N/mm (Yoganandan et al., 2000).
The material properties of CoCrMo, E = 210 GPa and ν = 0.3 (Lui et al., 2005; Udofia et al., 2004), were used for the model. A dynamic friction coefficient of 0.07 was taken from an experiment conducted by Soncini et al. (2003). The principle contact stresses given by an analytical Hertz stress analysis of a sphere in a spherical cavity (giving values from 12 to 21 MPa) were used to validate those from the FEM (8-20 MPa) for the interface between the superior half and the inferior half of a neutrally orientated implant with a vertically applied force. This showed that the FEM results are sufficiently accurate to be used to evaluate the simulation chain, because its sole purpose is to provide contact pressure and slip velocity inputs to the wear modelling.
Using the finite element models, the contacting pressure and slip velocity histories were determined for two sets of inputs, i.e. the ISO wear test protocol and those provided by the kinematic simulation.
Wear Estimation
Wear rate is proportional to the contact pressure and slip velocity (Link et al., 2004; Shigley et al., 2004; Fregly et al., 2005; Yan et al., 2006; Field et al., 2007), which are estimated as described above. The proportionality constant is called the steady state wear factor and for CoCrMo it is approximately 1.3 10-6 mm3/N-1m-1 (Yan et al., 2006; Field et al., 2007).
The wear rate ideally should be calculated on the contacting surface at each nodal point of the FEM at each time increment in the 107 ROM cycles (the number prescribed by the ISO protocol). By integrating the wear rate over time, the wear depth at each nodal point for the entire lifetime of cycles can be estimated. Redistribution of contact stresses will occur because the contact stresses will gradually reduce in areas of high wear rates, leading to increased stresses in other areas. To account for this effect, the FEM analysis would have to be repeated at regular intervals during the life time, each analysis using a geometry updated with the wear up to that time step. Since the focus of the work presented in this paper was to investigate the viability of the chain of simulations, and because of the long run time required for each finite element analysis, redistribution was not taken into account in the results presented here. A further reason for not taking the redistribution into account is that the available validation data was for the accumulated wear over the whole contact surface. The average wear rate is less sensitive (compared to local wear rate values) to the redistribution of contact pressure since a reduction in wear rate at one position (due to a reduction in contact pressure), will necessarily be compensated for partially by an increase somewhere else (where the contact pressure had to increase).
The preceding procedure gives the wear depth at each node on the contacting surface, thus giving non-uniform values at discrete points. To quantify the volumetric wear distribution over the entire contact surface, the following method was used: The contact surfaces for pre-wear and post-wear were triangulated by inputting the nodal coordinates of the contact surface (after subtracting the wear from the Z-coordinate for the post-wear triangulation case) into a triangulation algorithm of reverse engineering software. Using equations for triangular based pyramids, the volume of each triangular element was calculated for pre-wear and post-wear, the difference giving the volumetric wear.
The simulated ISO test's wear rate distribution of the inferior surface is shown in Figure 7 with the Z-changes proportional to the initial wear rate. To better visualise the initial wear rate patterns, the initial phase nodal wear depths were imported into Patran as displacement data and merged with their respective contact nodes. Figure 8 shows the resulting contour map.
No published experimental wear test data suitable for validation of the simulation chain's wear predictions was available. The closest currently available test data is that of Metcalf et al. (2004) and Anderson et al. (2004). The former's data is for a hip implant (with different load conditions), but with materials and surface diameters similar to the simulation presented here. They report a volumetric wear rate of 1 mm3 per million cycles. By comparison, the simulation's wear volume (using the ISO test protocol's inputs and the initial wear rate for the whole life), was 1.18 mm3 per million cycles. Anderson et al. (2004) investigated cervical disc implants, but with titanium on polyurethane contact. They found a percentage weight loss of 1.76%, while the simulation gave 1.92%. Given the differences between the experimental and simulation conditions, no firm conclusions about the simulation's accuracy can be made, but the comparison still shows that the simulation results are realistic. conditions as well as the wear tests, and can assist designers in understanding the loads, motions and wear mechanisms that spinal disc implants must contend with. They also provide a quicker route (compared to testing) for comparing alternative implant designs, even though numerical simulations involve various approximations.
This study demonstrates that the chain of simulations presented here can provide additional insight into the loads that a cervical disc implant is subjected to in vivo and can predict the wear that can be expected sufficiently accurate for comparative studies. The result of each simulation in the chain was transferred to the next simulation, thus demonstrating that the “linking” of the chain of simulations was possible.
The simulation of the kinematics of the implant is the link of the chain most reliant on judgement and the most difficult to calibrate, in particular in the "training" of the muscles, the selection of the IAR and the compliance of the joints. The modelling of the face joints' contribution to the kinematics also needs further consideration.
The FEM is the link of the chain that uses a well established approach and can therefore be used with confidence, subject to the usual element selection, convergence and constraint considerations. The wear model is also well established and can be expected to predict wear rates commensurate with the accuracy of the contact stresses and velocities. Even though redistribution of wear was not demonstrated in this paper, it is a straight forward addition to the procedure reported.
It should be noted that the effects of approximations or selections made in the previous simulations will accumulate in the wear predictions. Results from simulation starting with a selected ROM could therefore differ substantially from that obtained when simulating a wear test protocol.
It is not an objective of this study to attempt to find an alternative for the wear test protocols, but to find a way to reduce design related time and cost through comparative simulation studies. Also, the sensitivity of wear to effects such as the placement of the implant and body posture, can be investigated. The authors suggest that this objective was adequately met by the methods proposed in this paper.
In a similar study Fregly et al. (2005) used a chain of simulations to calculate the wear of a TKA (Total Knee Arthroplasty) and compared the results with that of a retrieved prosthesis. The methods used in their study resemble the methods used in this study. The results in both the TKA investigation and the work presented here showed good comparison to the retrieved prosthesis or published values, respectively. It should however be kept in mind that motions in the spine are more complex and it is not possible to define a definitive cyclical motion as was the case in the TKA study.
References
Ahn, H.S., 2005. A virtual model of the human cervical spine for physics-based simulation and applications. Ph.D. dissertation, The University of Tennessee, Memphis.
Anderson, P.A., Sasso, R.C., Rouleau, J.P., Carlson, C.S. and Goffin, J., 2004. The Bryan Cervical Disc: wear properties and early clinical results. The Spine Journal 4, 303s-309s.
ASTM, 2006. Standard Guide for Functional, Kinematic, and Wear Assessment of Total Disc Prostheses. ASTM F2423-05.
Bogduk, S. and Mercer, N., 2000. Biomechanics of the cervical spine. I: Normal kinematics. Clinical Biomechanics 15, 633–648.
Chin, K.H., Tan, K.W., Goh, J.C.H., Toh, S.L. and Lee, V.S.P., 2003. Flexibility testing of the human upper cervical spine under continuous loading and unloading. Summer Bioengineering Conference. Sonesta Beach Resort, Key Biscayne, Florida, 1185-1186.
De Jongh, C, 2007. Critical Evaluation of Predictive Modelling of a Cervical Disc Design. MScEng (Mechanical) thesis, Stellenbosch University, South Africa.
Fregly, B.J., Swayer, W.G., Harman, M.K. and Banks, S.A., 2005. Computational wear prediction of a total knee replacement from in vivo kinematics. Journal of Biomechanics 38, 305-314.
Field, S.K., Cooke, K., Jones, A.H.S. and Teer, D.G., Feb 2007. Sputtered Coatings for Medical Applications. Available at: http://www.teercoatings.co.uk/index.php?page=papers [13 April 2008]
Gilmour, L.J., 2006. Biomechanical Testing of Two-Level Cervical Disc Arthroplasty. MSc thesis, The University of Tennessee, Memphis.
Goel, V.K., Panjabi, M.M., Patwardhan, A.G., Dooris, A.P. and Serhan, H., 2006. Test protocols for evaluation of spinal implants. Journal of Bone and Joint Surgery 88, 103-109.
Goel, V.K., Vijay, K. and Clausen, J.D., 1998. Prediction of load sharing among spinal components of a C5-C6 motion segment using the Finite Element Approach. SPINE 23, 684-691.
Ha, S.K., 2006. Finite element modelling of multi-level cervical spinal segments (C3-C6) and biomechanical analysis of elastomer-type prosthetic disc. Medical Engineering & Physics 8, 534-541.
ISO, 2005. Implants for Surgery – Wear of Total Intervertebral Spinal Disc Prostheses – Part 1: Loading and Displacement Parameters for Wear Testing and Corresponding Environmental Conditions for Tests. Draft International Standard ISO/DIS 18192-1.
Jansen, T.H. and DeAngelo, D.J., (undated). Influences of the Location of Vertebral “IAR” on Cervical Spine Kinematics. School of Biomedical Eng., The University of Tennessee-Memphis.
Jasiewicza, J.M., Treleaven, J., Condie, P., and Jull, G., 2007. Wireless orientation sensors: Their suitability to measure head movement for neck pain assessment. Manual Therapy 12, 380–385.
Kumaresan, S., Yoganandan, N., Pintar, F.A. and Maiman, D.J., 1999. Finite element modelling of the cervical spine: role of intervertebral disc under axial and eccentric loads. Medical Engineering & Physics 21, 689-700.
Link, H.D., McAfee, P.C. and Pimenta, L., 2004. Choosing a cervical disc replacement. The Spine Journal 4, 294S-302S.
Liu, F., Udofia, I.J., Jin, Z.M., Hirt, F., Rieker, C., Roberts, P. and Grigoris, P., 2005. Comparison of contact mechanics between a total hip replacement and a hip resurfacing with metal-on-metal articulation. Journal of Mechanical Engineering Science 219, Part C, 727-732.
McGuan, S. and Friedrichs, A., 1998. Exploring Human Adaptation using Ultrasonic Measurement Techniques and Optimized, Dynamic Human Neck Models. Mechanical Dynamics Inc., Friedrich-Schiller-Universität Jena.
Metcalf, J.E.P., Cawley, J. and Band, T.J., 2004. Cobalt chromium molybdenum metal-on-metal resurfacing orthopaedic hip devices. Business Briefing, Medical Device Manufacturing & Technology.
Miura, T., Panjabi, M.M. and Cripton, P.A., 2002. A method to simulate in vivo cervical spine kinematics using in vitro compressive preload. SPINE 27, no. 1, 43-48.
Nightingale, R.W., Chancey, V.C., Ottaviano, D., Luck, J.F., Tran, L., Prange, M. and Meyers, B.S., 2007. Flexion and extension structural properties and strengths for male cervical spine segments. Journal of Biomechanics 40, 535-542.
Rezasoltani, A., Ylinen, J. and Vihko, V., 2002. Isometric cervical extension force and dimensions of semispinalis capitis muscle. Journal of Rehabilitation Research and Development 39-3, 423-428.
Shigley, J.E., Mischke, C.R. and Budynas, R.G., 2004. Mechanical Engineering Design. 7th ed, McGraw-Hill, New York, pp. 652-653.
Soncini, M., Redaelli, A. and Montevecchi, F.M., 2003. FE dynamic analysis of a knee joint prosthesis during the walking cycle. Summer Bioengineering Conference, Sonesta Beach Resort, Key Biscayne, Florida, 0547-0548.
Smith, L.K., Weiss, E.L. and Lehmkuhl, L.D., 1996. Brunnstrom’s Clinical Kinesiology. 5th ed., Philadelphia.
Udofia, J., Yew, A. and Jin, Z.M., 2005. Contact mechanics analysis of metal-on-metal hip resurfacing prostheses. Proc. IMechE Part H: Journal of Mechanical Engineering Science 218, 293-305.
Van De Graaff, K.M., 2002. Human Anatomy. McGraw-Hill, New York, pp. 158–257.
Wasserman, J., 2006. Personal Communication. Neurosurgeon, Wilgeheuwel Private Hospital, Randburg, South Africa.
Wei, R., Booker, T., Rincon, C. and Arps, J., 2005. High-intensity plasma ion nitriding of orthopedic materials Part I: Tribological study. Surface & Coatings Technology 186, 305-313.
White, A.A. and Panjabi, M.M., 1978. Clinical Biomechanics of the Spine. Philadelphia, p. 22.
Yan, Y., Neville, A. and Dowsen, D., 2006. Biotribocorrosion - an appraisal of the time dependence of wear and corrosion interactions: I. The role of corrosion. Journal of Physics D: Applied Physics 39, 3200-3205.
Yoganandan, N., Kumaresan, S. and Pintar, F.A., 2000. Biomechanics of the cervical spine Part 2. Cervical spine soft tissue responses and biomechanical modelling. Clinical Biomechanics 16, 1–27.
