Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic

ARTICLE IN PRESS

Journal of Biomechanics 38 (2005) 2440–2450 www.elsevier.com/locate/jbiomech www.JBiomech.com

Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic Sandra J. Shefelbine, Peter Augat, Lutz Claes, Ulrich Simon Institute for Orthopaedics and Biomechanics, Helmholtzstrasse 14, 89081 Ulm, Germany Accepted 7 October 2004

Abstract Trabecular bone fractures heal through intramembraneous ossification. This process differs from diaphyseal fracture healing in that the trabecular marrow provides a rich vascular supply to the healing bone, there is very little callus formation, woven bone forms directly without a cartilage intermediary, and the woven bone is remodelled to form trabecular bone. Previous studies have used numerical methods to simulate diaphyseal fracture healing or bone remodelling, however not trabecular fracture healing, which involves both tissue differentiation and trabecular formation. The objective of this study was to determine if intramembraneous bone formation and remodelling during trabecular bone fracture healing could be simulated using the same mechanobiological principles as those proposed for diaphyseal fracture healing. Using finite element analysis and the fuzzy logic for diaphyseal healing, the model simulated formation of woven bone in the fracture gap and subsequent remodelling of the bone to form trabecular bone. We also demonstrated that the trabecular structure is dependent on the applied loading conditions. A single model that can simulate bone healing and remodelling may prove to be a useful tool in predicting musculoskeletal tissue differentiation in different vascular and mechanical environments. r 2004 Elsevier Ltd. All rights reserved. Keywords: Metaphyseal fracture healing; Trabecular remodelling; Finite element analysis; Fuzzy logic; Tissue differentiation

1. Introduction Trabecular bone is found at the ends of long bones, in the vertebral bodies, and in flat bones such as the ribs and scapula. The microscopic structure of trabecular bone differs from cortical bone found in the diaphyseal shafts of long bones. In particular, trabeculae form a porous framework of bone spicules, often classified (and often quantified) as ‘‘rods’’ or ‘‘plates’’ (Hildebrand and Ruegsegger, 1997). The porous space between the trabeculae is filled with bone marrow, providing the bone with a rich vascular environment. In trabecular bone fractures the trabeculae collapse, break, or are crushed, reducing the strength of the structure (Hayes Corresponding author. Tel.: +49 731 500 33754; fax: +49 731 500 23498. E-mail address: [email protected] (S.J. Shefelbine).

0021-9290/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2004.10.019

and Carter, 1976). Trabecular bone fractures are common in the vertebrae and in the metaphyseal and epiphyseal regions of long bones, such as intratrochanteric or sub-trochanteric fractures of the proximal femur. Elderly people with osteoporosis are particularly at risk for trabecular bone fractures as the integrity of the trabecular structure is compromised in osteoporotic bone. Trabecular bone fractures heal through a different process than most diaphyseal fractures. In diaphyseal fractures that heal through endochondral ossification a callus forms around the fracture, soft tissue becomes cartilage and is then ossified. In trabecular bone fractures, bone forms directly through intramembraneous ossification with almost no cartilage intermediary (Jarry and Uhthoff, 1971; Uhthoff and Rahn, 1981). In intramembraneous ossification bone is laid down directly on the existing bone. In the first stage of

ARTICLE IN PRESS S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

2441

Fig. 1. Histology of canine trabecular bone osteotomy gap healing. (A) After 1 week there is formation of new bone around the old trabeculae. The fracture gap is filled with soft tissue. (B) After 3 weeks woven bone has started to form in the osteotomy gap. (Only half of the gap is pictured.) (C) After 4 weeks the osteotomy gap is completely filled with woven bone. This bone then remodels to form trabecular bone. (Histologies courtesy of Hans K. Uhthoff.)

trabecular fracture healing, there is resorption of necrotic tissue as soft tissue forms in the fracture (or osteotomy) gap. Osteoblastic activity begins around the existing trabeculae on either side of the fracture gap (Fig. 1A). The thickening of the trabeculae is often seen on radiographs as a bright band directly adjacent to a radiolucent band, the fracture gap. After the existing trabeculae are reinforced, osteoblasts begin to lay down bone in the fracture gap (Fig. 1B). This initial bone that forms is called woven bone, has no apparent trabecular orientation, and is relatively homogeneous. Radiographically, two bright bands of bone on either side of the gap approach each other, narrowing the gap. The bridging of the gap with woven bone (Fig. 1C) occurs before the bridging of the cortical shell. After the gap is stabilized and filled with woven bone, osteoblasts and osteoclasts start the remodelling process. The relatively isotropic, homogeneous woven bone is remodelled into an anisotropic, inhomogeneous trabecular bone with the trabeculae orientation dependent on the loading conditions. In radiographs, the bright band of bone in the gap slowly disappears as the bone is remodelled. Numerous previous studies have used numerical methods to simulate diaphyseal fracture healing. The models implement mechanobiological theories, proposing that the mechanical environment influences tissue differentiation. It has been proposed that the stress and strain (Pauwels, 1960; Carter et al., 1988; Claes and Heigele, 1999), strain and fluid velocity (Prendergast et al., 1997), or concentration of growth factors (BailonPlaza and van der Meulen, 2001) regulate tissue differentiation. Fuzzy logic has also has been used in combination with finite element models to simulate progressive tissue differentiation (Ament and Hofer, 2000; Simon et al., 2003). In addition to modelling fracture healing, numerical models have been used to simulate adaptation of trabecular bone density (Carter et al., 1987; Huiskes et al., 1987; Fyhrie and Carter, 1990; Weinans et al., 1992; Cowin et al., 1993) and

