Skip to main content

An adaptation model for trabecular bone at different mechanical levels

Abstract

Background

Bone has the ability to adapt to mechanical usage or other biophysical stimuli in terms of its mass and architecture, indicating that a certain mechanism exists for monitoring mechanical usage and controlling the bone's adaptation behaviors. There are four zones describing different bone adaptation behaviors: the disuse, adaptation, overload, and pathologic overload zones. In different zones, the changes of bone mass, as calculated by the difference between the amount of bone formed and what is resorbed, should be different.

Methods

An adaptation model for the trabecular bone at different mechanical levels was presented in this study based on a number of experimental observations and numerical algorithms in the literature. In the proposed model, the amount of bone formation and the probability of bone remodeling activation were proposed in accordance with the mechanical levels. Seven numerical simulation cases under different mechanical conditions were analyzed as examples by incorporating the adaptation model presented in this paper with the finite element method.

Results

The proposed bone adaptation model describes the well-known bone adaptation behaviors in different zones. The bone mass and architecture of the bone tissue within the adaptation zone almost remained unchanged. Although the probability of osteoclastic activation is enhanced in the overload zone, the potential of osteoblasts to form bones compensate for the osteoclastic resorption, eventually strengthening the bones. In the disuse zone, the disuse-mode remodeling removes bone tissue in disuse zone.

Conclusions

The study seeks to provide better understanding of the relationships between bone morphology and the mechanical, as well as biological environments. Furthermore, this paper provides a computational model and methodology for the numerical simulation of changes of bone structural morphology that are caused by changes of mechanical and biological environments.

Background

Bone is a living organ; it has the ability to adapt to mechanical usage or other biophysical stimuli in terms of its mass and architecture. This attribute is known as functional adaptation [1, 2]. Modeling by drifts and remodeling by groups of osteoclasts and osteoblasts that are organized into basic multicellular units (BMUs) determine both the bone mass and architecture, with the exception of longitudinal bone growth. Bone modeling works best during the growing years and works poorly on adult cortical bone. However, it works satisfactorily within trabeculae throughout life [3]. Bone remodeling consists of biologically coupled BMU activation, bone resorption by osteoclasts, and bone formation by osteoblasts. It occurs in all in vivo bone tissues and is an important way to renew bone.

The observed adaptations of bone mass and architecture to the mechanical usage show that a certain mechanism exists for monitoring the mechanical usage and controlling bone adaptation behaviors. A conceptual "mechanostat" theory developed in a previous study indicates that the mechanically generated bone strain signals have important roles in governing bone adaptation behaviors [4, 5]. Threshold ranges of such signals appear to reside in certain skeletal cells as genetically determined internal standards. Minimum effective strain (MES) has been used to describe the above internal standards (or "thresholds") [4]. Three thresholds, i.e., remodeling, modeling and microdamage thresholds (MESr, MESm and MESp, respectively) serve as boundaries of four zones describing different bone adaptation behaviors: the disuse, adaptive, overload, and pathologic overload zones.

Bone adaptation process is controlled by mechanical usage and biological factors, which have coupled contributions [6]. Aging, menopause, drug treatments, and so on, are known as biological factors. Researchers found that biological factors regulate the rate of bone turnover; however, the overall balance between bone formation and resorption is influenced by prevailing levels of mechanical usage [7].

Nowadays, heavy computing power has been mastered by human beings, making the quantitative simulation of bone adaptation process possible; this is driven by the goals of further predicting and explaining the bone adaptation behaviors, as well as the formation and maintenance of bone architecture under different mechanical and estrogen levels. In doing so, the method of conducting actual experiments has limited efficacy since certain difficulties existed in some experimental investigations. Simulation saves not only time, but also expenditures spent on research and development, it is also regarded as 'the third way of science' [8]. Many attempts have been made to gain quantitative insight into the bone modeling and remodeling processes:

Researchers begin with theoretical analyses before conducting numerical simulation on bone structures in cooperation with the finite element method. One of the fundamental theories ("theory of adaptive elasticity") for cortical bone was based on general continuum mechanical principles [9]. Based on this theory, a finite element model has been developed [10]. A different approach to predict bone adaptation behavior postulated that bone is a self-optimizing material able to adapt its orientation and density in response to its stress/strain state [11]. In another remodeling algorithm, the idea of a lazy or homeostatic zone has been included [12]. The essential idea is that bone mass increases above a certain level of strain or strain energy density; in addition excessive bone remodeling can be observed below a certain threshold. In between the two levels, the bone structure is maintained, representing the case when bone is under physiological loading conditions encountered in normal activities [13]. The idea of homeostatic zone is similar with the adaptation zone proposed by Frost [4, 5]. There are also a number of alternative remodeling algorithms that have been proposed and used in previous study [14–16]. All trabecular models simulated the outcome of coordinated osteoclastic and osteoblastic activities as either a net increase or decrease in density. A model has been developed that included separate descriptions of osteoclastic resorption and osteoblastic formation, enabling simulation of trabecular bone growth, adaptation and maintenance [17, 18]. It was assumed that osteocytes can transfer an osteoblast recruitment stimulus to the surface and then enhance bone formation; in addition, osteoclasts could resorb bone that is in disuse state or has been damaged in a spatially random manner [17]. In addition, using six bone resorption models, Tanck et al. investigated how using alternative strain-based local stimuli for osteoclasts in bone resorption would affect remodeling and adaptation of trabecular architecture [19]. Vahdati and Rouhi added the cellular accommodation effect into the model and took into consideration both microdamage and disuse on activation of resorption [20]. Nowadays, increased computational resources have made it possible for large-scale bone-remodeling simulations to have an element size that is as small as about 50 με. This demonstrates that bone remodeling at the tissue level can create a highly complex and optimized trabecular structure in terms of bone density and orientation [21–23].

