Level-wise differences in in vivo lateral bending moment are associated with microstructural alterations in bovine caudal intervertebral discs

Intervertebral discs (IVDs) serve the dual roles of transmitting large forces between adjacent vertebral bones while allowing flexibility of the vertebral column. During activities of daily living, lumbar IVDs are subjected to large axial forces and bending moments resulting from the carriage of upper body weight and muscle contractions. While these forces produce large fluid pressurization in the disc's inner nucleus pulposus (NP) (Wilke et al., 1999), fiber bundles in the outer annulus fibrosus (AF) are subjected to tensile stresses (Schmidt et al., 2007).

Bovine caudal IVDs are frequently used as a mechanical analog to that of the human lumbar spine (Demers et al., 2004). They have been used to study in vitro disc biology (Alini et al., 2008; Chan et al., 2011; Roberts et al., 2008), as a source of IVD cells (Horner et al., 2002; Kraus et al., 2017; Kraus and Lufkin, 2016) and as a pre-clinical model to assess treatments for disc degeneration (Chang et al., 2010; Likhitpanichkul et al., 2014; Teixeira et al., 2016). Unlike the complex and highly variable loads applied to human lumbar discs by a combination of body weight and muscle contraction, bovine caudal discs are loaded only by muscles. The articular processes on adjacent bovine caudal vertebrae do not contact each other to form facet joints, which have been shown to support loads in human (Stokes, 1988) and quadruped (Smit, 2002) lumbar spines. Additionally, the bovine tail is moved only in lateral bending (with a frequency of 5–40 events per minute; Dougherty et al., 1995) and dorsal extension (with a frequency of 0.01 events per minute; Miedema et al., 2011). While the tail may be moved in combined lateral bending and flexion–extension, it lacks the obliquely oriented muscles required to perform axial rotations. As bovine tails also lack large tendinous structures (such as those seen in rodent tails), loads applied to each IVD can be attributed primarily to muscles connecting its two flanking vertebrae. Tapering of the tail from proximal to distal ends thus suggests that the proximal discs are subjected to larger muscle loads than the distal discs. Despite similar composition and water content between tail levels (Demers et al., 2004), different tail levels within the same animal are thus likely to have significantly different load histories, resulting in level-wise differences in cell metabolism (Wiseman et al., 2005). Although prior research has uncovered level-wise differences in human (Busscher et al., 2009) and animal (Cunningham et al., 2010; Zimmerman et al., 1992) lumbar discs, and between animal lumbar and caudal discs (Barbir et al., 2010; Elliott and Sarver, 2004; Sarver and Elliott, 2005), the associations between load history, metabolism and mechanics represent a critically underexploited aspect of the bovine caudal model for IVD research.

Despite its low rate of biosynthesis, the IVD is in a constant state of remodeling, with cellular expression of matrix proteins and matrix-degrading enzymes. Numerous studies have confirmed that extracellular matrix homeostasis in the IVD depends on mechanical loading. At the whole-organ level, both immobilization (MacLean et al., 2003; Wang et al., 2018) and overload (Latridis et al., 2011; Walter et al., 2011) induce degenerative responses. Axial compression at moderate magnitude has been shown to maintain healthy matrix biosynthesis (Rosenzweig et al., 2016; Wuertz et al., 2009). At the cellular level, cyclic tensile stretch has been shown to upregulate collagen I, collagen II and aggrecan in the AF (Cho et al., 2011). At lower frequencies, however, the response of AF cells to cyclic stretch shifts from anabolic to catabolic (Gilbert et al., 2010). As internal strains vary throughout the AF under load (Costi et al., 2007), and strain rates may vary during different motions, it is expected that extracellular matrix protein expression will also vary spatially in response to repetitive joint-level loading.

