Development of a unifying theory for mechanical adaptation and maintenance of trabecular bone

Trabecular bone is a tissue with a complex 3D structure, consisting of struts and plates, which attains its mature morphology during growth in a process called ‘modeling’. In maturity, the tissue is renewed continuously by local bone resorption and subsequent formation in a process called ‘remodeling’. Both these metabolic activities are executed by bone-resorbing osteoclastic and bone-forming osteoblastic cells. It is known that bone mass and trabecular orientation are adapted to the external forces and that alternative loading conditions lead to adaptations of the internal tissue architecture. The question is how the characteristics of external loads are sensed in the bone, and how they are translated to structural adaptation of the tissue. The time scale of the underlying processes is on the order of months, or even years. This aspect makes bone a complex research topic. In this paper, we discuss the application of computer simulation to investigate the remarkable adaptive processes. We describe our developments of empirical models in the past 15 years, able to predict bone adaptation to external loads from a macroscopic level towards a cell-based level, in which the most important relationships of the cellular processes are captured. The latest model explains the morphological phenomena observed in trabecular bone at a microscopic level.


Introduction
Bone tissue, which forms the skeleton, is a remarkable material. Two different types of bone tissue are distinguished. The first type is cortical or compact bone, which is a low porosity solid material, mainly found in the shaft of long bones. The second type is trabecular bone. It has a complex 3D structure consisting of struts and plates and is mainly found near joint surfaces, at the end of long bones and within vertebrae. It is capable of optimizing its internal tissue structure under the influence of external forces to fulfill its primary function, mechanical load transfer. This paradigm was originally known as Wolff 's Law (Wolff 1892). How the bone actually produces its internal architecture to fulfill its task in an optimized sense is currently unknown, but it is obvious that mechanical feedback must be involved (Huiskes 2000). Trabecular orientation is in the direction of the principal mechanical loads. High loads, such as in intensive exercise, increase bone mass (Courteix et al. 1998), while reduced loads, as in inactivity or disuse, lead to decreased bone mass (Zerwekh et al. 1998). Excessive loss of bone (osteoporosis) leads to increased fracture risk. As the population gets older, metabolic bone diseases such as osteoporosis are becoming increasingly problematic. It often leads to bone fractures, most commonly in the femoral neck, vertebrae or the distal forearm. This results in immobility, social and emotional problems, and high societal costs. Improved understanding of bone cell biology is an important issue in order to prevent osteoporosis and to improve physical and pharmaceutical treatment methods, as well as prosthetic designs. Of course it also involves many complex biochemical processes of which much is unknown.
One approach to obtain insight in bone cell biology is by isolating the different components to unravel their individual (and often very complex) behaviour. Although this cannot fully explain the behaviour of bone as a whole organ, it can provide valuable information. It is, however, difficult to take all aspects of bone biology into account in realistic conditions. In vivo experiments are hard to control and they still have limited life spans; insufficient to study anything but the earliest stages of adaptive responses, as the time scales at which the relevant processes occur are on the order of months, or even years. An alternative approach to gather knowledge of the processes involved in bone remodeling is by computer simulation. In this paper, an overview is presented from work directed by the second author concerning computer simulations in bone biomechanics, starting with an empirical model that was able to predict density changes in bone under the influence of external forces, and its subsequent expansion towards a more mechanobiologically oriented model that is cell-based and explains bone behaviour and morphology on the scale of individual trabeculae as effects of force transmission. For a wider historical survey, refer to Hart and Fritton (1997).