structure (Mullender and Huiskes, 1995; Turner et al., 1997; Adachi et al., 2001; Koontz et al., 2001; Ruimerman et al., 2003). These models propose that the remodelling process is dependent on the mechanical environment. In particular, the strain energy density or the strain energy gradient was used as a stimulus for the remodelling process. Previous models, however, have not simulated the fracture healing process in trabecular bone. It has also not been demonstrated that a single model can simulate both fracture healing and trabecular bone formation. The objective of this study was to determine if trabecular bone fracture healing could be simulated using the same mechanobiological principles of those proposed for diaphyseal fracture healing. In particular, we wanted to simulate the de novo formation of woven bone in the fracture gap and remodelling of woven bone to load oriented trabecular bone. The numerical methods were based on a previous model that predicted endochondral ossification in diaphyseal fracture healing (Simon et al., 2003).

2. Methods 2.1. Model We modelled a fracture gap in trabecular bone using a micro-scale idealized cubic section (600 mm on a side) of the fracture gap consisting of approximately 38,000 elements (Fig. 2). We modelled it on the microscopic level in order to capture the microstructure of the trabeculae. The fracture gap was 500 mm in width and bordered by trabeculae spicules on both sides. The model was similar to previous mFEA voxel models of trabecular bone in which the individual trabeculae were modelled (Adachi et al., 2001; Koontz et al., 2001; Ruimerman et al., 2003). The trabeculae bordering the fracture gap were idealized as rectangular blocks four

ARTICLE IN PRESS 2442

S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

Fig. 2. Finite element model of a 0.5 mm osteotomy gap consisting of trabecular bone and soft tissue. The cross-section displays the trabeculae spicules that border the osteotomy gap.

elements wide and three elements long, giving a trabecular thickness of 100 mm. The space between trabeculae and in the fracture gap was filled with elements of soft tissue to represent marrow and the formation of the haematoma in the gap. Finite element modelling and calculations were performed using ANSYS (version 7.11, Canonsburg, PA). 2.2. Material properties All material properties were linear elastic, isotropic, and homogeneous. The initial model consisted of soft tissue and bone. Trabecular bone modelled on the tissue level has an apparent modulus of 500–5000 MPa (Carter et al., 1987; Van Rietbergen et al., 1995). However, we modelled the trabecular bone on the microscopic level, and material properties of a single trabecula are similar to cortical bone (Muller and Ruegsegger, 1995; Van Rietbergen et al., 1995). As the simulation progressed, the material properties for each element were updated. Each element was considered to have a mixture of soft tissue (E ¼ 3 MPa; n ¼ 0:3), bone (E ¼ 10; 000 MPa; n ¼ 0:3), and cartilage (E ¼ 10 MPa; n ¼ 0:49) (Hori and Lewis, 1982; Rho et al., 1995). Mixture rules were applied in order to determine Young’s modulus (E) and Poisson’s ðnÞ ratio for each element (Ament and Hofer, 2000; Lacroix et al., 2002). The elastic modulus was calculated as a cubic function of the tissue volume fractions (Carter and Hayes, 1977): E el ¼ E cart c3cart þ E tiss c3tiss þ E bone c3bone;

(1)

where ccart, ctiss, and cbone were the volume fractions of cartilage, soft tissue, and bone respectively in the element. Poisson’s ratio was calculated as a simple weighted average nel ¼ ncart ccart þ ntiss ctiss þ nbone cbone :

(2)

The volume fractions of each tissue were determined from the fuzzy logic output (described below) and were such that ccart þ ctiss þ cbone ¼ 1:

(3)

The material properties were updated for each iteration of the simulation. In addition to mechanical material properties, each element was assigned a vascular perfusion from 0 to 1, with 1 representing maximum vascularity. All elements were initially assigned maximum vascularity to represent well-vascularized trabecular bone. During the simulation the vascularity was a function of the mechanical environment (Table 1). 2.3. Loading The intricate structure of normal trabecular bone forms through complicated (and varying) loading conditions. As it is not possible to determine the in vivo load on a single trabecula over time, we assumed equal pressure to the top of each trabecula on one side of the fracture gap while holding the trabeculae on the other side of the fracture gap fixed. The soft tissue elements were not constrained. The applied pressure was assumed to vary during the consolidation of the fracture (Fig. 3). In this model we assumed rigid fixation during the initial healing phase and therefore stress shielding occurred at the fracture site. During the first stage of healing before the gap was bridged (iterations 1–15), the load was maintained at 0.05 MPa. The low pressure represented the initial healing phase, when most of the load in the bone was supported by the osteosynthesis (a rigid plate or fixator) that stabilized the fracture. Once the gap was filled with woven bone, the pressure on the trabeculae increased (iterations 16–20) and was maintained at 3 MPa (iterations 20–140). This represented increased the weight bearing on the fractured limb and