Localized deposition of extracellular matrix material is the mechanism underlying the development of residual strains in various biological tissues (Nelson, 2014). These strains, defined as deformations in the absence of external load, serve to decrease peak tissue stresses when load is applied (Chuong and Fung, 1986). Residual strains have been identified in the IVD, although their role and origin are unclear. While some measurements (Mengoni et al., 2017; Michalek et al., 2012) suggest adaptation of the AF solely to contain NP pressure, others (Duclos and Michalek, 2017) indicate complex strain patterns suggesting response to repetitive bending. Understanding how an IVD has adapted in response to its load history is essential for reducing the risk of injury when loading changes acutely. IVD pathology has been associated with sudden changes in loading, such as following adjacent segment fusion (Hashimoto et al., 2019) or returning from space flight (Johnston et al., 2010). Understanding the biomechanical effects of remodeling would also allow improved clinical recommendations to assist patients with recovery after injury, prior to, or in between, bouts of repetitive activity. This could potentially lead to training recommendations to physically prepare workers for potentially injurious tasks and inform the decisions around rehabilitating injured individuals prior to return to work.

We hypothesized that in vivo loading of bovine caudal IVDs is level dependent, and that exposure to repetitive bending moments would result in localized alterations in AF microstructure and composition. Level-wise in vivo bending moments in the bovine tail were estimated using two methods: analytical modeling of a dynamic lateral bending event and analysis of muscle cross-sections obtained via computed tomography (CT) scans. This parallel approach allowed us to compare estimates of how much bending moment could be applied (from muscle measurements) with estimates of how much bending moment would be necessary for a typical movement (from modeling). Microstructure was assessed by making localized measurements of the collagen fiber crimp period, and composition was assessed via quantitative histology, using Sirius Red and Safranin O/Methyl Blue.

Five proximal (levels c1c2 through c5c6, where cXcY refers to the disc between the Xth and Yth vertebrae counting up distally) bovine tails from skeletally mature animals (Bos taurus, Linnaeus 1758) were obtained from a local abattoir. All tails had skin removed and were frozen in their resting posture at −20°C upon receipt.

Muscle cross-sections were obtained by imaging frozen tails using a clinical CT scanner (LightSpeed VCT, GE Healthcare, Waukesha, WI, USA) with 0.3 mm in-slice pixel size and 0.45 mm spacing between slices. CT scanning was chosen for this study owing to its ability to accurately differentiate between muscle and other tissues, including fat (Engelke et al., 2018). Additionally, three-dimensional radiographs allowed more accurate placement of slices through the mid-planes of the discs than more traditional dissection techniques (Fig. 1A–C). The scans were analyzed using a suite of custom-written MATLAB (MathWorks, Natick, MA, USA) scripts. First a top-down projection was used to manually define lines running through the middle transverse planes of the c1c2 through c5c6 discs (Fig. 1D). Two-dimensional cross-section images were then extracted on these planes, rotated manually such that the sagittal axis was vertical, and stored (Fig. 1E,F). In the two-dimensional section images, IVD, muscle and perimuscular fat are distinguishable by decreasing image intensity. For each tail, the c1c2 slice was manually adjusted to determine upper and lower threshold values for muscle and disc tissue. The same threshold values were then used to produce binary muscle and disc images for each level within the same tail.

Following thresholding, each disc image was analyzed to calculate total area and center of area. Muscle area to either the right or left of the disc center of area was analyzed to calculate total area and center of area relative to the disc center (eccentricity). Maximum contractile force was calculated from the total muscle area and average muscle fiber diameter, as described below (Seow and Ford, 1991). Maximum force was multiplied by eccentricity to yield maximum lateral bending moment. Additionally, forces and moments at each level were normalized to the area or section modulus (ratio of second moment of area to radius), respectively, of the IVD.

Three muscle samples were taken from a bovine tail during dissection, representing lateral and posterior aspects of the proximal tail and lateral aspect of the distal tail. The samples were frozen, and 20 µm cryosections were taken perpendicular to the axis of the muscle. The cryosections were mounted on glass slides and imaged using an inverted phase-contrast microscope (Olympus IMT2, Olympus, Waltham, MA, USA) equipped with a digital camera. Ten muscle fibers from each cryosection were manually measured from the micrographs using ImageJ.

