A preliminary case study of the effect of shoe-wearing on the biomechanics of a horse’s foot

Subjects

One Thoroughbred healthy male adult horse (540 kg body mass) from the Royal Veterinary College (RVC) participated in the study. The horse had previously been trained for and participated in locomotor studies in the laboratory. Fifteen minutes of training were provided for the horse to adapt to the experimental setup. The study was reviewed and approved by the Royal Veterinary College’s Ethics and Welfare Committee (approval number URN 2011 1094).

Data collection

Each trial lasted two to four seconds, during which the horse was led across a custom-designed platform (Fig. 1A). A custom-designed platform rather than a treadmill was used for this study because our methodological approach for the in vivo estimation of the intersegmental forces for FEA required accurate ground reaction force (GRF) measurements. Accurate measurements of all GRF components could not be obtained using any available treadmill.

A Sony HDR (Sony, London, UK) high definition video camera was placed perpendicular to the platform to approximate walking speed (25 Hz). To obtain foot kinematics, two custom X-ray fluoroscopes (RSTechnics, Netherlands; refurbished Phillips systems, 36 cm intensifier; ≤110 kV, ≤3 mA) were retrofitted with two AOS high-speed digital cameras (AOS Technologies AG, Switzerland) to acquire biplanar fluoroscopy images at 250 Hz of the horse’s feet as it walked through an undistorted and calibrated capture volume (∼30 cm per cube edge) located on the forceplate.

Standard videos introduce geometric distortion to the images due to lens misalignments (Dobbert, 2005) and the fluoroscopes create nonlinear pincushion and spatial artefacts (Wang & Blackburn, 2000), which degrade image quality. To remove this distortion, we imaged a standardized sheet of perforated metal grid (1.6 mm thick with 3 mm diameter holes, spaced 4 mm apart in a 60° staggered pattern) (RS Components Ltd, UK) in front of the intensifiers and captured the grid images at 32 frames for 1 s. The exposure setting for the two intensifiers during the collection of the grid images was respectively set at 57 kv, 38 mA and 50 kv, 19 mA.

To correct the distorted video images and the calibration cube images, we created undistortion files from the standardised reference grid images using the XrayProject MatLab (Mathworks Inc. Natick, Massachusetts, USA) protocol (www.xromm.org). By using the reference grid images of each intensifier as templates, we created a uniform set of squares and applied it to all video images (Brainerd et al., 2010). The distortion correction algorithm compared the spacing between holes in the fluoroscope image with the idealized spacing and calculated a transformation matrix for correcting all video images. We used a local weighted mean (LWM) distortion correction algorithm, implemented in MatLab.

To calibrate the video images and generate a DLT algorithm in order to create “virtual cameras” in the Maya virtual 3D space for rotoscopy of the horse’s foot, we followed the calibration process described by Brainerd et al. (2010). We placed a custom-designed calibration cube with 64 steel 3 mm diameter spheres as calibration points within the field of view of the two intensifiers and recorded the calibration images at 32 Hz for 1 s (Fig. S1) (Brainerd et al., 2010). The calibration cube consisted of acrylic sheets of uniform thickness (5.42 mm), on which we drilled 16 holes of 3 mm diameter in a square pattern of 65 mm separation. Nylon pillars served as spacers between the tiers. The exposure settings for the two intensifiers during the collection of the calibration cube images were set at 46 kv, 39 mA and 71 kv, 11 mA.

The undistorted images of the calibration cubes for each camera were used as templates to calibrate all video images using protocols established by Brainerd et al. (2010). The undistorted calibration cubes were digitised by manually selecting the centroid of a minimum of 12 calibration cube beads that could be seen in both cameras using MatLab. Digitizing was repeated until a calibration coefficient of low residuals (<0.3) was given for each camera. The resulting direct linear transformations (DLT) were used to optimise the focal point and axes of the cameras. The camera data were then imported to Maya software (Autodesk, San Rafael, California, USA) to recreate a virtual scene with similar coordinates to the experiment (Brainerd et al., 2010).