Table 1 Fuzzy rules implemented in the fuzzy controller IF

% bone

% cartilage

1.

3.

Max vascularity neighbour

Max % bone neighbour

Destructive

5. 6. 7. 8. 9. 10.

Low

11.

Low Not high

Low

13.

Not low

14.

High

15. 16. 17. 18.

Not Not Not Not

19.

High

Low

20.

High

Low

low low low low

Negative Medium medium Positive Medium medium Negative high Not destructive Negative high Not destructive Negative Not medium destructive Medium Low Negative high Low Negative high Medium Negative Medium medium Negative Medium medium Negative Low medium Low Low

Vasculartiy

Bone

Cartilage

Decrease

Decrease

Decrease

Decrease

Decrease

Decrease

Decrease

Decrease

Decrease

No change

No change

No change

No change

Low

Not low

Increase

Medium

High

Increase

High

9 No change > = No changes > No change ; 9 > > > > > > > > > = Vascularity changes > > > > > > > > > ;

Increase Decrease

High

Not low

Increase

High

Not low

Increase

9 > = Intra-membraneous ossification > ; Increase Increase Increase

Not Not Not Not

low low low low

Not Not Not Not

low low low low

Increase Increase Increase Increase

Decrease Decrease Decrease Decrease

Not low

High

Increase

Decrease

Not low

High

Increase

Decrease

High

9 > > > > = Tissue distruction > > > > ;

Decrease

9 > > > > = Cartilage formation > > > > ; 9 > > > > > > > > > > > = Endochondrial ossification > > > > > > > > > > > ; o Remodelling

The rules are implemented as if/then statements. Each input and output has linguistic ranges. The inputs % bone, % cartilage, % vascularity have the possible ranges low, medium, high. Hydrostatic strain has the ranges negative destructive, negative high, negative medium, low, positive medium, positive destructive. Octahedral shear strain has the ranges low, medium, high, destructive. The biological event that occurs when a rule is activated is indicated on the right-hand side.

ARTICLE IN PRESS

Not negative destructive Not positive destructive Not negative Not high high Not negative Not high high Not negative Not high high Negative high

THEN

S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

4.

21.

Octahedral Vascularity shear strain

Positive destructive Negative destructive

2.

12.

Hydrostatic strain

2443

ARTICLE IN PRESS S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

2444 woven bone formation

trabeculae bone formation

applied pressure (MPa)

3

0.05 0

20

40

60

80

100

120

140

iteration number Stage: 1

2

3

Fig. 3. Applied loading conditions over time. During the first 15 iterations (Stage 1) woven bone filled in the gap. The bone matured during iterations 16–20 (Stage 2) as the load increased. Resorption of bone in these first two stages was not allowed. After 20 iterations the fracture was stable and the load was maintained at 3 MPa (Stage 3). During this stage resorption was allowed and remodelling of the woven bone into trabecular bone occurred.

removal of the stabilization device once the gap was stable (Augat et al., 1996). Similar load magnitudes have been applied in trabecular remodelling simulations (Huiskes et al., 2000). Loads were applied to the trabeculae and the hydrostatic (dilatational) and octahedral (deviatoric) strain were determined for each element in the model. A second load case was used to determine the effects of load direction on the remodelling of trabecular structure. This load case was equivalent to the first one until iteration 20 (stage 3) at which time a diagonal load was applied to three of the trabeculae spicules. The diagonally opposite trabeculae spicules were fixed in all directions. There was no load on the remaining trabeculae spicules or other elements. The magnitude of the total load was 0.4 N applied at a 451 angle to the z-axis, resulting in a total pressure of 1 MPa on the entire cube. 2.4. Fuzzy logic In order to determine tissue differentiation over time we used a fuzzy logic controller, consisting of a set of 21 linguistic rules to describe tissue destruction, angiogenesis, endochondral ossification (including chondrogenesis and calcification of cartilage to form bone), intramembraneous ossification, and remodelling (bone resorption) (Simon et al., 2003). In general, the fuzzy controller received input from the finite element analysis, determined the location of the strain state on the tissue differentiation diagram, and decided what

