Adaptive bone formation in acellular vertebrae of sea bass(Dicentrarchus labrax L.)

Mammalian bone is an active tissue in which osteoblasts and osteoclasts balance bone mass. This process of adaptive modelling and remodelling is probably regulated by strain-sensing osteocytes. Bone of advanced teleosts is acellular yet, despite the lack of osteocytes, it is capable of an adaptive response to physical stimuli. Strenuous exercise is known to induce lordosis. Lordosis is a ventrad curvature of the vertebral column, and the affected vertebrae show an increase in bone formation.

The effects of lordosis on the strain distribution in sea bass(Dicentrarchus labrax L.) vertebrae are assessed using finite element modelling. The response of the local tissue is analyzed spatially and ontogenetically in terms of bone volume.

Lordotic vertebrae show a significantly increased strain energy due to the increased load compared with normal vertebrae when loaded in compression. High strain regions are found in the vertebral centrum and parasagittal ridges. The increase in strain energy is attenuated by a change in architecture due to the increased bone formation. The increased bone formation is seen mainly at the articular surfaces of the vertebrae, although some extra bone is formed in the vertebral centrum.

Regions in which the highest strains are found do not spatially correlate with regions in which the most extensive bone apposition occurs in lordotic vertebrae of sea bass. Mammalian-like strain-regulated bone modelling is probably not the guiding mechanism in adaptive bone modelling of acellular sea bass vertebrae. Chondroidal ossification is found at the articular surfaces where it mediates a rapid adaptive response, potentially attenuating high stresses on the dorsal zygapophyses.

Bone is an active tissue. The activity of osteoblasts and osteoclasts normally balances bone mass in mammalian trabecular bone. Strain-related stimuli are thought to control extra bone formation in high strain regions(Huiskes et al., 2000). Osteocytes are considered to be the strain sensors that recruit new osteoblasts to deposit bone matrix (Burger and Klein-Nulend, 1999). Furthermore, shear and tensile stresses are thought to promote ossification of cartilage(Carter and Wong, 2003).

Bone tissue in the Actinopterygii can be divided into two types(Fig. 1): cellular bone and acellular bone. Cellular bone contains osteocytes in its matrix and is found in the lower taxa of the Actinopterygii. The bone of higher Actinopterygii,especially the neoteleosts, is devoid of osteocytes, as osteoblasts recede from the mineralization front and never become entrapped(Parenti, 1986). Although the phylogenetic dichotomy is clear, exceptions to the rule do occur. Notably within the Protacanthopterygii, the Salmonidae have cellular bone, as do the Thunnini within the Acanthopterygii(Kölliker, 1859; Moss, 1961, 1965). Furthermore, within the Ostariophysi, which normally have cellular bone, larval bone is reported to be acellular (Fleming et al.,2004).