Each vertebral body and surrounding tissue was modeled as a mass, m_i, located at the center of the bone. Masses were assigned proportionally to approximate conical volume (assuming an average outside diameter of 1.58× the diameter of each disc based on CT scans) of each segment and a density of 997 kg m⁻³ (water). Cubic functions were applied to all bone parameters to smooth out variability resulting from hand measurement of the vertebrae. Final model parameters for each level are provided in Table 1.

Eqn 2 was evaluated in MATLAB to determine bending moment at each joint through a full flexion and relaxation cycle with a 0.005 s time step. Additionally, applied force at each level was estimated by dividing bending moment by eccentricity, defined as the distance from the disc center at which it was applied. Based on the CT scans above, eccentricity at each level was defined as the radius of each disc times 2.63. Both peak and root mean square (RMS) values of moment and force were recorded at each level. RMS values were calculated to account for the possibility of non-sinusoidal solutions from sinusoidal motion. In addition to recording absolute values, forces and moments at each level were also normalized to the area or section modulus (ratio of second moment of area to radius), respectively, of the IVD.

Collagen fiber crimp period was measured as previously described (Duclos and Michalek, 2017) on levels c1c2 through c5c6 of five tails. Briefly, 30 µm cryosections were taken in the transverse plane of each disc at approximately mid-height. After mounting on glass slides, the sections were imaged using an inverted microscope equipped with crossed polarizing filters (Olympus IMT2). MATLAB was used to perform fast Fourier transform analysis on user-selected lines of interest parallel to the fiber direction producing a spatial period (Duclos and Michalek, 2017; Michalek et al., 2018). In the inner, middle and outer zones of the anterior, lateral and posterior regions of each disc (Fig. 3), three measurements were made and averaged to calculate the crimp period in each anatomical location.

Distributions of extracellular matrix in the bovine caudal disc at levels c1c2 through c5c6 from four tails were assessed via semi-quantitative histology. Transverse disc cryosections of either 10 µm or 20 µm thickness were placed on coated glass slides (MAS-GP, Matsunami, Osaka, Japan). The 10 µm sections were stained with Methyl Blue (5 mg ml⁻¹ in water; Sigma-Aldrich, St Louis, MO, USA) followed by Safranin O (1 mg ml⁻¹ in water; Sigma-Aldrich). The 20 µm sections were stained with Sirius Red (0.5 mg ml⁻¹; Alfa Aesar, Ward Hill, MA, USA) with aqueous picric acid (60 mg ml⁻¹; Ricca, Arlington, TX, USA) in water. Following staining and washing, the slides were rinsed with xylene, dried and coated with PolyGlass coverslipping medium (Polysciences, Warrington, PA, USA). Both sets of slides were imaged using an inverted microscope equipped with a digital camera. Safranin O/Methyl Blue slides were imaged in brightfield, yielding micrographs with collagen stained blue and proteoglycans stained red (Pedrini-Mille et al., 1992). Sirius Red slides were imaged with crossed polarizing filters, yielding micrographs in which thicker fibers appear red and thinner fibers appear yellow/green (Junqueira et al., 1979). In each imaging mode, care was taken to maintain illumination conditions and camera white balance for all micrographs.

Both sets of micrographs were analyzed using similar methods in MATLAB. First, a user-defined region of interest was used to isolate an area in a particular anatomical location. Within this region of interest, pixels that were either above a threshold value (for Safranin O/Methyl Blue) or below a threshold value (for Sirius Red) were identified as holes in the cryosections. For the remaining pixels, red channel intensity was compared with either blue (for Safranin O/Methyl Blue) or green (for Sirius Red) channel intensity. The number of pixels for which the red channel intensity was higher was divided by the total number of pixels to yield a red fraction within the range of zero to one. The resulting red fractions are representative of relative proteoglycan content in the Safranin O/Methyl Blue slides (Pastoureau et al., 2003) and relative fraction of thick fibers in the Sirius Red slides (Martins et al., 2018; Rich and Whittaker, 2005).