type of tissue to build. In addition to the strain state, surrounding tissue properties and vascularity are also critical to tissue differentiation, thus requiring the use of multiple rules and fuzzy logic. The fuzzy controller required seven input variables for each element in the model: (1) percentage of bone in the element, (2) percentage of cartilage in the element, (3) hydrostatic strain, (4) octahedral shear strain, (5) vascularity of the element, (6) the maximum percentage bone in a neighbouring element, and (7) the maximum vascularity of a neighbouring element. Each rule was formulated as an if/then statement: IF these conditions are present, THEN this tissue differentiation occurs (Table 1). The output of the fuzzy controller was the percentage change of bone, cartilage, and vascularity for each element. The Matlab Fuzzy Logic Toolbox (Mathworks, Inc., Natick, MA) was used for all fuzzy logic calculations. The fuzzy rules were adapted from previous studies that used fuzzy logic to determine tissue differentiation during diaphyseal healing (Ament and Hofer, 2000; Simon et al., 2003) and based on the tissue differentiation diagram of Claes and Heigele (1999) (Fig. 4). There were five main areas of the diagram: (1) When the hydrostatic or octahedral strain is too high (usually out of the physiological range), tissue destruction occurred. (2) With relatively high-hydrostatic compressive strain, cartilage formed. (3) With high-hydrostatic tensile or octahedral strain, soft tissue formed. (4) With small amounts of hydrostatic and octahedral strain, bone formed. (5) However, if the strain stimulus was too small, bone resorption occurred. The fuzzy rules were a linguistic description of the tissue differentiation diagram (e.g. regions of high, medium, and low strain). For example one rule (#10, Table 1) read, ‘‘IF the percentage cartilage is low, the hydrostatic strain is negative medium, the octahedral strain is medium, the vascularity is high, and the percentage bone in a neighbouring element is not low, THEN the percentage of bone in the element increases.’’ All conditions must be met in order for the rule to be activated. The regions were assigned numerical ranges in order to translate the finite element strain state into the appropriate region of the diagram (Fig. 5). For strain values between two regions (shaded areas in Fig. 5), multiple rules could be activated simultaneously, creating the ‘‘fuzzy’’ part of the logic. In effect the boundaries between the different regions in the tissue differentiation diagram were not strict thresholds, but allowed for the differentiation to go both ways. The amount a particular rule was activated was determined by where in the overlapping region the value fell. The final output was the sum of the outputs of all rules that were activated. The numeric values for the regions were chosen based on previous in vitro experiments, in vivo measurements, and modelling simulations (Rubin and Lanyon, 1987; Frost, 1992; Claes and Heigele, 1999; Kaspar et al.,

ARTICLE IN PRESS S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

2445

ε s octahedral shear strain

tissue destruction

destructive

connective tissue

high

cartilage medium

bone

bone resorption 0

ε h hydrostatic strain

positive high

positive destructive

tension

positive medium

low

negative high

negative medium

compression

negative destructive

low

Fig. 4. Tissue differentiation diagram based on Claes and Heigele (1999). The boundaries between the regions are ranges rather than thresholds, making it ‘‘fuzzy.’’ negative high

negative very high

negative near middle zero

positive middle

positive very high

probability

1

0.05

0.02

-0.001 -0.0001 0.0001 0.001

-0.008

-0.01

-0.02

-0.05

0

Hydrostatic strain low middle

high

very high

probability

1

0.17

0.15

0.06

0.05

0.0001 0.001

0

Octahedral shear strain

Fig. 5. Linguistic translation of hydrostatic and octahedral shear strain. Shaded areas correspond to the boundaries between the different tissue regions of Fig. 4. These shaded areas make the logic ‘‘fuzzy’’ as a value that lands in the shaded area (near a boundary) may activate multiple rules.

ARTICLE IN PRESS S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

2000). The output from the controller included percent changes in bone, cartilage, and vascularity for each element. In addition to the mechanical stimuli there are biological factors that play important roles in fracture healing. During healing the fracture gap is flooded with local and systemic growth factors (e.g. BMPs, TGF-b; PDGF, VEGF, FGF). These growth factors promote osteoblast proliferation and differentiation and angiogenesis and inhibit bone resorption (Gerstenfeld et al., 2003). We implemented three changes in the fuzzy rules between the early (stages 1 and 2) and late (stage 3) healing stages to reflect the influx of biological growth factors. During the early stage of healing the resorption rule was turned off so that the concentration of bone in any element was not allowed to decrease. This in effect removed the ‘‘resorption’’ region in Fig. 4. Studies have shown that TGF-b inhibits osteoclast formation, and therefore bone resorption (Chenu et al., 1988). Likewise, in trabecular bone fractures, resorption of the bone is not seen in the fracture gap before the gap is closed (Uhthoff and Rahn, 1981). After the fracture gap is bridged and the gap is stabilized, the concentration of the growth factors decreases, resorption and remodelling occurs, and a trabecular structure forms. At this time in the simulation, the resorption rule was activated so that areas with too little mechanical stimulation experienced bone resorption. Furthermore, to mathematically represent the influx of growth factors during healing, the rate of maximum bone apposition was 30% per iteration for the early stages of healing (until the fracture gap was bridged) and reduced to 10% per iteration during bone remodelling. Osteoblasts are found on the surface of existing bone, and therefore, bone formation occurs only where there is currently bone present. Very rarely are isolated ‘‘islands of bone’’ found in the fracture gap during healing. To reflect this in our model bone apposition was allowed only when the concentration of bone in a neighbouring element was above a threshold. In the early stage of healing the threshold was low (25%) to represent the increase in osteoblast activity due to growth factor influx. At the later stages of healing the threshold increased to 75%, and bone apposition only occurred when the current neighbouring bone concentration was high. All other fuzzy logic rules remained the same throughout all stages of the simulation. After the output from the fuzzy logic was determined, the material properties in the finite element model were changed accordingly and the next iteration began (Fig. 6). We continued the simulation until the structure did not change significantly for 50 consecutive iterations. As there was constant apposition and resorption, this homeostatic state was not static but was stable, and the overall trabecular structure remained the same.