Acellular bone has been regarded as `dead' bone(Moss, 1962), although osteoclasts do invade acellular bone(Glowacki et al., 1986; Sire et al., 1990; Weihs and Watabe, 1979; Witten, 1997). Nevertheless,the absence of canaliculi in acellular bone(Simmons et al., 1970)precludes an adaptive response of bone in the manner of the amniote cellular bone (Currey, 2002).

Despite the lack of strain-sensing osteocytes, acellular bone is capable of showing conformational changes under increased loading(Kihara et al., 2002; Kranenbarg et al., 2005). Strenuous exercise [evoked by, for example, high current velocity(Backiel et al., 1984) or caudal fin removal (Kihara et al.,2002)] is known to induce lordosis in several teleosts. Lordosis is characterized by extra-mineralized vertebrae, accompanied by an abnormal ventrad curvature of the vertebral column(Chatain, 1994; Divanach et al., 1997; Paperna, 1978). Similar deformities have been reported in elasmobranchs(Heupel et al., 1999). The conformation of the affected vertebrae is probably a response of the local tissue to a change in loading (Kranenbarg et al., 2005).

The response of bone to physical stimuli can be classified as a modelling response (growth-related changes in the gross shape of the bone) or a remodelling response (change of architecture rather than total volume of the bone) (Currey, 2002). The mechanisms governing the growth-related modelling process in the acellular vertebrae of lordotic European sea bass (Dicentrarchus labrax L.) are unknown (Kranenbarg et al.,2005).

Adult teleost vertebrae are composed of a vertebral centrum, a dorsal neural arch and a ventral haemal arch. The centra develop through membranous ossification (Bertin, 1958; Liem et al., 2001). The central part of the vertebral centrum is formed by ossification of the perinotochordal sheaths. Rostral and caudal parts are added to the centrum by sklerotomal membranous ossification(Grotmol et al., 2003).

Both the neural and the haemal arch generally develop through cartilage intermediates (see Bird and Mabee,2003), collectively called the arcualia(Liem et al., 2001). The cartilaginous precursor of the arches is maintained for a longer period of time in the lower orders of the teleosts, while ossification proceeds more rapidly in the higher orders(François, 1966). The zygapophyses and bony ridges connecting pre- and postzygapophyses develop through membranous ossification (Mookerjee et al., 1940).

The combination of acellular bone and membranous ossification seems to preclude adaptive modelling in the classical sense for vertebral centra of the neoteleosts.

We present a numerical analysis of the strain energy density distribution in normal and lordotic vertebrae under compression. Subsequently, we present an analysis of the spatial distribution of bone formation in adaptively modelled vertebrae. We discuss the potentials for strain-regulated modelling in acellular bone and provide new insights into the mechanism governing bone modelling in a species without osteocytes to act as strain gauges.

Sea bass (Dicentrarchus labrax L.) eggs were incubated and the larvae were subsequently reared at the Aquaculture institute of the Hellenic Center for Marine Research (HCMR), Crete (for details, see Kranenbarg et al., 2005). Samples for the present study were taken at approximately 35 mm and 45 mm total length (TL) and fixed in PBS buffered formalin. Based on X-radiographs (60 s exposure at 15 kV and 15 mA with an Enraf-Nonius Diffractis 581 X-ray generator; Delft, The Netherlands), four lordotic specimens (34.6-46.6 mm TL) and four non-lordotic specimens(35.8-47.7 mm TL) were selected for the finite element analysis (cf. fig. 1b,c from Kranenbarg et al.,2005). In addition to the eight specimens mentioned, two extra specimens (WAT38-S02, 28.0 mm TL; WAT38-S08, 27.3 mm TL)were selected for the bone volume distribution analysis. One normal(WCT33-S15, TL 46.5 mm) and one lordotic (WCT5-S13, TL 42.9 mm) specimen were selected for an initial histological analysis.

The range of ỹ is given in Table 1.

A large eccentricity causes a large (dimensionless) stress. The radius of gyration, r, can be interpreted as the distance from the neutral axis at which the entire cross-sectional area could be concentrated and still have the same second moment of area as the original area(Gere and Timoshenko, 1999). In case of a large eccentricity, the apposition of extra bone material can theoretically - through an increased r - resist the stress-augmenting effect of the large eccentricity. The dimensionless ratio e/r will be termed the `lordosis ratio' as it signifies the lordosis severity in terms of increased (dimensionless) stress due to the characteristic dorsoventral curvature of the vertebral column.

The eccentricity, e, was measured for each specimen, and the mean radius of gyration in the dorsoventral direction, r, was calculated in Matlab 7.0® (MathWorks, Inc., Natick, MA, USA) for vertebra number 15 of each specimen. The stack of bitmaps was hand-edited before the calculation to remove part of the neural and haemal spines where they did not connect to the vertebral body in cross section. Table 1 gives the parameter values.

The selected sea bass specimens were scanned with a Skyscan® 1072 micro-CT system (Antwerp, Belgium) at 80 kV and 100 μA and with a cubic voxel size of 6.08 μm. Vertebra number 15 (centre of lordosis) was scanned for each specimen. The scans were subsequently reconstructed, segmented and edited according to Kranenbarg et al.(2005) to obtain a stack of bitmaps representing a single vertebra.

The stack of bitmaps from the micro-CT analysis was resampled with a cubic interpolation algorithm to a resolution of 14 μm. Each voxel was converted to a an eight-node hexahedral element (type 7 in MSC.Marc2005; MSC Software Benelux BV, Gouda, The Netherlands) and the resulting surface was smoothed. Mesh generation was performed in Matlab 7.0®. The resulting output file was imported into MSC.Patran2005 for further pre-processing.