Exposure settings were set to 69 kV, 53 mA and 72 kV, 54 mA for the two sources, which were placed on the left side of the platform. Each intensifier was placed 2 m from its corresponding X-ray source on the right side of the platform such that the beams were oriented horizontally at approximately 90° to one another (Fig. 1A).

Kinetic data were collected simultaneously at a rate of 1,000 Hz using a forceplate (60 × 90 cm with Hall Effect sensors, 2,000 lb maximum vertical force; AMTI, Watertown, MA, USA). Prior to analysis, the forceplate data were low-pass filtered using a 4th order zero-lag Butterworth filter with a cutoff frequency of 15 Hz. All data were synchronized to the fluoroscope system.

The unshod horse was guided 344 times across the experimental platform over a period of two weeks. Following the end of the experiments for the unshod condition, the horse’s forefeet received mild trimming and were each fitted with a stainless steel shoe with toe clips (5 inches wide) and 6 nails. The identical procedure was then followed to guide the shod horse over the platform 65 times on a subsequent day. The difference in trial numbers between the unshod (344 strides) and the shod (65 strides) conditions was due to the large number of spatially incomplete data for the former. Strides that were spatially incomplete (i.e., the right forefoot only stepped partially within the capture volume) and/or unsteady (i.e., with evident deceleration and acceleration following observation of the video images during data collection) were excluded from further analysis. Four steps from the shod and four steps from the unshod right forefoot that were spatially complete and steady were processed using the markerless XROMM (X-ray Reconstruction of Moving Morphology: Brainerd et al., 2010; Gatesy et al., 2010) workflow to construct a model and obtain 3D joint rotations and translations. The limited number of steps per conditions is a limitation of the XROMM approach when used in live animal studies for species as large as a Thoroughbred horse. For a step to be valid during the XROMM procedure, the animal has to step within the calibrated field of view without deviation. If there were minimal deviations from the capture volume, we were unable to visualise the distal right forefoot in both cameras in order to extract the 3D joint kinematics.

Model construction

The horse was euthanized at the end of the experiment for unrelated studies and its right forelimb was removed and frozen (–20°C). Computed tomography scans (GE Lightspeed 16-detector unit; General Electric) were used to obtain the three dimensional (3D) skeletal geometry of the horse’s forefoot (slice thickness 0.625 mm, 0.460 pixels mm⁻¹, 512 × 512 pixel images, 620 slices). These data were then processed to extract solid 3D polygonal mesh objects in Mimics (version 16.0; Materialise, Inc., Leuven, Belgium) and then imported into Maya (Autodesk, San Rafael, California, USA) to construct the biomechanical models’ segments (Fig. 1B). Four segments were defined: the metacarpus (MC), first phalanx (P1), intermediate phalanx (P2) with the navicular bone and the distal phalanx (P3). An articulated skeleton was then created by hierarchically linking these segments (Gatesy et al., 2010) into a kinematic chain representing the metacarpophalangeal (MCP), proximal interphalangeal (PIP) and distal interphalangeal (DIP) joints (Fig. 2).

Joint orientations and positions were defined by first positioning all bone segments into a neutral anatomical pose (forefoot lying fully horizontally). Cylinders were then visually fit to the joint surfaces (i.e., the distal epiphyses of the MC, P1 and P2) to identify the axes of joint rotations. These locations were confirmed when manual manipulation of the virtual joint resulted in a natural motion where adjacent bones did not interpenetrate each other. Dissected cadaveric specimens and plastic models were used to further assess joint locations and positions. Transformations between coordinate systems were defined using a Z (green axis), Y′ (yellow axis) and X″ (red axis) Cardan rotation, respectively representing abduction-adduction, flexion-extension and long axis rotation (Fig. 2). Axes were defined so that positive joint angles represented adduction, flexion and external long axis rotation.

Markerless XROMM