Initial State material properties

vascularity

Finite Element Analysis s train s tate iteration

2446

Fuzzy Controller change in material properties

change in vascularity

Final State Fig. 6. Flowchart of simulation method. Simulation continued until a homeostatic structure was obtained.

3. Results The simulation predicted the stages of trabecular fracture healing: formation of woven bone, maturation of the woven bone, and remodelling of the woven bone to form trabeculae (Fig. 7). In the first stage of healing, woven bone formed in the fracture gap layer by layer starting on each side of the fracture gap and progressing toward the middle (iteration 0,6,11 in Fig. 7). For each iteration, the next layer of elements increased in concentration of bone, and therefore elastic modulus. The mechanical stimuli (hydrostatic and octahedral shear strain) in the gap were appropriate for bone formation throughout the entire gap. None of the elements had strains that led to cartilage or soft tissue formation. Bone formation was limited to elements that had neighbouring elements with the percentage of bone in the element over 25%. Therefore, islands of bone could not form in the middle of the gap but rather formed layer-by-layer. The modulus of the woven bone that formed was 250 MPa. After the bone formed, the strain in the new woven bone elements was too small to stimulate additional bone formation. After the fracture gap was filled with woven bone, the applied load was increased to simulate increased weight bearing on the healing bone (iterations 17–20). With increased load the strain was in the ‘‘intramembraneous bone’’ region of Fig. 4. The percentage of bone in each element increased resulting in an average elastic modulus of 2000 MPa in the elements in the fracture gap. The increase in elastic modulus represents the maturation of the woven bone. Because all elements had

ARTICLE IN PRESS S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

2447

Fig. 7. Cross-section through the middle of the model showing the change in the elastic modulus (E) over time. Iterations 1–16 represent the early stages with low loading and no bone resorption allowed. After the gap was bridged with woven bone the load increased in iterations 17–20. In subsequent iterations the load was maintained at 3 MPa and resorption was allowed. The woven bone was resorbed and formed trabecular structure.

over 25% bone already, the entire gap matured simultaneously and was not limited to layer-by-layer as in the first stage. The resorption of woven bone and formation of a trabecular structure occurred during the final stage (iterations 20–140 in Fig. 7). The elastic modulus first increased at the tips of the trabeculae on both sides of the fracture gap. The trabeculae then grew toward each other and progressively increased in elastic modulus. Simultaneously, between the trabeculae the mechanical stimulus was too low to maintain bone, bone was resorbed, and the elastic modulus decreased. The simulation continued until a stable structure (iteration 140) was maintained. This structure did not change significantly over 50 consecutive iterations and therefore was considered homeostatic. However, bone apposition and resorption continued in the elements at the surface of the trabeculae. Bone concentration and stiffness increased and decreased as the mechanical stimulus jumped back and forth between the ‘‘intramembraneous bone’’ and ‘‘resorption’’ regions of the tissue differentiation diagram (Fig. 4). This simulates the normal remodelling process in healthy bone. With a diagonal load applied during the late stage of healing the final trabecular structure formed diagonally across the cube (Fig. 8). The trabeculae that had no loading resorbed fully after iteration 40. After 120 iterations the structure reached a homeostatic state.

4. Discussion We have developed a model to simulate trabecular fracture healing, in particular the formation of woven bone and remodelling to form a trabecular structure. This model takes into account the influence of mechanical stimuli on the healing process, using finite element analysis and fuzzy logic. We used a micro-scale finite element model (mFE) with idealized geometry, material properties, and loading conditions in order to predict individual trabecular structure. Clearly in a clinical fracture or osteotomy, the trabeculae on either side of the gap would not be aligned in perfect columns. Subsequent models demonstrated that different initial trabecular structure bordering the fracture gap altered the final trabecular structure, but not the sequence of healing events (i.e. formation of woven bone and remodelling). Increasing the fracture gap size in the model also did not alter the final structure, but increased the number of iterations needed to reach a homeostatic state. Similar to previous models that predicted trabecular structure (Huiskes et al., 2000; Koontz et al., 2001; Ruimerman et al., 2003), loading conditions in our model were simplified. The application of pressure on the face of the trabeculae (vertical force) resulted in vertical trabeculae. With a diagonal load, trabeculae were built diagonally. Because no loads were applied in

ARTICLE IN PRESS 2448

S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

Fig. 8. With a diagonally applied load, the trabeculae in the final state formed diagonally. Only elements with E47000 MPa are illustrated.