The centre of mass (COM) was calculated in MSC.Patran®, and a new coordinate frame was defined with the origin in the COM and the axes along the principal axes of the vertebra. Subsequently, the dimensionless stress distribution,

⁠, was applied perpendicular to the free face of each element on the frontal surface of the vertebra (rim of the vertebral centrum and prezygapophyses). This procedure ensures the component of the compressive force parallel to the long axis of the vertebra to be equal to F, as defined in the previous section. It was not feasible to incorporate full biological complexity (e.g. exact positions of tendons and muscle fibres) in our current model and therefore our mode of application of the loads is necessarily simplified. Nevertheless, Saint-Venant's principle attenuates the effects of possible inaccuracies in our mode of application by stating that the stress distribution may be assumed independent of the actual mode of application of the loads (except in the immediate vicinity of the points of their application) (Beer and Johnston,1992).

As the range of ỹ is comparable for all vertebrae in the present study (see Table 1), the applied dimensionless stress is determined primarily by the lordosis ratio e/r. The mean lordosis ratio of the normal vertebrae is 2.29 and of the lordotic vertebrae is 4.17 (see Table 1). Fig. 3 shows the applied dimensionless stress distributions for a normal (blue line) and a lordotic(red line) vertebra.

As mammalian bones respond to mechanical loads by altering their architecture, rather than the quality of the bone tissue(Currey, 2002), both the normal and lordotic vertebrae were assumed to consist of the same material. The vertebra was specified as an isotropic linear elastic material with a tissue-level Young's modulus (E) of 1 GPa. As young bone is characterized by a smaller E than old bone (due to a lesser degree of mineralization; Currey, 2002),a value was chosen from the lower range of values reported for the Eof bony tissue (Cowin, 2001). The apparent E of three juvenile teleost species' whole vertebra is reported to be approximately 75 MPa(Hamilton et al., 1981). Back-calculation with FE-analysis of the tissue-level E from the apparent experimental E yielded a tissue-level Eapproximately 15 times the apparent experimental E(Hou et al., 1998) for human vertebral body specimens. This factor supports our value of 1 GPa for juvenile piscine bone. The Poisson ratio was chosen as 0.3(Cowin, 2001).

To analyze the effect of increased loading on normal vertebrae, we doubled the eccentricity for the normal vertebrae, thereby simulating a lordotic configuration for the normal vertebrae.

where u_t is the total strain energy density and u_d is the strain energy density due to direct compression. The results were post-processed in MSC.Patran.

No lordosis was observed in any specimen under 30 mm TL. The distribution of bone volume over the anterior-posterior axis of each vertebral centrum was resampled to 300 data points to obtain the bone volume as a function of position (%) along the vertebral centrum. Subsequently, the bone volume at each position was regressed on total length of the animal with a quadratic function. This procedure was performed for lordotic and non-lordotic animals separately, yielding two surfaces showing differences in bone volume over the length of the vertebral centrum. A 95% confidence interval was calculated for both surfaces. All calculations were performed in Matlab 7.0®.

The specimens were post-fixed for 1 week in Bouin's fixative(Romeis, 1968) (modified to 17 vol. % glacial acetic acid). The specimens were embedded in paraffin,sectioned parasagittally at 5 μm and stained according to Crossmon's protocol (Romeis, 1968).

A lordotic vertebral column is characterized by an eccentricity that is about twice as large as that of a normal vertebral column(Table 1). Fig. 4A,B shows the effect of a doubled eccentricity without an accompanying change in architecture. The mean ũ (which equals the mean total dimensionless strain energy of the vertebra, assuming an equal volume of all finite elements) is significantly increased in normal vertebrae in such a lordotic configuration(Fig. 4A). This is further illustrated in Fig. 4B, which shows a significantly larger percentage of elements with large ũ (>2000) in the normal vertebra with double eccentricity. The opposite is true for the percentage of elements with a small ũ (Fig. 4B).