Bone density controlled by mechanical factors
With increasing computer capacity, computer simulation of bone adaptation in whole bones became practical in the late 1980s. The early models, based on finite element analysis (FEA), calculated mechanical signals within the bone that were assumed to initiate changes in density or material properties, for which process numerical formulations of adaptive bone remodeling theories were applied. The models were empirical and related mechanical signals, like stress or strain, to bone adaptation, without direct consideration of the underlying cell-biological mechanisms. They were applied to complete bones and simulated adaptation of bone tissue on a macroscopic level.
In the work of Fyhrie and Carter (1986a), the bone tissue was considered as a continuum, characterized by particular apparent density value distributions. The bone was assumed to be a self-optimizing material with the objective to adapt its apparent density to an 'effective stress' eff . Based on strength or strain optimization criteria they derived where A and are constants, and the effective stress is determined from either a failure or an elastic energy criterion. In later publications (Carter et al. 1987, Fyhrie andCarter 1986b), they assumed ¼ 0.5 for the effective stress, and postulated where E is the apparent elastic modulus and U the apparent strain energy density (SED). By assuming a modulus-density relation of E ¼ c 3 (Carter and Hayes 1977), with c a constant, the optimization function transforms to where c 0 is a constant. Applying this optimization criterion made it possible to predict a density distribution of the proximal femur assumed optimal (Carter et al. 1987, Fyhrie andCarter 1986b). To test the theory, they visually compared the resulting density patterns to those in radiograms of a real femur, and observed some similarity in it. Frost (1964) suggested that internal and external remodelingy should be distinguished. Internal remodeling is the adaptation of the density of bone tissue, while external or surface remodeling is the apposition or removal of bone tissue on the bone surface. The two forms were actually separated by Cowin and associates and by Huiskes and co-workers. Cowin et al. (1981Cowin et al. ( , 1985 used the strain tensor as the mechanical signal to be the driving factor of bone adaptation, and assumed a quadratic relationship between strain and the rate of adaptation . Huiskes et al. (1987) used the strain energy density (SED) U [J/mm 3 ] as the signal that controls remodeling of the bone, with

Internal and external remodeling
The relation between adaptive rate and SED was linear. For external surface remodeling, the bone can either add or remove material, according to where dX/dt is the rate of bone growth perpendicular to the surface, U the SED, U n a homeostatic SED, and C x a proportionality constant. For internal remodeling the bone could adapt its density value. Assuming that the elastic modulus relates to the apparent density one can write where E [MPa] is the local elastic modulus, and C e a proportionality constant. Carter (1984) suggested that bone is 'lazy' in terms of reacting to mechanical signals. The concept of a so-called 'lazy zone', meaning that there are thresholds to yAlthough 'modeling' and 'remodeling' of bone were clearly defined by Frost (1964) as change in shape (as in growth) and the subsequent localized resorption and formation of bone, the term 'remodeling' is often used in biomechanics as if to mean 'adaptation'. In fact, 'modeling' is closer to adaptation than 'remodeling'. be exceeded before bone adaptation could occur (figure 1), were incorporated by Huiskes et al. (1987). The SED was compared to a homeostatic one, U h . For U > (1 þ s) U h or U<(1 À s) U h adaptive activity is initiated, whereby s is the threshold level that marks the borders of the lazy zone. The remodeling rule transforms to a new set of equations for internal remodeling and a similar set of equations for external remodeling. This theory was applied to predict the density distribution in the normal proximal femur. Starting from an initial uniform density distribution the simulation produced a configuration that was similar to the one obtained by Fyhrie and Carter (1986b), showing some similarity to the natural density distribution in the femur. The theory was also applied to predict cortical adaptation after hip-prosthetic implantation in a geometrically idealized model (Huiskes et al. 1987). The bone was represented as a hollow cylinder in which a cylindrical stem was placed and the effects of 'stress shielding' due to load transfer through the stem could be predicted. The results illustrated the potential applicability of computer simulations for periprosthetic bone adaptation and encouraged further investigation. Later the theory was successfully applied to evaluate bone adaptation in a realistic 3D femur model after implantation of prostheses. Three-dimensional finite element models were constructed from an experimental animal configuration, in which smooth press fitted stems were applied. The adaptive remodeling procedure was integrated with the model and predicted similar amounts of proximal bone loss and distal bone densification as found in animal experiments . The prosthetic stress-shielding theory was also validated relative to postmortem patient series (Kerner et al. 1999).

Discontinuous end configurations
In the remodeling rules discussed in the previous section, bone adapts towards reference states of deformation caused by single loads. To incorporate variable loading, a stimulus S was defined that took account of individual loading configurations by averaging their individual contributions (Carter et al. 1989), as where U i is the SED value for loading case i, n the total number of loading cases and the apparent density. In order to estimate the actual SED values in the trabeculae, the remodeling signal was approximated by U/ as the strain energy per unit of bone mass (Carter et al. 1989). In 1992, Weinans et al. applied this theory to a two dimensional finite element model of a proximal femur. Bone was represented as a continuum, capable of adapting its apparent density due to mechanical stimulation. The stimulus value S was measured per element, and the apparent density was adapted per element according to: with B and k constants and S the stimulus value in the element. This theory produced density distributions that showed good resemblance with those in a real proximal femur. However, only if these distributions were locally averaged from discontinuous density patterns in the femoral head, where the trabecular bone is located (figure 2). In fact, the underlying patchwork of either full or empty elements is inadmissible relative to the mathematical basis of FEA. However, Weinans et al. (1992) also showed that it is in the nature of the differential equations used to mathematically describe the adaptive remodeling process that the simulation produces discontinuous end configurations, a phenomenon called 'checker boarding'. Using an analytical two-unit model they showed that the differential equations have an unstable uniform solution and that the slightest non-uniformity causes stress shielding by the stiffer elements. Hence, in the ongoing process the stiffer elements will even become stiffer, while the others loose density, until they are virtually empty, so that a 'checkerboard' results (figure 2). This phenomenon is the result of capturing both sensation of the local mechanical stimulus and actual adaptation of bone mass in that same location in one equation, which make the elements--which are artifacts for biology--'work for themselves'.

Osteocytic mechanosensation
Besides being mathematically unstable, the previous theory was also biologically naı¨ve, because osteoblasts and osteoclasts-the bone producing and removal cells-do not reside within the bone, but only appear at the trabecular surface as an effect of mechanical signals, which must be sensed within the bone matrix. Hence, the sensation of mechanical stimuli and adaptation of material properties should be separated in the theory, as was shown by Mullender et al. (1994) and Mullender and Huiskes (1995). Osteocytes within the tissue matrix were assumed to be mechanosensitive and capable of translating signals to the bone surface to attract so-called basic multicellular units (BMUs) (i.e. osteoclasts and osteoblasts, the actor cells of bone biology) which control the net apposition or removal of bone tissue. The signal sent to the surface by an osteocyte was assumed to decay exponentially with increasing distance d [mm], according to where x is the location in which the signal strength is determined, x 0 is the location of the osteocyte concerned and D [mm] determines the decay in signal strength.
The relative density m [À] at location x was regulated by the BMUs. Relative to the total amount of stimulus P in that location, they adapt the bone density according to where [MPa À1 s À1 ] is a rate constant. All osteocytes N located within the region surrounding x contribute to the stimulus P, so that its value in x is determined by Here, x i is the location of osteocyte i, S(x i ,t) is the mechanical signal this osteocyte senses, and k is a reference value. Although hypothetical, the mechanosensory role of osteocytes is not controversial. Osteocytes are shown to be sensitive to mechanical loading. They derive from osteoblasts that are encapsulated in the tissue matrix they produce. They are interconnected by a network of canaliculi, that seem perfectly equipped for signal transmission (see also Burger and Klein-Nulend 1999). A two-dimensional computational model was developed for this theory and applied at a microscopic level to a square plate of 2 Â 2 mm 2 , divided in 80 Â 80 elements. The separation of sensation and density adaptation into two different functions prevented the so-called checker-boarding phenomenon. The theory produced configurations with characteristics that resemble actual trabecular bone structures, with average density correlating to the loading magnitude and trabecular orientation directly related to the external loading direction (figure 3). The influence parameter D in equation (10) had a prominent role and variations affected the resulting morphology. The smaller the influence distance D, the more refined the resulting architectures, with thinner trabeculae. The model was able to realign its trabeculae to alternative loads and so explained both modeling and adaptation. After the structure was optimized or adapted to the external loads, the activity of the BMUs stopped.

Separation of osteoclastic and osteoblastic activity
In actual trabecular bone however, the activity of both osteoclasts and osteoblasts continues throughout life. In the growth stage, or during adaptation, these cells actually shape the bone structure (modeling), while in homeostatic conditions (or adulthood) they maintain the bone structure by tissue renewal through constant osteoclastic resorption and subsequent refilling of the cavities by osteoblasts (remodeling). In metabolic diseases (e.g. osteoporosis), it is the activities of these cells that are deregulated. How are the separate actions of these cells controlled? How is coupling between activity of these cells established, and how does this relate to bone maintenance? In order to answer such questions, a more refined theory was required, in which the separate activities of osteoclasts and osteoblasts were incorporated. Such a theory was developed by Huiskes and associates , Ruimerman et al. 2001. The new theory is based on the regulation scheme depicted in figure 4 and the relationships involved were quantified. Tissue adaptation on a specific location x at time t was supposed to be the result of osteoclastic bone resorption and osteoblastic bone formation, hence What exactly initiates resorption is not known. In the theory, it was assumed that osteoclasts are attracted towards the bone surface by microdamage, and that microdamage occurs spatially random, so that bone resorption is determined according to where r cl [mm 3 /day] represents a stochastic function. Osteoblastic activity was controlled by an osteocytic bone formation signal P [mol mm À2 day À1 ]. If the stimulus Figure 4. Scheme of the proposed regulation mechanism: (i) osteoclasts are assumed to resorp bone on locations where microdamage occurs which is thought to be on stochastic locations on the trabecular surface. (ii) osteocytes locally sense a mechanical signal due to external load transfer through the architecture and locally recruit (iii) osteoblasts to form bone tissue . exceeds a certain threshold value k tr [mol mm À2 day À1 ], there is assumed to be tissue formation at the trabecular surface according to where [mm 5 mol À1 ] is a proportionality factor that determines the formation rate. The formation stimulus P in location x is determined by all osteocytes N located within the influence region, the signal R [J mm À3 s À1 ] sensed by each osteocyte i, with mechanosensitivity i [mol mm J À1 s day À1 ] and location x I , according to where f(x, x i ) determines the decay in signal strength according to equation (10). Until now, computational models of bone adaptation used static loads to evaluate the mechanical signals. It is, however, known that bone only reacts to dynamic loads. In the new theory, the stimulus sensed by osteocytes is assumed to be a 'typical' strain energy density rate (SED-rate) R(x, t) in a recent loading history. Cyclic loading conditions, characterized by frequency and magnitude, were imposed and it was assumed that osteocytes react to the maximum SEDrate during the loading cycle. It was shown that the maximal SED-rate is related to the SED value for some substitute static load and that it could be calculated by static finite element analysis (Ruimerman et al. 2001). The relationships were incorporated in a two-dimensional computer simulation model. Starting from different initial configurations the structures remodeled toward similar homeostatic configurations, in which the process of resorption and formation continued, but no more architectural changes occurred (figure 5). Trabeculae were aligned to the loading direction. The regulatory mechanism was able to adapt the structure to alternative loading conditions. Increasing the magnitude led to apposition of bone tissue while decreasing it led to bone loss. After rotating the direction of the external load, the trabeculae realigned their orientations accordingly ( figure 5). From the simulation model, it was concluded that the theory was able to explain both modeling (formation and adaptation) as well as remodeling (maintenance) of trabecular-like architectures as governed by external forces. Several additional conclusions could be drawn from the theory. Osteoblasts and osteoclasts are often assumed to work together in so-called BMUs. In our theory, coupling between these cells is established only implicitly through mechanics. The resorption cavities produced by osteoclasts cause local stress and strain concentrations. Consequently, osteocytes within the tissue matrix sense elevated mechanical signals and locally recruit osteoblasts to form bone tissue until the mechanical signal has decreased to normal levels, i.e. until the cavity is refilled ( figure 6). This does not imply that biochemical pathways are irrelevant (Udagawa et al. 1999), but it does show that mechanical feedback can be a potent coupling factor for the relevant biochemical processes to take place. Another conclusion came from the observation that both modeling (growth and adaptation) and remodeling (maintenance), can be explained by one unifying theory. Indeed, both processes are based on the same cellular mechanisms and it could be concluded that both modeling and remodeling are similar in nature and just different stages of the metabolic cascade.

Recruitment stimulus
Obl Elevated strain Ocl Ocy Figure 6. In the theory, osteoblast-osteoclast coupling is established as illustrated. Load transfer around osteoclastic resorption cavities results in elevated mechanical signals as sensed by osteocytes. Consequently, osteoblasts are recruited to form bone tissue until the mechanical signals have decreased to normal levels, i.e. when the resorption cavity is filled. During bone formation, some osteoblasts are entrapped in the tissue matrix they produce and differentiate to osteocytes. The remaining osteoblasts become inactive lining cells that completely cover the trabecular surface.

Towards three-dimensional structures
Applied in a 2D computer simulation model, the theory could explain morphological adaptations observed in trabecular bone as related to load-induced cell activities.
Only, however, in a conceptual sense. To actually validate the theory in a quantitative sense, and to apply it for valid predictions that are useful and falsifiable, its application in a 3D geometry is required. Computer capacity is a limiting factor for this requirement and it is currently not possible to simulate bone modeling and remodeling in complete bones at the trabecular level. In order to test whether our theory can produce trabecular-like structures in three dimensions, with morphological characteristics in a quantitative realistic range, for actual trabecular bone, we developed a 3D computer simulation model that can be applied to a small domain of bone tissue. Preliminary results show that the theory mimics the realistic morphological expressions of bone cell metabolism in a robust way, in growth, adaptation and remodeling, when applied to a small cube (2 Â 2 Â 2 mm 3 ) of bone tissue with external loads imposed. When started from an initial porous configuration, the model produced structures with trabecular alignment in the loading direction (figure 7), while volume fraction correlated to the magnitude of the external loads.
The model described the growth as well as adaptation of complex 3D structures that are similar to actual trabecular bone. The architectures produced are remarkably stable, while the remodeling process proceeds.

Discussion
Initially, theories for bone adaptation related bone density changes directly to mechanical strain-derived variables as effects of external forces. The computational simulation models based on these theories proved useful tools for prosthetic design. They were empirical models on a macroscopic level and the underlying biological processes involved were not considered. In the past decade, the focus has shifted from biomechanics to mechanobiology, as it is likely that actual intervention in metabolic processes-to cure bone diseases-will be based on understanding the relationship between mechanics and cell biology. Understanding of bone biology has improved enormously in the past 15 years. This, combined with increased computer capacities and efficient finite element algorithms, has enabled the application of realistic theories, including the main aspects of bone biology. The theory relates cellular activity in bone modeling and remodeling to external forces and explains adaptive behaviour of trabecular bone at a microscopic level. Although the most important relationships are captured in the theory, it should still be considered as a framework. The relationships are captured in simple mathematical equations, but actually encompass complex biological processes, the biochemical components of which are largely unknown. Nevertheless, the computational models based on the theory will enable us to investigate the morphological consequences of alternative loading conditions, metabolic disorders and their pharmaceutical interventions. They are useful tools for the development and investigation of hypotheses and for efficient design of experiments. Potentially, they are useful for prosthetic design, as the adaptation of bone tissue to alternative mechanical environments can be pre-clinically tested. They can also be used in the development of treatment methods of osteoporosis, either based on physical exercise or in combination with pharmaceutics. In the near future, the 3D computational model will be used in order to investigate questions like: can the morphological changes seen in senile osteoporosis be explained by reduced physical activity alone? Can anti-resorptive drug administration prevent loss of bone tissue and preserve trabecular connectivity, and at what stage should such a treatment start? Another area where the model can be applied is in postmenopausal osteoporosis. Can the theory mimic the phenomena observed in postmenopausal osteoporotic bone (loss of bone density and trabecular connectivity at a higher remodeling rate) and what is the role of mechanical coupling between osteoclasts and osteoblasts for this condition? Can osteoporosis be an effect of reduced osteocyte mechanosensitivity? Can the morphology of osteopetrotic bone be predicted as an effect of dies-functioning osteoclasts? Can the loss of trabeculae in older ages be prevented? Could the effects of physical activity be replaced by high frequency vibrations? Can trabecular and cortical (osteonal) bone modeling and remodeling processes be described with the same theory, and can the theory explain the transformation of trabecular into cortical bone from mechanobiological stimuli, as it occurs in growth? These are all important questions in bone biology and pharmacology.