the xy plane, no horizontal trabeculae were built. The loads in vivo on the trabeculae are multi-directional and changing, resulting in a much more complex trabecular structure. Additional models demonstrated that by changing the boundary conditions and loading parameters, we could alter resulting trabecular structure. The magnitude of the applied loads affected the width of the trabeculae that were built. With increased loading, the trabeculae were thicker. In the models shown, the loads were such that strains were within the ‘‘intramembraneous bone’’ and ‘‘resorption’’ regions of the tissue differentiation diagram (Fig. 4). Therefore, only rules in the fuzzy controller for intramembraneous ossification and bone resorption were activated. Additional models demonstrated that with increased loading cartilage formed in the fracture gap. With very high-loading tissue destruction occurred and all tissue became soft tissue, indicating that all rules in the controller could be activated with the appropriate conditions. The applied loads varied over time in three distinct stages (Fig. 3). The number of iterations for each stage was chosen such that the model reached a homeostatic state at the end of each stage before the new loading condition was applied. The number of iterations in each stage, however, could be increased without affecting the results. In particular, the load in stage 2 simulated increasing force on the fracture, representing increased weight bearing on the fractured limb and removal of the osteosynthesis. A smaller increase in load per iteration during stage 2 would not result in a different homeostatic bone structure.

The number of iterations required to reach the final structure was dependent on mesh density. With a finer mesh, it required more iterations to form the final structure. The time scale (number of days/iteration) was arbitrary, but it is important to note that it is dependent on the mesh size. The final trabecular structure was relatively independent of mesh density. Previous studies have shown that mFE models of trabecular structure should have 4–6 elements per trabecular cross-section in order to accurately predict cancellous bone apparent density (Homminga et al., 2001). In our model, the initial trabeculae thickness was four elements and the final trabeculae thickness was 4–5 elements. A finer mesh required significantly longer computation time and did not alter the final structure. Previous models of trabecular bone remodelling have demonstrated a dependence on mesh density (Jacobs et al., 1992; Weinans et al., 1992; Mullender et al., 1994; Weinans and Prendergast, 1996). Weinans et al. (1992) proved analytically that models using strain energy or stress as the stimulus for remodelling create a positive feedback system and produce mesh-dependent discontinuous solutions. The elements become either maximal or minimal density, creating a fractal pattern (Weinans and Prendergast, 1996). The fractal pattern has some resemblance to trabecular structure, but the discontinuity of the solution at the tissue level violates the continuum assumption on which finite element analysis is based. Subsequent models have proposed means to solve this problem in remodelling simulations. Jacobs et al. (1992) suggested using quadratic elements and averaging the solution over surrounding elements.

ARTICLE IN PRESS S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

Mullender et al. (1994) proposed separating the stimulus sensors from the mesh so that the solution was not mesh dependant. This method has been shown to produce trabecular structures in 2-D (Mullender and Huiskes, 1995; Huiskes et al., 2000), 3-D (Ruimerman et al., 2003), and predict bone loss in osteoporosis (Mullender et al., 1998). Because the stimulus in our method is strain instead of stress or energy, it does not exhibit the same positive feedback and discontinuous behaviour. Using strain as a stimulus is, in itself, a smoothing or averaging solution. Elements parallel to each other with different stiffnesses may have extremely different stresses, but will have relatively equal strains. We believe that the strain is a more physiological representation of the stimulus on the cell than stress or strain energy density. In vivo the cells are embedded in a matrix that is much stiffer than the cells themselves (Wang et al., 1993). Therefore the surrounding matrix has higher stress than the cells within the matrix. The strain in the cells and the matrix, however, are similar. The strain on the cell is a more accurate reflection of the mechanical environment surrounding the cell. Though the cell feels very little stress, the strain may be significant enough to stimulate a response. Additionally there is biological evidence of stretch activated ion channels in bone cells (Ypey et al., 1992). In this model we assumed all materials were linear elastic, homogeneous, and isotropic. Other tissue level models have used biphasic material properties and proposed that solid matrix shear stress and fluid flow over the cells stimulate tissue differentiation (Prendergast et al., 1997; Lacroix and Prendergast, 2002; Lacroix et al., 2002). By neglecting fluid flow in our model we may be underestimating the amount of shear strain the cells experience. We feel our assumption provided an appropriate first approximation of material properties, as well as conserving inordinate amounts of computer time. The boundaries between the regions in the tissue differentiation diagram (Fig. 4) certainly affect the final results. In particular in this model, the boundary between ‘‘resorption’’ and ‘‘intramembraneous bone’’ formation was critical. Previous studies have shown that strains under 100 mstrain leads to bone resorption (Rubin and Lanyon, 1987). In our model the ‘‘resorption’’ region had hydrostatic and octahedral shear strains below 100 mstrain. In vitro and in vivo studies have demonstrated that with strain values of 1000 mstrain bone apposition occurred (Rubin and Lanyon, 1987; Kaspar et al., 2000). The region between 100–1000 mstrain in our model was the fuzzy region, where either apposition or resorption could occur depending on the magnitude of the strain and other stimuli. The other boundaries in this model were not as critical as the applied loads did not result in strains in other regions. Future studies, however, will explore the sensitivity of theses boundaries.

2449