For force and moment estimated from CT, a one-way analysis of variance (ANOVA) tested the effect of joint level. For crimp period and red fraction for both Sirius Red and Safranin O/Methyl Blue staining, a three-way ANOVA was used to test the effects of region, zone and level as variables as well as their interactions. A Tukey–Kramer post hoc test was used to make relevant pairwise comparisons following all ANOVAs. All statistical tests were performed using MATLAB.

The average muscle fiber diameter was 53±12 µm, with no difference between biopsy locations, indicating an isometric contractile force of 150 kN m⁻² at all locations (Seow and Ford, 1991). Total muscle cross-sectional areas ranged from 24.2±1.7 cm² at c1c2 to 9.1±1.1 cm² at c5c6. As there was no laterality observed, right and left side muscles were pooled for lateral bending moment estimates. Estimated maximum contractile forces during lateral bending motions showed a significant (P<0.001) decrease from ∼150 N at the c1c2 level to 50 N at the c5c6 level (Fig. 4). Estimated maximum bending moments trended significantly (P<0.001) downward from ∼3 N m at c1c2 to 0.75 N m at c5c6. Disc area decreased significantly (P=0.002) from 6.3±0.4 cm² at c1c2 to 3.9±0.3 cm² at c5c6.

The analytical model predicted a positive bending moment at all points in the cycle (Fig. 5), suggesting that there is not a braking action at the top of the swing. Given the sinusoidal input to the model, moment versus time was approximately sinusoidal at proximal levels. At more distal levels the bending moment peak was flatter and wider. Non-zero moments at time=0 s are attributed to inertial effects, as this is the point of maximum acceleration. The peak moments occur at time=0.5 s, when acceleration is zero. Peak and RMS bending moments were 4.12 N m and 2.76 N m, respectively, at the proximal-most joint level. Peak and RMS axial forces were 133 N and 88.8 N, respectively, at the proximal-most joint level. At the distal end of the tail, peak and RMS moments decreased to 5.4 N mm and 3.4 N mm, respectively, and peak and RMS force decreased to 0.83 N and 0.52 N, respectively (Fig. 6A,C). Downward trends remained after normalizing bending moment by disc section modulus and force by disc area (Fig. 6B,D).

Sensitivity of the model to elastic and inertial terms was assessed at the sacrocaudal (sc1) and mid-tail (c8c9) joints (Fig. 7). At both levels, bending stiffness of the IVDs and mass of the tail had similar magnitude effects.

In all regions of all levels, there were visible differences in microstructural and compositional micrographs going from outer to inner AF (Fig. 8). In general, the collagen fiber crimp period decreased from the outer AF to the inner AF, as shown previously (Duclos and Michalek, 2017). Sirius Red staining decreased in overall intensity from the outer to inner AF, along with a slight decrease in the amount of green/yellow staining. Safranin O/Methyl Blue staining showed the expected trend from blue (collagen) in the outer AF to red (proteoglycan) in the inner AF.

There were significant decreases in the AF fiber crimp period from inner to outer zones (P<0.001) and from proximal to distal levels (P<0.001), but no effect of region (P=0.490) (Fig. 9). There was a clear downward trend in crimp period only in the outer AF zone in all regions going from proximal to distal levels. Additionally, there was a significant interaction between zone and level (P<0.001) on crimp period. This is evident from a decreasing difference in outer versus inner crimp period moving from proximal level c1c2 to distal level c5c6 (Fig. 4). There were no significant interactions between region and level (P=0.435) or between zone, region and level (P=0.602). Although the effect of region was not significant, post hoc analysis found significant (P<0.05) level-wise differences in the outer AF crimp period only in the anterior and posterior regions.

Relative intensities of red staining for Safranin O/Methyl Blue and Sirius Red are provided in Table 2. ANOVA of Sirius Red found a significant (P<0.001) increase in relative red staining from outer to inner AF zones, indicating increasing thickness of collagen fibers. There was no effect of region (P=0.149) or level (P=0.271), nor any significant interactions. Although not significant at any particular level or region, ANOVA of relative density of Safranin O/Methyl Blue staining, indicating relative proteoglycan content, showed a significant (P<0.001) increase from the outer to inner AF zones. There was no significant effect of level (P=0.722), although inner AF red staining was slightly lower in the anterior and posterior regions of c1c2 than in other levels.