Fig. 5 shows the distribution of ũ over the centra of the vertebrae (vertebral body). Normal vertebrae of both the 35 mm TL stage and the 45 mm TL stage show only limited regions of high ũ(Fig. 5A,D). The same vertebrae show extensive regions of high ũ when the eccentricity is doubled (Fig. 5B,E). Both stages show high ũ values in the struts connecting the neural arch with the parasagittal ridges, and the parasagittal ridges themselves. In addition, the 35 mm TL stage shows high ũ values in the dorsal and ventral part of the vertebral centrum sensu stricto.

A comparison of the average ũ between normal vertebrae(with normal eccentricity) and lordotic vertebrae (with a characteristically large eccentricity) shows a significantly increased value in the lordotic vertebrae (Fig. 4C), although the difference is less pronounced than in Fig. 4A. This is illustrated again in Fig. 4D, which shows a significantly larger proportion of elements in the lordotic vertebrae with a high ũ, although the proportion of elements with a low ũ only tends to be smaller in lordotic than in normal vertebrae.

The distribution of ũ over the vertebrae shows clear differences in magnitude between normal and lordotic vertebrae(Fig. 5A,C,D,F), with the parasagittal ridges being the high ũ regions. These quantitative differences between normal and lordotic vertebrae are present at both the 35 mm TL stage and the 45 mm TL stage.

By comparing the normal vertebrae with double eccentricity with the lordotic vertebrae, we analyzed the effect of architectural changes on the ũ distribution. The eccentricity is similar in both groups,only the architecture of the vertebrae differs(Kranenbarg et al., 2005). A comparison of average ũ shows no significant differences,although the ũ in lordotic vertebrae tends to be smaller(Fig. 4E). The same trend is observed in the histogram analysis of Fig. 4F.

The distribution of ũ does show clear differences between the lordotic vertebrae and normal vertebrae with double eccentricity,especially in the 35 mm TL stage(Fig. 5B,C,E,F). Although the parasagittal ridges show high ũ values in both groups, the regions of high ũ in the vertebral centrum sensu stricto in the normal vertebrae with double eccentricity(Fig. 5B,E) have disappeared in the lordotic vertebrae (Fig. 5C,F).

The absolute value of the ũ distribution reaches peak values of well over 2000 and is nearly 200 on average in the normal vertebrae with normal eccentricity. This reveals that even the relatively small eccentricity of a normal vertebral column(Table 1) causes strain energy densities many times larger than those in a compressed, perfectly straight column.

Fig. 6A shows the distribution of bone over the length of the vertebral centrum in lordotic(red) and normal (blue) vertebrae during growth. The amount of bone increases with TL in both groups. The difference in bone volume between the normal and lordotic vertebrae is shown in Fig. 6B. The dark green region indicates where the 95% confidence intervals of the surfaces of Fig. 6A do not overlap. Although some additional bone is deposited in the central part of lordotic vertebrae, most of the extra bone material in lordotic vertebrae is deposited at about 20% and 80% vertebral centrum length, which are the articular areas of the vertebrae.

Fig. 7 shows a parasagittal section through a normal (A) and a lordotic vertebra (B). The intervertebral ligament connects the two consecutive normal vertebral centra. These vertebral centra, as well as the projecting zygapophyses, consist of acellular bone(Fig. 7A). The intervertebral ligament is similarly present in the lordotic vertebra, and the vertebral centrum also consists of acellular bone(Fig. 7B). The zygapophyses of the lordotic vertebrae, however, are highly proliferated and the tissue in between the postzygapophysis of one vertebra and the prezygapophysis of the next ranges from densely fibrous to chondroid tissue. Chondroid tissue was not observed in the normal specimen.

We performed a dimensionless stress-strain analysis on normal and lordotic vertebrae of sea bass juveniles. A dimensionless analysis had two important advantages over an analysis based on actual stress and strain values. First,the analysis was rendered dimensionless by scaling to the direct compressive stress component σ_d. This made the analysis independent of the magnitude of the direct compressive stress, which proved difficult to estimate reliably (although see later in the Discussion). Second, the dimensionless analysis prompted the introduction of the lordosis ratio (see Materials and methods). The lordosis ratio enabled us to define a typical normal vertebra (lordosis ratio 2.29) and a typical lordotic vertebra(lordosis ratio 4.17), including their respective loading patterns(Fig. 3), irrespective of body size. The lordosis ratio is closely related to the eccentricity ratio(Gere and Timoshenko,1999).