In this model all elements had an initial maximum vascularity to represent the rich vascular environment of the trabecular tissue. The mechanical conditions were such that the vascularity stayed maximal throughout the simulation. As angiogenesis is critical to the formation of bone, diminished vascularity could inhibit or prevent formation of bone. A previous model examined how diminished revascularization rates affect the healing of diaphyseal fractures, where initial vascular perfusion in the callus is low (Simon et al., 2003). This study indicated that the vascularity was critical in controlling the type of tissue formed. Without vascularity, cartilage formed instead of bone and drove the healing process toward endochondral ossification. Because angiogenic growth in this model is dependent on the mechanical environment, much higher strains in our simulation would reduce vascularity and inhibit bone formation. With a reduced initial vascularity or a decrease in vascularity due to mechanical conditions, the model presented in this study has the potential to predict the healing process in avascular or necrotic trabecular bone. Numerous models have been proposed to predict bone formation, density patterns, and remodelling. However, until now no single model was used to predict tissue differentiation in diaphyseal fracture healing, metaphyseal fracture healing, and bone remodelling. This study adapted a model that had previously been used to predict diaphyseal fracture healing (Simon et al., 2003) and predicted the sequence of events in metaphyseal fracture healing. A consistent model that explains bone formation, maintenance, healing, and degeneration on both a tissue and micro structural level could contribute towards understanding how processes during bone metabolism and tissue differentiation are influenced by vascular and mechanical environments.

Acknowledgements The authors would like to acknowledge the National Science Foundation International Research Post Doctoral Fellowship (INT-0202562) for funding this research. References Adachi, T., Tsubota, K., Tomita, Y., Hollister, S.J., 2001. Trabecular surface remodeling simulation for cancellous bone using microstructural voxel finite element models. Journal of Biomechanical Engineering 123 (5), 403–409. Ament, C., Hofer, E.P., 2000. A fuzzy logic model of fracture healing. Journal of Biomechanics 33 (8), 961–968. Augat, P., Merk, J., Ignatius, A., Margevicius, K., Bauer, G., Rosenbaum, D., Claes, L., 1996. Early, full weightbearing with flexible fixation delays fracture healing. Clinical Orthopaedics 328, 194–202.

ARTICLE IN PRESS 2450

S.J. Shefelbine et al. / Journal of Biomechanics 38 (2005) 2440–2450

Bailon-Plaza, A., van der Meulen, M., 2001. A mathematical framework to study the effects of growth factor influences on fracture healing. Journal of Theoretical Biology 212, 191–209. Carter, D.R., Hayes, W.C., 1977. The compressive behavior of bone as a two-phase porous structure. Journal of Bone and Joint Surgery— American Volume 59 (7), 954–962. Carter, D.R., Fyhrie, D.P., Whalen, R.T., 1987. Trabecular bone density and loading history: regulation of connective tissue biology by mechanical energy. Journal of Biomechanics 20 (8), 785–794. Carter, D.R., Blenman, P.R., Beaupre, G.S., 1988. Correlations between mechanical stress history and tissue differentiation in initial fracture healing. Journal of Orthopaedics Research 6 (5), 736–748. Chenu, C., Pfeilschifter, J., Mundy, G.R., Roodman, G.D., 1988. Transforming growth factor beta inhibits formation of osteoclastlike cells in long-term human marrow cultures. Proceedings of National Academy of Sciences of the United States of America 85 (15), 5683–5687. Claes, L.E., Heigele, C.A., 1999. Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing. Journal of Biomechanics 32 (3), 255–266. Cowin, S.C., Arramon, Y.P., Luo, G.M., Sadegh, A.M., 1993. Chaos in the discrete-time algorithm for bone-density remodeling rate equations. Journal of Biomechanics 26 (9), 1077–1089. Frost, H.M., 1992. Perspectives: bone’s mechanical usage windows. Bone and Mineral 19 (3), 257–271. Fyhrie, D.P., Carter, D.R., 1990. Femoral head apparent density distribution predicted from bone stresses. Journal of Biomechanics 23 (1), 1–10. Gerstenfeld, L.C., Cullinane, D.M., Barnes, G.L., Graves, D.T., Einhorn, T.A., 2003. Fracture healing as a post-natal developmental process: molecular, spatial, and temporal aspects of its regulation. Journal of Cellular Biochemistry 88 (5), 873–884. Hayes, W.C., Carter, D.R., 1976. Postyield behavior of subchondral trabecular bone. Journal of Biomedical Materials Research 10 (4), 537–544. Hildebrand, T., Ruegsegger, P., 1997. Quantification of bone microarchitecture with the structure model index. Computer Methods in Biomechanics and Biomedical Engineering 1 (1), 15–23. Homminga, J., Huiskes, R., Van Rietbergen, B., Ruegsegger, P., Weinans, H., 2001. Introduction and evaluation of a gray-value voxel conversion technique. Journal of Biomechanics 34 (4), 513–517. Hori, R.Y., Lewis, J.L., 1982. Mechanical properties of the fibrous tissue found at the bone–cement interface following total joint replacement. Journal of Biomedical Materials Research 16 (6), 911–927. Huiskes, R., Weinans, H., Grootenboer, H.J., Dalstra, M., Fudala, B., Slooff, T.J., 1987. Adaptive bone-remodeling theory applied to prosthetic-design analysis. Journal of Biomechanics 20 (11–12), 1135–1150. Huiskes, R., Ruimerman, R., van Lenthe, G.H., Janssen, J.D., 2000. Effects of mechanical forces on maintenance and adaptation of form in trabecular bone. Nature 405 (6787), 704–706. Jacobs, C., Levenston, M., Beaupre, G., Simo, J., 1992. A new implementation of finite element-based remodelling. Proceedings of the International Symposium on Computer Methods in Biomechanics and Biomedical Engineering. Swansea, Wales. Jarry, L., Uhthoff, H.K., 1971. Differences in healing of metaphyseal and diaphyseal fractures. Canadian Journal of Surgery 14 (2), 127–135. Kaspar, D., Seidl, W., Neidlinger-Wilke, C., Ignatius, A., Claes, L., 2000. Dynamic cell stretching increases human osteoblast proliferation and cicp synthesis but decreases osteocalcin synthesis and alkaline phosphatase activity. Journal of Biomechanics 33 (1), 45–51.