In the current quantitative bone adaptation models, the amount of osteoblastic bone formation is proportional to the difference between the mechanical stimulus and a certain threshold (setpoint). According to Frost's mechanostat theory, the changes of bone mass in different zones, as calculated by the difference between the amount of bone formation and that of bone resorption, should be different, resulting in different bone adaptation behaviors. This means that the coupling relationship between bone formation and resorption may be different. In this paper, an adaptation model for trabecular bone at different mechanical levels is proposed based on a number of experimental observations and numerical algorithms in the literature, in which the amount of bone formation and the activation of bone adaptation behaviors have been proposed in accordance with the mechanical levels. It is convenient to determine which zone a bone surface lies in as well as calculate the changes of its mass. Hence, more insights about the relationships between bone morphology and the mechanical/biological environments may be gained.

Methods

This section presents an adaptation model for trabecular bone at different mechanical levels. Seven numerical simulation cases under different mechanical conditions are then analyzed as examples by incorporating our adaptation model with the finite element method on a simplified two-dimensional finite element model of a trabecular bone.

Effects of strain energy density (SED) on osteoblastic formation and the probability of osteoclastic activation

Osteocytes act as sensors of mechanoreceptors and regulators of bone mass by mediating the actor cells (osteoblasts and osteoclasts) [24, 25]. In mathematical models governing the bone adaptation process, the strain energy density has been suggested as the mechanical signal that the osteocytes appraise [12]. The total mechanical stimulus at location x on the trabecular surface (P(x,t)) is the contribution of all the osteocytes made in the forms of SED, relative to their distance from x [26]:

(1)

where μ i is the mechanosensitivity of the osteocyte i, R i (t) is the SED of the osteocyte i, and f i (x) is the spatial influence function, which describes the influence of osteocyte i on the osteoblasts and osteoclasts at location x [15]:

(2)

where d i (x) is the distance between osteocyte i and location x. The parameter D represents the distance from an osteocyte, at which location the effect has been reduced to 0.36788.

There are about 1600 osteocytes per square millimeter. An osteocyte is located in the center of each bone-filled element [27]. The influence parameter D was equal to 100 μm [15].

For the osteoclastic activation, experimental and numerical studies show that apart from disuse-activated resorption, microdamage-stimulated resorption also correlates with the SED value in bone tissue, i.e., accumulated osteocyte signal [20, 28].

Based on the aforementioned investigations, several SED thresholds were used to represent different mechanical levels for osteoblast formation, as well as the probability of osteoclast activation at the tissue level. Then, the adaptation model for the trabecular bone at different mechanical levels was presented, resulting in different adaptation behaviors at the apparent level as described by Frost's mechanostat theory [4, 5]. The symbols and expressions of the SED thresholds are listed in Table 1.

Table 1 Symbols and expressions of the SED mechanical thresholds, and the values in our numerical simulations

Osteoclast activation and bone resorption

The amount of bone resorbed by osteoclasts at a trabecular surface patch is denoted as r oc , where the affected surface patches is determined by chance (Huiskes et al. [17]). These authors introduced a probability of osteoclast activation per surface site at any time, which is assumed to be regulated either by the presence of microcracks or by disuse. Vahdati and Rouhi [20] proposed that the probability of microdamage-stimulated resorption is a quadratic function of the accumulated osteocyte signal, and that below a critical load value, microdamage would not initiate remodeling, thus include a region similar to the lazy zone. Based on the aforementioned investigations, the probabilities of osteoclast activation per surface site at different mechanical levels, i.e., in different mechanical zones are described below.

  1. (1)

    In the disuse zone, the probability of osteoclast activation was regulated by disuse. The response of the probability of the osteoclast activation to the osteocyte signals is set to be sigmoidal, similar to the responses found in pharmacological applications [29].

  2. (2)

    Under physiological loading conditions encountered in normal activities, the probability of osteoclast activation is equal for all surface sites. This assumption is based on the observations that remodeling continuously renews the skeleton daily, and that the dynamic forces of daily living produces microcracks. The distribution of microcracks has been found to be spatially random [17, 30].

  3. (3)

    When overloaded, the probability of osteoclast activation is assumed to be regulated by microdamage as the first phase of remodeling to repair damaged regions. Thus, the quadratic function of the accumulated osteocyte signal, as suggested by Vahdati and Rouhi [20], is used to describe the probability of osteoclast activation in this zone.

Thus the resorption probability (p(x,t)) as a function of stimulus value in the form of SED (P(x,t)) can be expressed as:

(3)

where p max, a, b, c, d, and g are constants. In the numerical examples presented in this paper, the values of these constants are p max = 0.2, a= 90.0, b= 0.0714, c= 0.0275, d= 176.0, and g = 19.506.