Both the analytical model and analysis of muscle cross-sections estimated similar magnitudes of bending moment in the proximal portion of the tail. Analysis of muscle cross-sections, however, indicated a steeper downward trend moving distally. Bending moments calculated from muscle cross-sections assumed that all muscles to one side or another of the sagittal plane contract during lateral bending. It is possible that some muscles are only activated during dorsal extension movements; however, these muscles are likely to be the ones closest to the sagittal plane and thus have the smallest influence on our calculations. Whereas the muscle study was performed on five tails of similar size, the analytical model was based on bone morphometry from a single, smaller subject. Both approaches neglect the resting posture of the proximal tail (forward flexion as seen in Fig. 1A,D) and the associated vertebral wedging (as seen in Fig. 1D). Although gravity acting on the tail applies both tension (∼12 N) and a flexural moment (∼1.5 N m) to the proximal joints, it is currently unknown how much is borne by the IVDs rather than the interspinous ligaments. Application of both methods to subjects with a larger range of ages and sizes, with consideration for loads and motions (obtained via in vivo tracking) in all three dimensions, is an important avenue of future research. Additional refinement of the model will also require measurement of lateral bending and flexion–extension stiffnesses, which have previously been shown to not be the same (Michalek and Iatridis, 2012), of all joint levels as well as inclusion of joint viscosities.

It has previously been suggested that, based on a swelling pressure of 0.3 MPa, compressive stresses in bovine caudal discs are comparable to those in human lumbar discs (Oshima et al., 1993). In the present study, peak force in the c1c2 disc, normalized by disc area, was estimated to be 0.23 MPa based on muscle cross-section analysis and 0.31 MPa based on the analytical model. These applied stresses correspond (Vergroesen et al., 2014) to NP fluid pressures of 0.5 MPa and 0.67 MPa, respectively; however, normalized RMS forces (Fig. 6B) would result in ∼0.3 MPa fluid pressure. As swelling acts to maintain disc height over time, it stands to reason that swelling pressure would be closer to average (RMS) load-induced pressure than peak. Additionally, these values are the range of intradiscal pressures measured in humans during typical activities (Wilke et al., 1999). In vivo moments in the human lumbar spine during flexion have been reported as high as 6 N m (Rohlmann et al., 2014). Normalizing by section modulus (∼4.8×10⁻⁶ m³ for the human lumbar IVD based on Zhou et al., 2000) yields a value of 1.25×10⁶ N m⁻², which is within the range estimated from CT muscle images (Fig. 4C). It should be noted that CT estimates of force and moment are based on a relationship between muscle fiber diameter and maximum isometric force, whereas in vivo motions are dynamic. The reported values are thus representative of an upper bound (Hill, 1938). Although numerous studies have quantified frequency of bovine tail motions (Dougherty et al., 1995; Miedema et al., 2011), measurements of mammalian tail velocity are less common (Matherne et al., 2018) and muscle shortening velocity is unknown. There is, however, not expected to be a large range in shortening velocities from the proximal to distal tail. Additionally, the analytical model predicts that the peak bending moment occurs at the point of maximum bending, where shortening velocity is zero.

Assembly of collagen fibers into a regular crimp pattern (Diamant et al., 1972) allows the AF tissue to undergo large strains during spinal motions and dictates its nonlinear elastic behavior (Freed and Doehring, 2005). Crimp period in the AF has been shown to correlate with residual strain, i.e. internal deformations in the absence of external load (Chuong and Fung, 1986; Nelson, 2014), with a longer crimp period indicative of a fiber that is resting under tension (Duclos and Michalek, 2017). It is presumed that residual strains arise from non-uniform extracellular matrix deposition in response to load (Cho et al., 2011). The observation that the outer AF crimp period was highest in discs that had been exposed to the highest magnitude repetitive load in vivo confirms this presumption.