Trajectories for each joint were quantified using established protocols (Brown University, USA; www.xromm.org) for scientific rotoscoping (markerless XROMM) (Gatesy et al., 2010; Baier & Gatesy, 2013; Baier, Gatesy & Dial, 2013; Nyakature & Fischer, 2010). In brief, markerless XROMM is a technique that allows one to quantify 3D motion by animating model segments (i.e., 3D polygonal mesh objects) to match postures observed in experimental X-ray video images (Fig. 1B). For each experimental X-ray trial, the horse foot model was aligned with the bone X-ray silhouettes in undistorted and calibrated video images using the anatomical features of each bone as reference guides (Movies S1 and S2). Joint transformations (i.e., joint rotations and translations) were then extracted from the model for the MCP, PIP and DIP joints (Fig. 2). The MCP kinematic data were excluded from further analysis because the midshaft and proximal epiphysis were out of the field of view for most of the stride.

All steps for shod (n = 4) and unshod (n = 4) conditions were used to measure joint kinematics and foot kinetics, but a single representative step was selected for each condition (shod n = 1; unshod n = 1) for the subsequent mid-stance inverse dynamics calculations and FEA used to estimate bone stresses in the right forefoot digit.

Inverse dynamics

Inverse dynamics methods were used to calculate the intersegmental forces (required for estimating bone stresses) at mid-stance for the two steps selected for FEA (shod and unshod conditions). Larger loads can occur at other points during the stance phase but the focus of this study is methodological and for simplicity we modelled solely mid-stance.

To perform the analysis, the skeletal model created for scientific rotoscoping was recreated in Software for Interactive Musculoskeletal Modeling (SIMM; Musculographics, California, USA). To create an exact replica of the original model, bone geometry was exported directly from Maya and imported into SIMM, where the geometry was reassembled by reproducing the original joint structure (i.e., number of joints and degrees of freedom). The model was then imported into OpenSim (Delp et al., 2007), which has a built-in routine to perform inverse dynamics.

Mid-stance was defined as the point halfway between foot strike and toe-off gait events, which were determined from the vertical GRF data. Mid-stance joint angles were exported directly from the XROMM workflow and used to position the model in OpenSim. Mid-stance GRFs, obtained from the synchronized forceplate data, were transformed into the same reference frame as the OpenSim foot model and then applied to the distal phalanx (P3). Data integrity between the motion and force data were verified by visually inspecting the location of the centre of pressure (CoP) (from forceplate data) relative to the foot placement (from XROMM kinematics) using OpenSim. Because the distal segments are small (i.e., low mass) and walking is a slow motion (i.e., low joint accelerations), gravity and inertial forces were assumed to be negligible relative to the forces resulting from the applied mid-stance GRFs. Thus, mass and inertial properties for all segments were not included in the analysis. OpenSim’s inverse dynamics and joint reaction force routines were then used to calculate the intersegmental joint forces and moments acting on the segments at each joint. These data were expressed in the local frame of the segment and used as inputs into the FEA.

Finite element analysis

For each phalanx (P1, P2 (including the navicular bone) and P3), a separate finite element model was created in Abaqus/CAE, software version 6.13 (Dassault Systemes Simulia Corp, Providence, Rhode Island, USA). The corresponding intersegmental forces from the inverse dynamics analysis were applied to each bone independently and stress was determined using the Abaqus/Standard implicit direct default solver.

Speed data

The mean walking speed of the shod and unshod conditions was at 0.72 ms⁻¹ and 0.76 ms⁻¹ respectively (Table 1). This corresponded to a Froude number—which is a dimensionless representation of movement speed (Alexander & Jayes, 1983; Fr = velocity^{2 ∗} (9.81 ms^{−2 ∗} hip height)⁻¹)—of 0.05 for the shod condition and 0.06 for the unshod condition (Table 1), indicative of a slow walk. The footfall patterns also maintained the typical lateral sequence found in walking (Data S1).

Kinetic data