These parameter values have been chosen as follows: since c represents the resorption probability in the adaptive zone, it has been given a small value of 0.0275. For the resorption probability in the disuse zone, the form is sigmoidal and the maximum value is p max; it has been set to 0.2 [20]. Parameters a and b are coefficients that define the slope and inflection point of the curve, respectively. The inflection point has been chosen as K OB , and a is determined from the condition of continuity at P(x,t) = K AD 1. The resorption probability is a quadratic function when overloaded. By assuming the resorption probability of 0.25 when P(x,t) = 0.122436 MPa and that the condition of continuity is at P(x,t) = K AD 2, the values of d and g can thus be obtained. As a result, the resorption probability is 1 when P(x,t) = 0.149 MPa. For a large range of mechanical stimulus from 0 to 0.149 MPa, the average resorption probability was 0.2148, which is of the order of trabecular bones that are regenerated per year [17]. The plot of resorption probability as a function of stimulus value is shown in Figure 1.

Figure 1
figure 1

Relationship between resorption probability and mechanical stimulus in the numerical simulation.

When the mechanical stimulus remains lower than K RD , i.e., the threshold for resorption drift of bone modeling, the bone surface is subject to resorption by osteoclasts alone, which is the case for bone resorption drift in the bone modeling process.

Osteoblast recruitment and bone formation

Osteoblast recruitment at a surface site is determined by the formation stimulus in the form of SED. The proposed amount of bone formation per surface site (r ob ) at different mechanical levels is expressed below:

  1. (1)

    If the total stimulus is less than the osteoblast formation threshold, bone formation will not occur.

    (4-1)
  2. (2)

    If the stimulus is higher than the osteoblast formation threshold (K OB ), but is lower than the lower threshold for the adaptive zone (K AD1 ), bone formation will occur at the surface site. This is proportional to:

    (4-2)

where Ï„ 1 is constant.

  1. (3)

    When the stimulus is between the two thresholds for the adaptive zone, the amount of bone formed per surface site (r ob ) will be equal to the amount of bone resorbed (r oc ):

    (4-3)

resulting in no net change of bone mass in the adaptive zone.

For Equations (4-2) and (4-3), under the condition of continuity, r ob = r oc when P(x,t) = K AD 1, thus .

  1. (4)

    If the stimulus is higher than the higher threshold for adaptive zone (K AD 2), bone formation at the surface site will be proportional to:

    (4-4)

where Ï„ 2 is a proportionality constant that had been chosen as Ï„ 2 = 1 MPa-1 in the numerical simulation presented in this paper [31].

Figure 2 shows the schematic diagram of the relationship between the amount of bone formation (r ob ) and the mechanical stimulus P(x,t).

Figure 2
figure 2

Schematic diagram of the relationship between the amount of bone formation ( r ob ) and the mechanical stimulus P ( x,t ).

The local change of the relative bone density m (1.0 for fully mineralized tissue) in location x can be expressed as:

(5)

The Young's modulus E(x,t) at each location is density-dependent:

(6)

where E max was the maximum tissue level Young's modulus: E max = 500 MPa [32], and r = 3 [33].

Numerical approach - a simplified two-dimensional finite element model of a trabecular bone

The bone adaptation simulation was performed on a simplified two-dimensional finite element model of a 2 mm × 2 mm portion of bone volume with a thickness of 0.02 mm; it was loaded at the edges by a tensile and compressive principal stress with a uniform density distribution of m(x,t 0) = 0.7, in which t 0 represented the beginning of the simulation (Figure 3). The orientation of the principal stress was φ = 30°. The bone plate model was meshed with 80 × 80 four-node bilinear elements. Boundary conditions were considered important since only a small volume of bone was to be analyzed. Side plates with 2-element thicknesses were added to the bone volume to create even deflection at the edges. The side plates were of constant maximal density and did not remodel. As suggested in [31], surfaces close to the edge of the mesh, particularly surface corners, have a reduced osteocytic environment to receive stimuli compared with the elements located in the centre of the mesh. A correction of stimulus was made by a factor in the regions close to the edge of the mesh, as shown in [31]. The structure was loaded by a sinusoidal stress cycling between 0 and 2 MPa at 1 Hz; it can be calculated using the static finite element analysis [31]. The static finite element analysis combined with the bone remodeling equations proposed by Zhu et al. [16] was performed using the ANSYS software to obtain the homeostatic architecture of the bone model. The volume fraction of this homeostatic architecture was 0.613, similar to that obtained in [31]. Mechanical strain in the bone tissue was about 1041 με, which was within the adaptive zone defined by the mechanostat theory of Frost [3, 4]. Hence, the reference value for mechanical stimulus in this simulation was chosen as the average of the two SED thresholds defining the adaptive zone, i.e., K AD 1and K AD 2, in the numerical simulations that follow. The values for the SED thresholds in Table 1 are chosen so that the numerical results in this paper may agree with the experimental observations.

Figure 3
figure 3

The initial finite element model. The 2 mm × 2 mm portion of bone volume was meshed with 80 × 80 four-node bilinear elements, representing the bone adaptation area; side plates with 2-element thicknesses were added to the bone volume to create even deflection at the edges, which did not remodel.