Koontz, J.T., Charras, G.T., Guldberg, R.E., 2001. A microstructural finite element simulation of mechanically induced bone formation. Journal of Biomechanical Engineering 123 (6), 607–612. Lacroix, D., Prendergast, P.J., 2002. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading. Journal of Biomechanics 35 (9), 1163–1171. Lacroix, D., Prendergast, P.J., Li, G., Marsh, D., 2002. Biomechanical model to simulate tissue differentiation and bone regeneration: application to fracture healing. Medical & Biological Engineering Computing 40 (1), 14–21. Mullender, M.G., Huiskes, R., 1995. Proposal for the regulatory mechanism of wolff’s law. Journal of Orthopaedics Research 13 (4), 503–512. Mullender, M.G., Huiskes, R., Weinans, H., 1994. A physiological approach to the simulation of bone remodeling as a selforganizational control process. Journal of Biomechanics 27 (11), 1389–1394. Mullender, M., van Rietbergen, B., Ruegsegger, P., Huiskes, R., 1998. Effect of mechanical set point of bone cells on mechanical control of trabecular bone architecture. Bone 22 (2), 125–131. Muller, R., Ruegsegger, P., 1995. Three-dimensional finite element modelling of non-invasively assessed trabecular bone structures. Medical Engineering & Physics 17 (2), 126–133. Pauwels, F., 1960. Eine neue theorie u¨ber den einfluX mechanischer reize auf die differenzierung der stu¨tzgewebe. Zeit & Chrift fur Antomie und Entwicklungsgeschichte 121, 478–515. Prendergast, P.J., Huiskes, R., Soballe, K., 1997. Esb research award 1996. Biophysical stimuli on cells during tissue differentiation at implant interfaces. Journal of Biomechanics 30 (6), 539–548. Rho, J.Y., Hobatho, M.C., Ashman, R.B., 1995. Relations of mechanical properties to density and ct numbers in human bone. Medical Engineering & Physics 17 (5), 347–355. Rubin, C.T., Lanyon, L.E., 1987. Kappa delta award paper. Osteoregulatory nature of mechanical stimuli: function as a determinant for adaptive remodeling in bone. Journal of Orthopaedics Research 5 (2), 300–310. Ruimerman, R., Van Rietbergen, B., Hilbers, P., Huiskes, R., 2003. A 3-dimensional computer model to simulate trabecular bone metabolism. Biorheology 40 (1–3), 315–320. Simon, U., Augat, P., Utz, M., Claes, L., 2003. Simulation of tissue development and vascularisation in the callus healing process. Transactions of 49th Annual Meeting Orthopaedic Research Society, New Orleans. Turner, C.H., Anne, V., Pidaparti, R.M., 1997. A uniform strain criterion for trabecular bone adaptation: do continuum-level strain gradients drive adaptation? Journal of Biomechanics 30 (6), 555–563. Uhthoff, H.K., Rahn, B.A., 1981. Healing patterns of metaphyseal fractures. Clinical Orthopaedics (160) 295–303. Van Rietbergen, B., Weinans, H., Huiskes, R., Odgaard, A., 1995. A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models. Journal of Biomechanics 28 (1), 69–81. Wang, N., Butler, J.P., Ingber, D.E., 1993. Mechanotransduction across the cell surface and through the cytoskeleton. Science 260 (5111), 1124–1127. Weinans, H., Prendergast, P.J., 1996. Tissue adaptation as a dynamical process far from equilibrium. Bone 19 (2), 143–149. Weinans, H., Huiskes, R., Grootenboer, H.J., 1992. The behavior of adaptive bone-remodeling simulation models. Journal of Biomechanics 25 (12), 1425–1441. Ypey, D.L., Weidema, A.F., Hold, K.M., Van der Laarse, A., Ravesloot, J.H., Van Der Plas, A., Nijweide, P.J., 1992. Voltage, calcium, and stretch activated ionic channels and intracellular calcium in bone cells. Journal of Bone and Mineral Research 7 (Suppl. 2), S377–S387.