Lordosis can be characterized as a ventrad curvature (double eccentricity when compared with the normal situation; see Table 1) of the vertebral column with an accompanying local response of the tissue. The effect of the ventrad curvature is an extensive increase in strain energy stored in the vertebrae (due to the excessive induced bending moment), if the architecture of the vertebrae would not change (Fig. 5A,B,D,E).

The local response of the tissue to the increased loading regime effectively removes the high strain regions from the vertebral centrum sensu stricto, possibly preventing the strains from exceeding the yield point. High strain regions remain, however, in the parasagittal ridges(Fig. 5C,F). This effect is especially clear in the 35 mm TL stage, as in the older stage, the normal development of the vertebrae already accomplishes a similar effect of redirecting the high strain regions to the parasagittal ridges.

Fig. 6 shows that the articular surface areas are the main sites of bone apposition in lordotic vertebrae. Kranenbarg et al.(2005) and Fig. 5 show that the apposition results in large articular areas in lordotic vertebrae. This effectively reduces the pressure on the zygapophyses.

Although a similar phenomenon is seen in dogs, where the articular surfaces of caudal vertebrae are disproportionately larger in large breed dogs(Breit, 2002), articular surface areas are generally found to be not phenotypically plastic in mammals on an individual basis (Lieberman et al.,2001).

Fig. 6 also shows a limited amount of extra bone being formed in the central region of lordotic vertebrae. Fig. 5(arguably) shows that in lordotic vertebrae more struts are formed that connect the neural arch with the dorsal postzygapophysis. The extra bone in the central region does increase the second moment of area in this region(Kranenbarg et al., 2005) and thus helps strengthen the entire vertebra. Nevertheless, the extra struts cannot prevent regions of high strain occurring in the parasagittal ridges of lordotic vertebrae, and the mean ũ remains significantly larger than in normal vertebrae with normal eccentricity.

The adaptive modelling of the articular surfaces of lordotic vertebrae is only partially effective in counteracting the overall strain-increasing effect of the ventrad curvature. This incomplete compensatory response is further illustrated by the larger value of the lordosis ratio in lordotic vertebrae(Table 1), which signifies the lordosis severity, as explained in the Materials and methods.

Adaptive modelling of bone is well documented in mammals(Cowin, 2001; Currey, 2002), yet (acellular)teleost bone is also capable of an adaptive response(Glowacki et al., 1986; Witten and Villwock, 1997; Huysseune, 2000). Strain-related stimuli are known to induce adaptive (re)modelling in bone(Cowin, 2001; Currey, 2002). In mammals,osteocytes are considered strain gauges that recruit osteoblasts to form bone(Burger and Klein-Nulend,1999).

Frost (1983) introduced the`minimum effective strain (MES)' principle, stating that `strains smaller than the MES would not evoke bone architectural adaptations, but those larger than the MES would'. The MES was found to be around 1000μϵ(Cowin, 2001). It should be noted that the MES principle is under much debate(Currey, 2002) and has been modified over the years. Lanyon(1987) proposed to generalize the principle to the `minimum effective strain-related stimulus (MESS)' to include strain distribution and strain rate in the stimulus.

Despite the lack of consensus, a generalized comparison of the MES with our results is useful. Assuming linear elasticity in our model and a direct compressive stress component of ∼22 kPa(Cheng et al., 1998), we can recalculate the MES to a minimum effective ũ of approximately 2000. Fig. 5 shows that,especially in the normal vertebrae with doubled eccentricity(Fig. 5B,E) and the lordotic vertebrae (Fig. 5C,F),extensive regions occur with a ũ higher than the minimum effective ũ (red regions). As the value of 22 kPa for the direct compressive stress component is not experimentally tested, we note that, although variations in the actual value may either reduce or expand the potential modelling area, our conclusion that extensive regions with a ũ higher than the minimum effective ũ occur remains valid.