The homeostatic architecture served as the starting point for the following simulation cases. The computer simulation was conducted as an iterative process in computer time, in which the relative density m per element was regulated between 0.01 and 1.0 (Eq. (5)). The bone adaptation scheme for the 2 D plate model was performed using the ANSYS Parametric Design Language (APDL). Each surface site was given a stochastic number by random function in the APDL for every iterative time step. If the stochastic number was smaller than the resorption probability p(x,t) in Eq. (3) for the same surface site, this site will be subjected to osteoclastic resorption. Once the osteoclasts were activated in a surface site, they resorbed a fixed amount of bone tissue with a relative density of 0.38, which was derived from the following aspects. The area within an osteonal cement line was 2.84 × 10-2 mm2 [34]. However, in the trabecular bone, BMUs resorbed and refilled trenches rather than tunnels; hence, the bone remodeling area was assumed to be 1.42 × 10-2 mm2. The resorption period was about 60 days [35]. The area resorbed by osteoclasts was 2.367 × 10-4 mm2/day. Each finite element area in the model was 6.25 × 10-4 mm2, indicating that a relative density of 0.38 was resorbed by osteoclasts.

The trabecular architecture in the following cases was simulated by incorporating the proposed bone adaptation model with the finite element method. This was done to verify whether the proposed bone adaptation simulation method is capable of simulating the adaptation behavior of the trabecular bone at different mechanical levels.

Simulation Cases

Simulation 1: Disuse

The mechanical loading condition of disuse was applied to the homeostatic architecture within a period of eight years (four years before menopause and four years after menopause). The mechanical loading on the homeostatic model decreased according to the curve in Figure 4 to simulate mechanical disuse [2]. The percentage change in mechanical loading was chosen so that the predicted change in bone mineral density (BMD) will correspond with the clinical data that measures the lumbar spine BMD every six months throughout the menopausal transition by Recker et al. [36].

Figure 4
figure 4

Pattern of disuse. The mechanical loadings decreased according to the curve to simulate mechanical disuse. The values on the curve show the percentages in the mechanical loadings to their magnitudes for the homeostatic structure. B describes the percentage of mechanical loadings in the end of menopause.

Figure 5(a) shows the homeostatic structure of the bone model, which is the initial architecture of simulations 1, 2, 4, 5, 6 and 7.

Figure 5
figure 5

The homeostatic architecture of the bone model, the initial architectures and loading conditions of the seven simulation cases, and their simulation results. a The homeostatic structure of the bone model, which is the initial architecture of simulations 1, 2, 4, 5, 6 and 7. b The initial architecture of simulation 3: two struts were artificially disconnected (inside the two circles), while the external stress was maintained. c The initial architecture and loading condition of simulation 4. The initial structure was the homeostatic structure. The orientation of the applied stress was changed from 30° to 0°. d Result of simulation 1. e Result of simulation 2. f Result of simulation 3. g Result of simulation 4. h Result of simulation 5. i Result of simulation 6. j Result of simulation 7.

Simulation 2: Overloading

A 20% increase in the magnitude of the mechanical loading was imposed on the homeostatic structure to investigate the effect of overloading.

Simulation 3: Artificially disconnected trabeculae

Two struts in the homeostatic architecture were artificially disconnected, while the same externally applied loads were maintained (Figure 5(b)), as proposed by Vahdati and Rouhi [20].

Simulation 4: Rotation of the external load

The orientation of the applied stress on the homeostatic architecture was changed from 30° to 0° (Figure 5(c)).

Simulation 5: Some increase in the external load within physiological loading condition

In this simulation, an 8% increase in the magnitude of the mechanical loading was imposed on the homeostatic structure.

Simulation 6: Some decrease in the external load within physiological loading condition

Here, an 8% decrease in the magnitude of the mechanical loading was imposed on the homeostatic structure.

Simulation 7: Effect of menopause

In this simulation, the mechanical thresholds in Table 1 were arbitrarily increased by 20%, given that the objective of this simulation case was to gain qualitative insight of the effect of menopause. This would result in the rightward movement of the resorption probability curve in Figure 1, which qualitatively means that the resorption probabilities were minutely increased for the mechanical stimuli originally in the adaptive zone, and even more increased for those originally in disuse condition. The mechanical loading condition was the same as that used in simulation 6.

There are numerous physiological bases for this simulation case. It is well known that bone loss at menopause is associated with estrogen deficiency. Researchers found that estrogen directly acts on osteoclasts and regulates the lifetime of osteoclasts [37]. Bone resorption is regulated by both the number of osteoclasts and the activity of each osteoclast. In another bone-remodeling simulation scheme, the bone loss patterns associated with menopause, similar to the clinical observations, have been obtained by increasing the birthrate of new BMUs over the perimenopausal period [38, 39]. Moreover, Nyman et al. increased the mechanostat set point over 3 years to simulate the effect of estrogen withdrawl that produces bone loss observed during menopause [40]. Based on the above investigations, the number of osteoclasts appears to be a more dominant parameter regulating bone resorption, compared with the activity of each osteoblast at menopause. In this study, by increasing the mechanical thresholds in the proposed adaptation model, the increase of resorption probability at menopause can be simulated.

Results