The level-wise trend in outer AF fiber crimp period was only significant in the anterior and posterior regions. In human IVDs, fiber tensile strain in the outer AF has been shown to be highest on the lateral aspects during flexion (Heuer et al., 2008). This suggests that the spatial pattern of AF microstructure observed in this study is in response to tissue tensile strains in the AF fiber direction.

It should be noted that the micrographs used to measure fiber crimp in this study were taken in the transverse plane, while the fibers themselves are oriented helically. The fibers were thus not lying flat in the plane of the section, and measured fiber crimp is an underestimate. Fiber helical angle in bovine caudal discs has been shown to be lower in the outer AF than in the inner AF (Michalek, 2019), resulting in a potential underestimate of 9% in the outer AF and 19% in the inner AF. This difference in error is much smaller than the differences in crimp period (nearly twice as high in the outer AF than in the inner AF at proximal levels) observed in this study. Furthermore, these prior measurements showed no effect of tail level on fiber angle at any anatomical location. The level-wise trend observed in the present study is thus unlikely to be the result of an artifact.

Histological estimates of composition are in agreement with previous literature. As expected, proteoglycan content (relative to collagen content) increased from the outer AF to inner AF (Bezci et al., 2019). Previous investigation of proteoglycan content and tail level (Demers et al., 2004) showed no effect with tissue pooled from the entire AF and entire NP. This is confirmed by the present study at all anatomical locations. Birefringence color of collagen fibers stained with Sirius Red is indicative of fiber thickness (Junqueira et al., 1979), which has been associated with both fiber type (Dayan et al., 1989) and maturity (Arun Gopinathan et al., 2015). Whereas green/yellow staining indicate either type II collagen or immature type I collagen, red staining primarily indicates mature type I fibers. Our results show a predominance of red (thick) fibers throughout the AF, which is consistent with previous observations in non-degenerate human IVDs (Martins et al., 2018).

Previous investigation of residual strains in the AF (Duclos and Michalek, 2017) suggested a process in which high stresses in the inner and middle AF stimulate increased collagen I production, resulting in an increase in fiber length. As the fiber endpoints are fixed (Michalek, 2019), this results in residual compression along the length of the fiber and subsequently crimp period shortening and lamellar thickening. Static equilibrium requires tissue adjacent to these areas of residual compression to be in tension, as indicated by fiber crimp period lengthening. An alternative mechanism (Bezci et al., 2015; Yang and O'Connell, 2019) is that residual inner AF compression is the result of proteoglycan-driven swelling. The present results suggest that both of these processes are contributing to the residual strain phenomenon. It is broadly confirmed that, moving from inner to outer AF, the collagen fiber crimp period (which we have previously correlated with residual strain; Duclos and Michalek, 2017) goes up as relative proteoglycan content goes down. However, although not statistically significant, Safranin O staining was reduced in the inner AF in the anterior and posterior regions of the c1c2 disc, adjacent to where the outer AF crimp period was highest. Furthermore, variation in the outer AF crimp period with both level and anatomical region is not accompanied by any similar variations in proteoglycan content.

Bovine proximal caudal IVDs could be a good substitute for cadaveric tissue in studies of the effects of remodeling in the IVD. The predictable variation in disc microstructure would allow analysis of injury and risk patterns that could eventually inform clinical treatment guidelines for primary and secondary injury prevention. Although our observations suggest localized adaptation to repetitive loading, the time scale of this effect is unknown and warrants further investigation.

Our findings indicate that bovine proximal caudal IVDs are subjected to substantial bending moments during frequent lateral bending events, and that the magnitude of bending moment decreases in distal levels. Although the previously reported lack of compositional difference between levels in the bovine tail is confirmed, the differences in microstructure suggest that level-wise effects may be an important consideration when using bovine tail discs in biomechanical studies.

The authors thank Dr Peter Edic of GE Global Research for assistance with CT scanning. Specimens were provided by Tri-Town Packing of Brasher Falls, NY. Additional assistance was provided by the Huntley Family Farm and Casey Havekes of the Cornell Cooperative Extension North Country Regional Ag Team.