For a stimulus to be effective, however, a sensor should be present. As sea bass bone is acellular, no osteocytes are present and the bone matrix lacks a strain gauge as is present in mammalian bone. A comparison of Figs 5 and 6 reveals that regions in which the highest strain energy densities are found (especially the parasagittal ridges) do not correlate spatially with the regions in which the highest amount of extra bone is deposited (especially the articular surfaces). A strain (related) stimulus is therefore unlikely to be the sole explanatory variable for adaptive modelling in the sea bass vertebral articulations.

The articular surfaces show extensive formation of extra bone, although the zygapophyses develop as extensions from the acellular vertebral centrum. Histology reveals that this articular surface modelling is mediated by chondroid tissue (cf. Beresford,1981). Chondroid tissue is also found in articular surfaces of other joints (Beresford, 1981)and is hypothesized to be especially suited to fulfil a fast reconstruction demand (Huysseune and Verraes,1986; Huysseune et al.,1986; Huysseune,2000). Apparently, in the central part of the vertebrae,conditions are not met to develop chondroid, while compression of the intervertebral disk causes chondroid tissue to develop. The chondroidal ossification subsequently results in a fast adaptive modelling of the articular surfaces.

Experiments on mammals have shown that bone material deposited in a modelling response is very similar to the original bone(Currey, 2002). We assumed all bone material detected by the CT scanner had equal material properties. The actual mechanical properties of chondroid bone are unknown, although Meunier and Huysseune (1992) suggest mechanical excellence to be inferior to speed of reconstruction. Further testing is needed to assess the material properties of lordotic sea bass vertebrae.

Despite the lack of osteocytes, a limited modelling response is found in the high-strained central part of lordotic sea bass vertebrae. Osteoblasts and fibroblasts are also known to be strain sensitive(Buckley et al., 1988; Chiquet et al., 2003) and might induce the extra bone deposition that does occur in the central part of lordotic vertebrae.

While the central part of lordotic vertebrae shows possibly strain-mediated modelling, the volumetrically most important modelling response at the articular surfaces does not correlate spatially with a high strain energy density. Chondroid tissue appears to meet the demand for an accelerated growth rate and the demand for a shear-resistant support(Huysseune, 2000; Meunier and Huysseune, 1992). Cartilage growth and ossification is known to be promoted by shear stresses(Carter and Wong, 2003). Cartilage mechanobiology might therefore provide interesting new insights into the mechanisms of chondroidal ossification and modelling in acellular bone.

Our dimensionless analysis enabled us to define the load on a `typically normal' and a `typically lordotic' vertebra, as illustrated in Fig. 3. Furthermore, the analysis prompted the introduction of a lordosis ratio as a dimensionless factor describing the severity of the lordosis irrespective of body size.

The present paper further shows that, although some extra bone is formed in the central part, the articular surface areas are the main modelling sites in lordotic vertebrae. The highest strain energy densities are found in the central part of the vertebrae. Mammalian-like strain-regulated bone modelling is therefore unlikely to be the sole mechanism guiding sea bass bone modelling. The adaptive response at the articular surfaces is mediated by chondroidal ossification, possibly induced by other physical stimuli (shear stresses). The increased articular surface area will effectively reduce the pressure on the dorsal zygapophyses.

Despite the adaptations, the overall strain energy of lordotic vertebrae remains higher than that of normal vertebrae in our numerical analyses. Apparently, the morphological adaptations cannot fully compensate for the effect of an increased eccentricity in terms of strain energy.

List of symbols and abbreviations
 
• A

  cross-sectional surface area

 
• COM

  centre of mass

 
• e

  eccentricity

 
• e/r

  lordosis ratio

 
• E

  Young's modulus

 
• F

  compressive force

 
• M

  induced bending moment

 
• MES

  minimum effective strain

 
• I

  second moment of area

 
• r

  radius of gyration

 
• TL

  total length

 
• u_d

  strain energy density due to direct compressive stress

 
• u_t

  total strain energy density

 
• ũ

  dimensionless strain energy density

 
• y

  spatial variable

 
• ỹ

  dimensionless spatial variable

 
• μϵ

  microstrain (equals strain × 10⁶)

 
• σd

  direct compressive stress

 
• σi

  induced stress

 
• σt

  total stress

 
• \({\tilde{{\sigma}}}\)

  dimensionless stress

The authors gratefully acknowledge the valuable comments of David Lentink.