Trabecular bone in the bird knee responds with high sensitivity to changes in load orientation

Since the late 19th century, it has been suggested that the thin struts and plates in trabecular bone adapt dynamically to meet the demands of mechanical loading (von Meyer, 1867;Wolff, 1982). Specifically, Wolff(1892) proposed that trabecular struts grow to align with the trajectory of the predominant loads experienced within a joint. Although this hypothesis, often termed `Wolff's law', has been widely employed, few controlled experiments have tested if the orientation of trabeculae changes as a consequence of altered loading patterns. Here, we present an experimental test of Wolff's law for trabecular orientation using a within-species, age-matched design.

Previous studies have correlated trabecular architecture with known or predicted orientations of loading using a diverse array of data and methods,with mixed results. Several modeling studies have shown that the arcuate trajectories of trabeculae in the human femur and calcaneus tend to be aligned with the orientation of peak compressive stresses predicted from finite-element models (e.g. Carter and Beaupre, 2002; Gefen and Seliktar, 2004). In addition, strain gauge studies that examined trabecular orientations relative to strain orientations measured in vivo, indicate a close correspondence between the preferred alignment of trabeculae and the orientation of peak compressive strains in the calcaneum of poteroos (Biewener et al.,1996) and the zygomatic arch of pigs(Teng et al., 1997). Additionally, comparative and observational studies indicate the predominant trabecular orientation in bones such as the human patella, the cervid calcaneus, and the equine ischium match well with their likely pattern of loading (Currey, 2002; Skedros et al., 1999). Finally, a number of more clinically focused studies have examined the effect of pathologies and surgically based, tissue engineering approaches on trabecular architecture. In general, these studies find that the mechanical environment of a joint has considerable effects on mechanisms of tissue repair relevant to trabecular bone(Hollister et al., 2001;Smith-Adeline et al., 2002; Goldstein, 2004; Papaloucas et al., 2004).

One limitation with these studies, however, is that a correspondence between known and predicted trabecular trajectories does not address to what extent trabecular bone responds dynamically to its mechanical environment or how other factors, such as growth, affect the spatial distribution of trabecular bone. Bertram and Swartz(1991) articulated the limitations of the available empirical evidence supporting Wolff's law for trajectorial orientation, arguing that the processes that organize trabeculae in growing or fractured bone might well differ from those in non-pathological adult bone. The lack of empirical evidence linking changes in mechanical environment to subsequent changes in trabecular architecture in healthy adults(Bertram and Swartz, 1991; Cowin, 2001) suggests Wolff's law for trabecular bone warrants further scrutiny.

We have experimentally tested the trajectorial theory component of Wolff's law in the distal femur by altering the orientation of joint forces in the knee under normal conditions of loading while holding constant other potentially confounding variables such as species differences, age, activity level and growth. Birds running on an inclined treadmill used a more flexed knee during stance phase than those exercising on a horizontal treadmill,thereby changing the presumed orientation of compressive forces in the knee. If, as proposed by Wolff and others, trabeculae adapt to their mechanical environment dynamically to align with the orientation of compressive force,then this change should result in a corresponding shift in trabecular architecture. We have tested the prediction that the between-group difference in the orientation of the thickest, most numerous, and/or densest trabecular struts in the distal femur is equivalent to the presumed difference in the orientation of peak compressive force.

Beginning at 14-days post-hatching, 13 guinea fowl Numida meleagris Linnaeus 1766 were divided into three treatment groups. Two groups were exercised on a treadmill (10 min/day; 6 days/week; 45-day treatment period) at either a 0° (i.e. horizontal; N=6) or a 20° (N=5) incline (Fig. 1A). Running speed was increased as the birds grew to maintain a constant Froude number (Alexander and Jayes, 1983) of 0.7-0.8 (∼0.85 m s^-1 by day 45). A third control group (N=2) was not exercised. All birds were housed in a 9 m² room, with ad libitum food and water. IACUC approval was obtained for all procedures prior to the study.

On day 45 of the study, birds from the Level (N=2) and Incline(N=3) groups were run over trackways set at pitches that matched exercise conditions (Level: 0°, Incline: 20°) in order to measure knee flexion and ground reaction forces (GRF) during stance phase. During these trials the positions of the greater trochanter, the tibiotarsus-tarsometatarsus joint (TTj), and the tarsometatarsus-phalangeal joint (TPj) were recorded using a high-speed infrared camera system(Qualysis®; Qualysis Motion Capture Systems, Gothenburg, Sweden) operating at 240 Hz. Anatomical landmarks were determined by palpation and small (5 mm)reflective markers were adhered to the overlying skin in order to track their position. Simultaneous ground force was recorded at 1000 Hz with a custom-built force plate, embedded in the trackway, which uses four rosette strain gauges (two fore, two aft) to measure vertical ground force. Normal(i.e. perpendicular to the trackway) GRF was measured by summing traces from fore and aft gauges. Force traces were passed twice through a Butterworth filter and the resulting traces summed to produce the normal ground reaction force (nGRF) trace.

To determine the angle of knee flexion at peak nGRF, force plate and kinematic recordings were compared. Because only a small number of strides(N=8) yielded clean nGRF traces and adequate kinematic data, we performed a kinematic analysis to determine which gait element coincided with peak nGRF (Fig. 2). A sample of 16 strides provided clean nGRF traces with simultaneous marker data for the distal limb. For unimodal force traces (N=10 strides) and one of the peaks in bimodal force traces (N=6 strides), foot-off, defined as the first kinematic frame in which the TPj rises from its lowest point(Fig. 2A), was found to be coincident with peak nGRF (mean difference = 0 s; Fig. 2B). Foot-off was therefore used as the indicator of peak nGRF in determining joint angles at peak joint reaction force for trials in which nGRF data was inadequate,increasing the number of strides that could be included to measure the mean angle of knee flexion (N=26). The position of the greater trochanter and TTj at foot-off were used to determine the angle of flexion at the knee trigonometrically, using femur and tibiotarsus lengths measured with digital calipers after birds were sacrificed and these elements harvested.

The surface of the femoral condyle that is in contact with the proximal tibiotarsus is a function of knee flexion: with greater flexion, the tibial plateau slides posterior-inferiorly along the femoral condyles. The orientation, relative to the distal femur in the parasagittal plane, of compressive joint reaction forces (JRF) transmitted via the tibial plateau must therefore be a function of knee flexion. Thus, although we were unable to measure JRF in vivo, we assumed that any between-group difference in knee angle at peak nGRF resulted in an equivalent difference in the orientation of peak compressive JRF relative to the femoral condyle in the parasagittal plane (Fig. 1B). We therefore expected any difference in knee angle to result in a similar difference in trabecular orientation in the femoral condyle. The limitations of this approach for estimating between-group differences in JRF trajectory are discussed below.

At the end of the treatment period, birds were sacrificed and their hindlimb bones were removed and cleaned using dermestid beetles and a 1%bleach solution. Computed microtomography scans of the femora were then obtained using a Skyscan 1072 100 kV microtomograph (Aartselaar, Belgium). Whole distal femora were scanned at 19.43 μm isotropic resolution at 100 kV and 98 μA with three-frame averaging. Raw data were reconstructed using a cone-beam algorithm. Reconstructed cross-sectional images were then filtered using a two-dimensional median filter (ImageJ, NIH) to reduce background noise.

We restricted our analysis of trabecular architecture to the dense,relatively homogenous, ∼3 mm thick layer of spongiosa underlying the articular surface (see Fig. 5B)because it makes up the majority of the epiphysis and because JRFs transmitted from the tibiotarsus must first traverse these trabeculae. The trabeculae in the spongiosa are formed within the epiphysis as it ossifies, and are generally arrayed radially from the center of the epiphysis outward (see Fig. 5B).

Established three-dimensional techniques used to measure trabecular orientation directly from the micro-computed tomography (micro-CT) volume,such as volume orientation, mean intercept length (MIL), and star-volume analysis (Odgaard, 1997, 2001) could not be used effectively here for two reasons. First, the resolution of the scan (19.43μm) did not permit the accurate binarization of the highly diffuse spongiosa in the microCT volume because of partial-volume effects, preventing established methods (Odgaard, 1997, 2001) from producing reliable measures of orientation. Second, and more importantly, the subchondral layer of spongiosa being examined does not present a spherical region of investigation. Instead, as in many hinge-like joints, these spongiosa lie in a distinct, nearly two-dimensional subchondral layer, bound proximally by a region of thicker trabeculae and distally by the articular surface (see Fig. 5B). The thickness of this layer is markedly smaller than the radius of curvature of the articular surface (see Fig. 5B), and the arc of the articular surface in the sagittal plane is far greater than in the coronal plane. Creating a spherical volume of interest to examine this spongiosa layer is therefore problematic; spheres large enough to encompass all of the relevant spongiosa inevitably include unwanted regions that compromise analysis.

Instead, the distribution of trabecular density within the spongiosa relative to the long-axis of the femur in the parasagittal plane was measured using radon transform analysis (RTA; Deans,1983). The two-dimensional RTA treats an image as a matrix of grayscale values (black=0, white=255). The values for all pixels in a given row are summed, creating a column of row-intensity values; this process is repeated iteratively as the image is rotated 179 times by 1°(Fig. 3A). The resulting columns of row-intensity values are arrayed by angle of rotation(Fig. 3B) and summed, producing a density distribution (Fig. 3C). Because the highest row-intensity values are produced when row trajectory corresponds with the orientation of the most optically dense objects in an image (Fig. 3A,B), the peaks in the density distribution reflect the orientation of maximum bone density (Fig. 3C). This measure of trabecular density combines strut orientation, size, number and localized mineral density, which are traits that determine the strength of a trabecular volume(Odgaard, 2001). Whereas other techniques are available for assessing trabecular orientation in planar images(e.g. Kuo and Carter, 1991;Oxnard, 1993; Herrera, 2001),the radon transform is optimal for assessing the orientation of peak trabecular density in a radial array of struts, and has the advantage of being particularly robust to noise within an image(Deans, 1983). As shown in Fig. 3, RTA correctly determines the orientation of the thickest struts when they are radially arrayed (Fig. 3D) and is robust to induced noise in the image and the removal of the image center(Fig. 3E,F), but will not produce a false signal if no primary orientation exists(Fig. 3G).

In addition to testing the RTA for two-dimensional images with known strut orientation (Fig. 3D-G), we compared results from the RTA method against those calculated via the mean intercept length (MIL) method (Odgaard,1997). Trabecular orientation in a sample of human femoral micro-CT volumes (N=4) was measured independently using the MIL and RTA methods. Orientations given by these two methods were similar(N=8, r²=0.97, P<0.01) when the image analyzed via RTA was roughly coplanar with the main axes of trabeculae (Fig. 4). However,because the RTA method employed here is inherently two-dimensional, for images in which the predominant orientation is perpendicular to the plane being analyzed, MIL and RTA methods do not agree(Fig. 4). Whereas the RTA method cannot detect the primary orientation of trabeculae orthogonal to the plane of analysis, it accurately and reliably determines the orientation of trabeculae in a given plane, and is comparable to other methods of analysis when that plane is coincident with the primary orientation of trabeculae. It is therefore ideal for the present study, as the plane of interest(parasagittal) was determined a priori, and the majority of trabeculae in the condyle lie roughly in this plane.

To measure trabecular orientation, or, more specifically, the orientation of peak trabecular density (OPTD), the isotropic μCT image stack was loaded into a volume analysis program (VGStudio Max 1.2, Volume Graphics®,Heidelberg, Germany) and oriented into true mediolateral view. Four parasagittal 2-dimensional images, spaced 0.97 mm apart, were then extracted from the lateral condyle and examined using RTA. Because RTA requires a circular image (Deans, 1983), a half-circle region of interest (ROI) with a perimeter that matched the articular surface of the condyle (Fig. 5B) was taken from each two-dimensional slice, and a duplicate image, rotated 180°, was added to the original. The center of this circular ROI, comprising a small number of larger trabeculae, was removed in order to restrict the analysis to the spongiosa. Analyzing such ring-shaped images requires analysis of only the central region of the intensity-angle matrix in order to calculate the density distribution(Fig. 6). The long axis of the distal femur, determined for each parasagittal slice as the anterior profile of the distal right femoral shaft in lateral view, was used to establish 0° for RTA analysis.

Four intensity-angle matrices, one from each parasagittal slice, were summed to produce a composite density distribution for each subject. This summation requires a control for bone density between slices, as differences in bone density affect row-intensity values, creating the possibility that slices with the greatest trabecular density would be weighted more heavily in the composite density distribution. In order to normalize bone density between slices, intensity-angle matrices were converted to grayscale and then equalized in Adobe Photoshop 7.0®. Doing this scales the intensity values within the image so that the highest values are white (255) and the lowest values are black (0). Scaling the intensity-matrices via this method prevents slices with greater bone density from being weighted more heavily while preserving the signal strength (peak intensity/background intensity) for each image. Similarly, in producing a mean density distribution for each group(grayscale), equalized intensity-angle matrices from each slice (N=4 slices × group size) were summed. Mean density distributions for each treatment group were calculated by summing the density distributions for all subjects in each group, thereby weighting the contribution of each individual by the strength of its signal.

Upon examining the micro-CT volumes for each subject, it was apparent that one bird in the incline group had experienced bone resorption in the femoral condyle, perhaps as a result of calcium reuptake related to egg production. The resulting condyle showed several large pockets within the layer of trabeculae being analyzed. Although including this individual in the analysis does not significantly affect the incline group OPTD value, the subject was removed from analysis in order to restrict comparisons to individuals not undergoing bone resorption.

Controlling for size between subjects is unnecessary when intensity-angle matrices are compared, as these matrices are bounded between 0-180°regardless of size. By contrast, controlling for differences in signal strength (i.e. signal amplitude) of the density distribution is useful in determining a mean OPTD for each condition, because peaks of lower-amplitude signals are more likely to be affected by noise in the image and other sources of measurement error. Peak amplitude is therefore a measure of signal quality. In order to weight each subject appropriately by signal amplitude, each density distribution was divided by peak density for that subject so that peak density for all subjects equaled 1. Mean density distributions for each treatment group were then calculated by summing the density distributions for all subjects in each group, thereby weighting the contribution of each individual by the strength of its signal. The resulting group density distributions were used to compare joint angle determined by kinematics (see above) to group OPTD (Fig. 7).

As expected, knee angle at peak nGRF in the Level group(76.3±1.33° mean ± s.d., N=16 strides) was greater than in the Incline group (62.6±3.52°, N=10, P<0.01 Student's two-tailed t-test, Fig. 1). This result is consistent with previous work reporting that incline running results in more-flexed limb postures during stance phase(Jayne and Irschick, 1999). Assuming the trajectory of the JRF relative to the proximal tibiotarsus was the same in both groups, we predicted this 13.7° difference in knee flexion would result in a corresponding difference in the orientation of the thickest trabecular struts in the lateral femoral condyle. That is, we predicted the orientation of peak trabecular density (OPTD) would be∼14° more inferocaudal with respect to the distal femur, in the sagittal plane.

As predicted by Wolff's law, mean OPTD (the peak of the summed density distributions of each group) reflected differences in joint angles between treatment groups (Fig. 4A,B). Mean OPTD calculated for the incline group was 99.6° relative to the long axis of the distal femur, significantly more acute, by 13.6°(P<0.01, Student's t-test, Table 1) compared with 112.9° for the level and 116.0° for the control groups(Fig. 7A,B). Remarkably, the 13.6° difference in OPTD between the incline and level groups is nearly identical to the 13.7° difference observed in joint angle at peak GRF. The small 3.1° difference in mean OPTD between the Level and Control groups is consistent with the prediction that these groups would have similar trabecular orientation (Fig. 7A,B). OPTD values correspond to the angle between the long axes of the femur and tibiotarsus at peak GRF (Fig. 1C) such that they were perpendicular to the tibial plateau during peak compressive loading, suggesting the thickest, densest and/or most numerous trabecular structures were aligned with the orientation of peak compressive loads in the distal knee.

Individual variation in OPTD was apparent within each treatment group(Table 1). Although within-group variation in OPTD was less in the Incline group than in the Level and Control groups (Fig. 5),this difference was not significant (P>0.05, F-test). Although direct comparisons with previous studies is complicated by the use of different methodologies, the degree of individual variation in OPTD in this sample appears to be consistent with that reported for primary trabecular orientation in other species (e.g. Ryan and Ketchum, 2005).

Results from this study indicate that between-group differences in trabecular orientation in the distal femur matched differences in knee-flexion at peak JRF, suggesting that trabecular architecture is both responsive and highly sensitive to its mechanical environment. Although some recent studies have questioned many previously assumed relationships between form and function in the skeleton (Cowin,2001; Lieberman et al.,2003) and their underlying mechanisms(Bertram and Swartz, 1991),these results support the hypothesis that trabecular struts grow to align with the orientation of peak compressive forces, at least in growing juveniles. Further, these results lend support to previous studies proposing the orientation of trabeculae in adults is a function of repeated loading during ontogeny (Hert, 1992), and that trabecular orientation reflects the orientation of peak compressive stress (Biewener et al., 1996; Teng et al., 1997; Carter and Beaupre, 2002). How bone struts sense the orientation of strains and then respond appropriately,however, remains unknown, although there is some indication from tissue engineering studies that the struts are particularly responsive to tensile strains (e.g. Smith-Adeline et al., 2002; Goldstein, 2004). Similarly, the radon transform approach used here does not enable us to distinguish the relative contributions of trabecular density and thickness to OPTD.

Because this study examined growing juveniles, the results cannot address the sensitivity of trabecular architecture in adults. Previous work suggests that the ability of adult bone to remodel in response to applied loads is diminished not only in cortical bone(Lieberman et al., 2003; Pearson and Lieberman, 2004),but also in trabecular bone (Christiansen et al., 2000; Keaveny and Yeh,2002; Knopp et al.,2005). Similarly, these results may not be directly applicable to instances of fracture or other pathology, as subjects were healthy and experienced physiologically normal levels of strain.

Another important limitation of this study was the use of gross differences in joint kinematics to establish between-group differences in peak JRF. The construction and validation of a finite-element model of the guinea fowl knee was beyond the scope of this study. Thus, while the differences in knee flexion between groups was marked and consistent(Fig. 1B), it is not possible in the present analysis to demonstrate conclusively that JRF trajectory differed to a similar extent. Similar studies employing FEM to determine JRF trajectory, such as those of Carter and Beaupre(2002) and Gefen and Seliktar(2004) will no doubt further elucidate the relationship between JRF trajectory and trabecular bone architecture.

Because our study used a within-species design comparing age-matched,non-pathological subjects experiencing physiologically normal levels of bone strain, our results strongly suggest that differences in trabecular architecture can result solely from differences in loading. This result, in turn, supports previous studies suggesting that interspecific differences in trabecular architecture (Fajardo and Muller, 2001; Currey,2002; Ryan and Ketchum,2005) result from differences in locomotor loading regimes. Still,potential differences in the response of trabecular bone in juveniles versus adults must be considered. Locomotor behavior, and presumably loading regimes,may change considerably throughout ontogeny, especially in species with a large suite of habitual locomotor behaviors (see Doran, 1997). Trabecular architecture in adults may reflect, to some extent, habitual loading experienced during ontogeny.

One interesting possibility raised by these results is that OPTD may help assess the orientation of habitual peak joint reaction forces, and hence habitual loading patterns, in species for which joint kinematics during locomotion are unknown. This method may therefore provide a useful complement to strain-gauge studies, as gauges are difficult to employ in joints in vivo. Similarly, the radon transform technique employed here may be useful for establishing limb postures and locomotor behaviors in extinct species in which well-preserved fossils are available. For example, this approach could shed light on current debates regarding whether fossil bipeds such as Tyrannosaurus rex(Hutchinson, 2004) and Australopithecus afarensis (Ward,2002) used extended or flexed hindlimbs. Such analyses may benefit by focusing on subchondral spongiosa as considered here, particularly for joints such as the knee in which larger, secondary struts not in direct contact with the articular surface constitute much of the joint volume. Subchondral spongiosa, although technically challenging to analyze,experiences joint forces most directly and should therefore exhibit the strongest architectural response to differences in loading regimes.

While the above results are encouraging, future research is needed to determine how aging influences the dynamic relationship evident between the orientation of compressive loading and trabecular distribution, and the extent to which various aspects of trabecular architecture (e.g. strut thickness,connectivity and orientation) change in response to loading. In addition, a better understanding of how trabecular bone adapts itself to peak compressive loading will require more detailed information on the distribution of in vivo loads within the joint along with insights on the mechanotransductive mechanisms by which the cellular populations responsible for bone maintenance, formation and resorption sense and respond to strain orientations.

List of abbreviations
 
• JRF

  joint reaction force

 
• micro-CT

  micro-computed tomography

 
• MIL

  mean intercept length

 
• nGRF

  normal ground reaction force

 
• OPTD

  orientation of peak trabecular density

 
• ROI

  region of interest

 
• RTA

  radon transform analysis

 
• Tpj

  tarsometatarsus-phalangeal joint

 
• Ttj

  tibiotarsus-tarsometatarsus joint

This study was supported by funding from Harvard University and the American School of Prehistoric Research to D.E.L., a NSF Graduate Fellowship to H.P. and an Izaak Walton Killam Memorial Scholarship (Alberta Prov. CIHR)to D.M.L.C. We thank E. Tytell and P. Madden for help in developing the RTA method, R. Main for constructing the force plate, and P. Martinez for animal care. The 3D Morphometrics Center (Calgary) is supported by NSERC grant 238992-02, CFI grant no. 3923 and Alberta Innovation and Science grant no. URSI-01-103-RI to B.H. We thank two anonymous reviewers for useful comments on a previous draft.