The GRF data from the shod (n = 4) and unshod (n = 4) conditions during the stance phase of locomotion are shown in Fig. 3 and Data S2. In all directions the force pattern was quite similar between the shod and unshod conditions (Fig. 3). The results from the cross-correlation analysis assessing the correlation between the unshod and shod horse in overall GRF patterns showed that there was a high positive correlation between the craniocaudal GRF patterns (with a maximum correlation coefficient of 0.994 with a 2% lag of the shod pattern). These results also showed a high positive correlation between the vertical GRF patterns (with a maximum correlation coefficient of 0.996 with a 0% lag of the shod pattern). While there also seemed to be a large similarity in mediolateral GRF patterns, the strength of the positive correlation was lower than in the other directions (with a maximum correlation coefficient of 0.747 with an 11% lag of the shod pattern).

The craniocaudal GRF at the beginning of the stance phase was directed caudally, which then shifted cranially from mid-stance until the end of the stance phase of locomotion (Fig. 3). The maximum cranial GRF for the shod horse was shown at approximately 75% of stance (497 N) and for the unshod horse at 78% of stance (396 N). Between 75–78% of stance, the shod horse showed on average a 21% higher craniocaudal GRF than the unshod horse, yet this difference was not assessed for significance.

The mediolateral GRF for the shod condition was directed medially throughout the whole stance phase and reached an average maximum force of 85 N after mid-stance. Contrastingly, the mediolateral GRF for the unshod horse was directed medially during the first 10% of the stance phase, shifted laterally until late mid-stance, and then was redirected medially until the end of the stance phase. The highest mediolateral GRF for the unshod horse was laterally directed and occurred before mid-stance (∼40–45%), reaching approximately 100 N.

There was a strong similarity in the vertical GRF pattern between the shod and unshod conditions during the stance phase. However, at mid-stance, the vertical GRF of the shod condition (3,195 N) was approximately 10% higher than the unshod condition (2,888 N).

Kinematic data

The kinematic data for the shod (n = 4) and unshod (n = 4) conditions during the stance phase of locomotion are shown in Fig. 4 and Data S3. The results from the cross-correlation analysis between the unshod and shod conditions showed a high positive correlation for the PIP mean joint angle (with a maximum correlation coefficient of 0.989 with a 0% lag of the shod pattern). A strong positive correlation was also found for the DIP mean joint angle (with a maximum correlation coefficient of 0.975 with a 0% lag of the shod pattern).

We describe some differences in kinematic patterns here for the shod vs. unshod conditions but it is very important to note that none of these have been tested for true statistical significance, because of the small sample sizes. Overall, in both the shod and unshod conditions, the PIP joint was in a flexed posture early in the stance phase (1040%) and ended in an extended posture before mid-stance (45% of stance) until the end of the stance phase (90%). Between 10–90% of the stance phase, the PIP joint for the shod condition had an average range of motion (ROM) of 14 degrees. The average ROM for the unshod condition’s PIP joint was estimated at an average of 13 degrees.

The maximum difference between the DIP joint angle in the shod and unshod conditions was at 40% of stance, when the unshod horse DIP approximates a neutral posture (∼=0°), while the shod horse DIP exhibited a more flexed posture. At mid-stance and until the end of the stance phase, the DIP extended in both the shod and unshod conditions (Table 2).

Finite element analysis (FEA)

The intersegmental forces obtained from the inverse dynamics technique and applied to the different segment models (P1, P2 and P3) during the FEA are presented in Table 3. Our FE analysis showed that the shoe increased the concentration of von Mises and maximum principal stresses on the dorsal (Fig. 5 and Fig. S3) and palmar (Fig. 6 and Fig. S4) aspects of the distal (P1, P2, P3) bones of the horse’s forefoot; however this possible increase has not been assessed for significance due to the small sample size. Comparisons of the principal stress patterns showed that in both conditions the dorsal (Figs. S3 and S5) and palmar aspects (Figs. S4 and S6) of all phalanges respectively underwent tension and compression.

The sensitivity analysis on the constraints for the P1 showed that fixing the constraints on the distal joint surface reduced von Mises stresses compared to fixing the centroid nodes and constraining the rest of the nodes to align with the long axis of the bone. However, in both cases bending occurred craniocaudally (Fig. S7).

Our sensitivity study of the P3 constraints showed that an increase in the surface area of the constrained nodes reduced von Mises stresses dorsally (Fig. S8).