Figures 5(d) to Figure 5(j) show the predicted trabecular architectures for simulations 1 to 7. Figure 5(d) is the numerical result for simulation 1, i.e., disuse mechanical loading condition. Figure 6 gives the percent changes in bone mass from 4 years before, to 4 years after menopause. A 20% end reduction in the external loading (B = 0.8 in Figure 4 is shown to reduce trabecular thickness to a loss of bone mass by 10.48%. The clinical and simulated data are highly correlated (R2 = 0.95). The sum of the squared error (SSE) is 11.68. Hence, our simulation results are consistent with the bone loss patterns clinically observed at menopause [36].

Figure 6
figure 6

Percent changes in bone mass from 4 years before, to 4 years after menopause. The data points show the clinical data measured by Recker et al. [36]. The curves representing percent changes in bone mass are the simulated outcomes. A 20% end reduction in the external loading (B = 0.8 in Figure 4) leads to the lowest sum of the squared error (SSE) from the clinical data.

In simulation 2, increasing the load amplitude by 20% gradually thickened the trabeculae (Figure 5(e)). As a result, the bone mass increased by 16.45%. Accordance is observed with several experiments showing that high-impact gymnastics or stronger muscles increases bone mass [41, 42].

In the case of artificially disconnected trabeculae (simulation 3), two struts in the homeostatic architecture were artificially disconnected (Figure 5(b)). As a result of modeling and remodeling, the unloaded trabeculae disappeared and the adjacent ones thickened with little change in bone mass (Figure 5(f)).

The direction of loads applied to the homeostatic architecture changed in simulation 4 (Figure 5(c)). After rotating the external loads from 30° to a perpendicular orientation, the trabeculae gradually realigned to adapt to the new loading directions (Figure 5(g)).

In the next two simulations, we investigated the bone adaptation behaviors under physiological loading conditions. Visibly, an 8% increase in the external loads only caused 3.12% increase in bone mass with little change in trabecular architecture (Figure 5(h)). Similarly, an 8% decrease in the external loads only resulted in 1.33% decrease in bone mass. The trabecular architecture changed insignificantly (Figure 5(i)).

In the final simulation, the effect of menopause on the mechanical thresholds was investigated. Given that most aging adults lose momentary muscle strength and the intensity of their physical exercise becomes weaker [43], the external loads on the homeostatic architecture were the same with those in simulation 6, i.e., the external loads fell toward the lower end of physiological loading condition. For the mechanical thresholds, a 20% increase resulted in a 7.08% decrease in bone mass. Figure 7 shows the percent changes in bone mass for simulations 6 and 7, from which the effect of mechanical thresholds on bone mass may be clearly seen. The numerical outcome of simulation 7 is shown in Figure 5(j).

Figure 7
figure 7

Percent changes in bone mass for simulations 6 and 7, respectively.

Discussion and Conclusion

The objective of this paper is to understand the contributions of mechanical loading and menopause on the adaptation behaviors of trabecular bone. Here, we developed a numerical bone adaptation model based on bone physiology, in which the amount of bone formation and bone remodeling activation were in accordance with the mechanical levels.

We investigated our bone adaptation algorithm using seven simulation cases, and confirmed that the algorithm is able to predict the adaptation behaviors of trabecular bone under different mechanical conditions. In this study, the formulations on the amount of bone formation and bone remodeling activation probability were proposed based on a number of experimental and numerical investigations. For bone remodeling activation probability, Vahdati and Rouhi [20] suggested that bone remodeling is governed by SED when loading was below a certain threshold value, and by damage when loading exceeded the threshold. In the bone metabolic model of Huiskes et al. [17], the probability of osteoclast activation is considered to be regulated, either by the presence of microcracks within the bone matrix caused by the dynamic forces of daily life, or by disuse since osteoclast activation is enhanced by lack of loading. Tanck et al. [19] investigated the effects of osteoclastic resorption characteristics on trabecular bone remodeling. We proposed the characteristics of bone remodeling activation at different mechanical levels: in disuse zone, the probability of osteoclast activation was controlled by disuse; under physiological loading conditions, it was equal for all surface sites; and in overload zone, it was assumed to be regulated by microdamage. For bone formation, there was more bone formation in highly strained areas, which was similar to the previous numerical algorithms, but the relationship between the amount of bone formation and the mechanical stimulus was not taken as a unique linear form, but a sectioned linear form. The phenomenon that there was no change of bone mass under physiological loading conditions in normal activities was also taken into account. This model could be extended to simulate age-related bone loss. The rate of age-related bone loss is generally between 0.3% and 1.1% per year [44]. In our bone adaptation model, the amount of formed bone can be modeled as a little bit less than the amount of bone resorbed by osteoclasts, resulting in the so-called formation deficit.

In our numerical implementations, the computer simulation was conducted as an iterative process in computer time, in which the relative density per element was regulated between 0.01 and 0.1 (Eq. (5)). The number of iterations was not the actual time (e.g., in the unit of days), but the relationship between the number of iterations and the actual time could be estimated from the beginning and end equilibrium bone masses (see Figure 6 for example) [16]. All simulation results are consistent with the clinical observations and experimental investigations. Hence, by incorporating the proposed bone adaptation model with finite element analysis, trabecular bone adaptations can be simulated at different mechanical levels.

The simulated bone adaptation mechanism is also in accordance with the mechanostat theory of Frost, who proposed that the mechanostat is a biological mechanism that fits skeletal mass and architecture to the needs of normal physical activities of an individual. The mechanical loading made the system correct errors between the mass and its loading. The behaviour of the mechanism is similar to a home thermostat, which would turn "on" in response to error and turn "off" in its absence, hence, the name "mechanostat" [4, 5, 45]. In the present numerical study, the initial trabecular architecture for simulations 1, 2, 5, 6, and 7 was the homeostatic architecture (Figure 5(a)), which was not suitable to the mechanical and/or biological conditions in those simulation cases. This is why the bone adaptation processes have been initiated. Figure 8 shows the average mechanical strains of bone elements in the initial and end trabecular architectures for simulations 1, 2, 5, 6, and 7 under their mechanical loading conditions, respectively. The average mechanical strain in the homeostatic architecture (Figure 5(a)) under its normal mechanical loading was 1041με. When the load decreased by 20% (simulation 1), the average strain in this architecture decreased to 833 με, which was far from that in the homeostatic environment, so the bone adaptation algorithm removed a few bone material, until the average strain rose to 990 με, i.e., the average strain in the new architecture (Figure 5(d)) went back to the adaptive zone. When the external load increased by 20% (simulation 2), the average strain in the homeostatic architecture increased to 1250 με; thus, the bone adaptation algorithm continued to add bone material, until the average strain fell to 1031με to make the new architecture (Figure 5(e)) return to the adaptive zone. Due to the existence of the so-called lazy zone, i.e., a physiological range in which alteration of strain will produce no significant adaptation response, an 8% increase or decrease of external load in simulations 5 and 6, respectively, caused little change of architecture (Figures 5(h) and 5(i)) and strain (Figure 8). Significant attention should thus be given to the comparison between simulations 6 and 7. The external loads in these simulations were the same, but the mechanical thresholds were increased due to menopause in simulation 7. As a result, the average strain in the end configuration (Figure 5(j)) increased by 10.72%, in comparison with simulation 6 (Figure 5(i)).

Figure 8
figure 8

The average mechanical strains of bone elements in the initial and end trabecular architectures for simulations 1, 2, 5, 6, and 7 under their mechanical loading conditions, respectively.

A number of limitations associated with this study may have contributed to our simulations. In the proposed adaptation model, the adaptation behaviors in the pathological overload zone were not included, because in this case, bone tissue showed nonlinear mechanical properties until fracture, which requires further study using nonlinear and fracture mechanics. Another limitation was that the adaptation processes were simulated on a 2 D model, therefore, estimations of adaptation for 3 D samples of bone cannot be obtained. Despite of the limitations mentioned above, the adaptation model for trabecular bone at different mechanical levels proposed in this paper has a number of potential applications. The model may help us gain a better understanding of the relationships between bone morphology and the mechanical, as well as biological environments. Recently, state-of-the-art computational techniques in both hardware and software have been utilized to handle millions of finite elements in the PC base (even hundreds of millions in the supercomputer base) [21–23]. Hence, a further prospect of the bone adaptation model proposed in this paper is to perform a detailed 3 D bone-remodeling simulation. By incorporation with finite element method, the adaptation of bone structure to the mechanical and biological stimuli, e.g., vertebral body, proximal femur, and tibia, which are trabeculae rich skeletal sites that are most at risk of bone loss in elderly population, can be simulated.

Furthermore, this paper provides a methodology for the numerical simulation of changes of bone structural morphology caused by changes of mechanical environment under different circumstances such as orthopaedic surgery, bone internal and external fixations, artificial joints or other implants, as well as changes of biological environment caused by aging, menopause, or pharmacological therapies for osteoporosis. Usually, under such circumstances, bone morphology will be changed to adapt to the new mechanical and biological environments. Nevertheless, approximately more than half a year is required to observe such changes in clinics. Evaluations may then be made as to whether these changes are beneficial for the bone in the long run. To simulate the long-term bone morphological changes by incorporating bone adaptation model with finite element analysis is efficient, thereby decreasing the research period for the related problems and save a large amount of experimental expenditures. Given that the changes of mechanical and/or biological environments will induce some bone tissue to switch to other zones with different adaptation behaviors, our bone adaptation model would be fundamental to such applications.

References

  1. Mullender MG, Huiskes R: Proposal for the regulatory mechanism of Wolff's law. J Orthop Res 1995, 13: 503–512. 10.1002/jor.1100130405

    Article  Google Scholar 

  2. Gong H, Fan YB, Zhang M: Numerical simulation on the adaptation of forms in trabecular bone to mechanical disuse and basic multi-cellular unit activation threshold at menopause. Acta Mech Sin 2008, 24: 207–214. 10.1007/s10409-007-0134-x

    Article  MATH  Google Scholar 

  3. Frost HM: Why do bone strength and "mass" in aging adults become unresponsive to vigorous exercise? Insights of the Utah paradigm. J Bone Miner Metab 1999, 17: 90–97. 10.1007/s007740050070

    Article  Google Scholar 

  4. Frost HM: Proposed general concepts for skeletal adaptation to mechanical usage. Calcified Tissue Int 1987, 45: 145–156.

    Google Scholar 

  5. Frost HM: From Wolff's law to the Utah paradigm: insights about bone physiology and its clinical applications. Anat Rec 2001, 262: 398–419. 10.1002/ar.1049

    Article  Google Scholar 

  6. Rodan GA: Mechanical loading, estrogen deficiency, and the coupling of bone formation to bone resorption. J Bone Miner Res 1991, 6: 527–530. 10.1002/jbmr.5650060602

    Article  Google Scholar 

  7. Westerlind KC, Wronski TJ, Ritman RL, Luo ZP, An KN, Bell NH, Turner RT: Estrogen regulates the rate of bone turnover but bone balance in overiectomized rats is modulated by prevailing mechanical strain. Proc Natl Acad Sci USA 1997, 94: 4199–4204. 10.1073/pnas.94.8.4199

    Article  Google Scholar 

  8. Kelly J: The third culture. Science 1998, 279: 992–993. 10.1126/science.279.5353.992

    Article  Google Scholar 

  9. Cowin SC, Hegedus DH: Bone remodeling I: theory of adaptive elasticity. J Elasticity 1976, 6: 313–326. 10.1007/BF00041724

    Article  Google Scholar 

  10. Hart RT, Davy DT, Heiple KG: A computational method for stress analysis of adaptive elastic materials with a view toward applications in strain-induced bone remodeling. J Biomech Eng 1984, 106: 342–350. 10.1115/1.3138503

    Article  Google Scholar 

  11. Fyhrie DP, Carter DR: A unifying principle relating stress to trabecular bone morphology. J Ortho Res 1986, 4: 304–317. 10.1002/jor.1100040307

    Article  Google Scholar 

  12. Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Slooff TJ: Adaptive bone-remodeling theory applied to prosthetic-design analysis. J Biomech 1987, 20: 1135–1151. 10.1016/0021-9290(87)90030-3

    Article  Google Scholar 

  13. Beaupre GS, Orr TE, Carter DR: An approach for time-dependent bone modeling and remodeling: theoretical development. J Orthop Res 1990, 8: 651–661. 10.1002/jor.1100080506

    Article  Google Scholar 

  14. Weinans H, Huiskes R, Grootenboer HJ: The behavior of adaptive bone-remodeling simulation models. J Biomech 1992, 25: 1425–1441. 10.1016/0021-9290(92)90056-7

    Article  Google Scholar 

  15. Mullender MG, Huiskes R, Weinans H: A physiological approach to the simulation of bone remodeling as a self-organizational control process. JBiomech 1994, 27: 1389–1394.

    Article  Google Scholar 

  16. Zhu XH, Gong H, Zhu D, Gao BZ: A study of the effect of non-linearities in the equation of bone remodeling. J Biomech 2002, 35: 951–960. 10.1016/S0021-9290(02)00028-3

    Article  Google Scholar 

  17. Huiskes R, Ruimerman R, van Lenthe GH, Janssen JD: Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 2000, 404: 704–706. 10.1038/35015116

    Article  Google Scholar 

  18. Ruimerman R, Van Rietbergen B, Hilbers P, Huiskes R: A 3-dimensional computer model to simulate trabecular bone metabolism. Biorheology 2003, 40: 315–320.

    Google Scholar 

  19. Tanck E, Ruimerman R, Huiskes R: Trabecular architecture can remain intact for both disuse and overload enhanced resorption characteristics. J Biomech 2006, 39: 2631–2637. 10.1016/j.jbiomech.2005.08.018

    Article  Google Scholar 

  20. Vahdati A, Rouhi G: A model for mechanical adaptation of tranecular bone incorporating cellular accommodation and effects of microdamage and disuse. Mechanics Research Communications 2009, 36: 284–293. 10.1016/j.mechrescom.2008.10.004

    Article  MATH  Google Scholar 

  21. Jang IG, Kim IY: Computational simulation of trabecular adaptation progress in human proximal femur during growth. J Biomech 2009, 42: 573–580. 10.1016/j.jbiomech.2008.12.009

    Article  Google Scholar 

  22. Jang IG, Kim IY: Computational study of Wolff's law with trabecular architecture in the human proximal femur using topology optimization. J Biomech 2008, 41: 2353–2361. 10.1016/j.jbiomech.2008.05.037

    Article  Google Scholar 

  23. Tsubotaa K, Suzuki Y, Yamada T, Hojo M, Makinouchi A, Adachi T: Computer simulation of trabecular remodeling in human proximal femur using large-scale voxel FE models: Approach to understanding Wolff's law. J Biomech 2009, 42: 1088–1094. 10.1016/j.jbiomech.2009.02.030

    Article  Google Scholar 

  24. Cowin SC, Moss-Salentijn L, Moss ML: Candidates for the mechanosensory system in bone. J Biomech Eng 1991, 113: 191–197. 10.1115/1.2891234

    Article  Google Scholar 

  25. Lanyon LE: Osteocytes, strain detection, bone modeling and remodeling. Calcif Tissue Int 1993, 53(S1):S102-S106. 10.1007/BF01673415

    Article  Google Scholar 

  26. Mullender MG, Huiskes R: Proposal for the regulatory mechanism of Wolff's law. J Orthop Res 1995, 13: 503–512. 10.1002/jor.1100130405

    Article  Google Scholar 

  27. Mullender MG, Huiskes R, Versleyen H, Buma P: Osteocyte density and histomorphometric parameters in cancellous bone of the proximal femur in five mammalian species. J Orthop Res 1996, 14: 972–979. 10.1002/jor.1100140618

    Article  Google Scholar 

  28. Nagaraja S, Couse TL, Guldberg RE: Trabecular bone microdamage and microstructural stress under uniaxial compression. J Biomech 2005, 38: 707–716. 10.1016/j.jbiomech.2004.05.013

    Article  Google Scholar 

  29. Hazelwood SJ, Martin RB, Rashid MM, Rodrigo JJ: A mechanistic model for internal bone remodeling exhibits different dynamic responses in disuse and overload. J Biomech 2001, 34: 299–308. 10.1016/S0021-9290(00)00221-9

    Article  Google Scholar 

  30. Burr DB, Forwood M, Fyhrie DP, Martin RB, Turner CH: Bone microdamage and skeletal fragility in osteoporosis and stress fractures. J Bone Miner Res 1997, 16: 6–15. 10.1359/jbmr.1997.12.1.6

    Article  Google Scholar 

  31. Ruimerman R, Huiskes R, van Lenthe, Janssen JD: A computer-simulation model relating bone-cell metabolism to mechanical adaptaion of trabecular architecture. Comp Meth Biomech Biomed Eng 2001, 4: 433–448. 10.1080/10255840108908019

    Article  Google Scholar 

  32. van Rietbergen B, Weinans H, Huiskes R, Odgaard A: A new method to determine trabecular bone elastic properties and loading using micro-mechanical finite-elements methods. J Biomech 1995, 28: 69–81. 10.1016/0021-9290(95)80008-5

    Article  Google Scholar 

  33. Currey JD: The effect of porosity and mineral content on the Young's modulus of elasticity of compact bone. J Biomech 1988, 21: 131–139. 10.1016/0021-9290(88)90006-1

    Article  Google Scholar 

  34. Hernandez CJ, Beaupré GS, Marcus R, Carter DR: A theoretical analysis of the contributions of remodeling space, mineralization, and bone balance to changes in bone mineral density during alendronate treatment. Bone 2001, 29: 511–516. 10.1016/S8756-3282(01)00613-5

    Article  Google Scholar 

  35. Parfitt AM: The physiologic and clinical significance of bone histomorphometric data. In In Bone Histomorphometry: Techniques and Interpretation. Edited by: Recker RR. Boca Raton, FL: CRC Press; 1983:143–223.

    Google Scholar 

  36. Recker R, Lappe J, Davies K, Heaney R: Characterization of peri-menopausal bone loss: a prospective study. J Bone Miner Res 2000, 15: 1965–1973. 10.1359/jbmr.2000.15.10.1965

    Article  Google Scholar 

  37. Nakamura T, Imai Y, Matsumoto T, Sato S, Takeuchi K, Igarashi K, Harada Y, Azuma Y, Krust A, Yamamoto Y, Nishina H, Takeda S, Takayanagi H, Metzger D, Kanno J, Takaoka K, Martin TJ, Chambon P, Kato S: Estrogen prevents bone loss via estrogen receptor α and induction of fas ligand in osteoclasts. Cell 2007, 130: 811–823. 10.1016/j.cell.2007.07.025

    Article  Google Scholar 

  38. Hernandez CJ, Beaupré GS, Carter DR: A theoretical analysis of the relative influences of peak BMD, age-related bone loss and menopause on the development of osteoporosis. Osteoporos Int 2003, 14: 843–847. 10.1007/s00198-003-1454-8

    Article  Google Scholar 

  39. Hernandez CJ, Beaupré GS, Carter DR: A theoretical analysis of the changes in basic multicellular unit activity at menopause. Bone 2003, 32: 357–363. 10.1016/S8756-3282(03)00037-1

    Article  Google Scholar 

  40. Nyman JS, Yeh OC, Hazelwood SJ, Martin RB: A theoretical analysis of long-term bisphosphonate effects on trabecular bone volume and microdamage. Bone 2004, 35: 296–305. 10.1016/j.bone.2004.03.007

    Article  Google Scholar 

  41. Kalsson MK, Johnell O, Obrant KJ: Bone mineral density in weight lifters. Calcif Tissue Int 1993, 52: 212–215. 10.1007/BF00298721

    Article  Google Scholar 

  42. Courteix D, Lespessailles E, Loiseau Peres S, Obert P, Benhamou CL: Effects of physical training on bone mineral density in prepubertal girls: a comparative study between impact-loading and non-impact-loading sports. Osteopor Int 1998, 8: 152–158. 10.1007/BF02672512

    Article  Google Scholar 

  43. Frost HM: Why do bone strength and "mass" in aging adults become unresponsive to vigorous exercise? Insight into the Utah paradigm. J Bone Miner Metab 1999, 17: 90–97. 10.1007/s007740050070

    Article  Google Scholar 

  44. Van Der Linden JC, Verhaar JA, Weinans H: A three-dimensional simulation of age-related remodeling in trabecular bone. J Bone Miner Res 2001, 16: 688–696. 10.1359/jbmr.2001.16.4.688

    Article  Google Scholar 

  45. Jee WSS: Integrated bone tissue physiology: anatomy and physiology. In Bone Mechanics Handbook. Edited by: Cowin SC. Florida: CRC Press; 2001:1–68.

    Google Scholar 

Download references

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Nos. 10832012, 10872078 and 10972090) and Scientific Advancing Front and Interdiscipline Innovation Project of Jilin University (No. 200903169).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Dong Zhu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

HG carried out the study design, the numerical simulations, and drafted the manuscript. DZ helped to perform the study design and draft the manuscript. JG participated in the numerical simulations. LL participated in the study design and manuscript draft. XZ coordinated the study, participated in the study design and manuscript draft. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Gong, H., Zhu, D., Gao, J. et al. An adaptation model for trabecular bone at different mechanical levels. BioMed Eng OnLine 9, 32 (2010). https://doi.org/10.1186/1475-925X-9-32

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1475-925X-9-32

Keywords