Influence of fracture gap size on the pattern of long bone healing: a computational study

ARTICLE IN PRESS

Journal of Theoretical Biology 235 (2005) 105–119 www.elsevier.com/locate/yjtbi

Influence of fracture gap size on the pattern of long bone healing: a computational study M.J. Go´mez-Benitoa, J.M. Garcı´ a-Aznara, J.H. Kuiperb, M. Doblare´a, a

Group of Structural Mechanics and Material Modelling, Arago´n Institute of Engineering Research (13A), University of Zaragoza, Marı´a de Luna 3, 50015, Zaragoza, Spain b Unit for Joint Reconstruction, The Robert Jones and Agnes Hunt Orthopaedic Hospital, Oswestry, Shropshire SY10 7AG, UK Received 12 August 2004; received in revised form 22 December 2004; accepted 27 December 2004 Available online 22 February 2005

Abstract Following fractures, bones restore their original structural integrity through a complex process in which several cellular events are involved. Among other factors, this process is highly influenced by the mechanical environment of the fracture site. In this study, we present a mathematical model to simulate the effect of mechanical stimuli on most of the cellular processes that occur during fracture healing, namely proliferation, migration and differentiation. On the basis of these three processes, the model then simulates the evolution of geometry, distributions of cell types and elastic properties inside a healing fracture. The three processes were implemented in a Finite Element code as a combination of three coupled analysis stages: a biphasic, a diffusion and a thermoelastic step. We tested the mechano-biological regulatory model thus created by simulating the healing patterns of fractures with different gap sizes and different mechanical stimuli. The callus geometry, tissue differentiation patterns and fracture stiffness predicted by the model were similar to experimental observations for every analysed situation. r 2005 Elsevier Ltd. All rights reserved. Keywords: Bone healing; Fracture callus; Mechano-biology; Gap size

1. Introduction Fracture healing is a natural process that regenerates bone to its original state and function. It is a complex process, involving proliferation, migration, and/or differentiation of mesenchymal stem cells (MSCs), fibroblasts, cartilage cells and bone cells, which synthesize different tissues such as granulation tissue, fibrous tissue, cartilage and bone (Einhorn, 1995). Several factors influence these bone healing events, such as genetic, cellular and biochemical factors, age, the type of fracture, interfragmentary motion and fracture Corresponding author. Tel.: +34 976761000X5111; fax: +34 976762578. E-mail addresses: [email protected] (M.J. Go´mez-Benito), [email protected] (J.M. Garcı´ a-Aznar), [email protected] (J.H. Kuiper), [email protected] (M. Doblare´).

0022-5193/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2004.12.023

geometry (Hadjiargyrou et al., 1998; Goodship et al., 1993). Much interest in this research field has focussed on mechanical factors, in particular movement at the fracture site (Claes et al., 1995; Kenwright and Gardner, 1998; Gardner and Mishra, 2000; Claes et al., 1997). This interest follows a long-standing debate in orthopaedics about motion versus rest as best treatment for fractures (Connolly, 2001; Sarmiento and Latta, 1995). Not only movement at the fracture is important but also the geometry of the fracture. One of these geometrical aspects is the size of the fracture gap, that together, with the rest mechanical conditions, will determine the strain in the gap newly formed tissue. A single experimental study exists which addresses the combined effects of gap size and applied strain (Claes et al., 1994, 1997; Augat et al., 1998). These authors used a sheep model to compare healing with a small, medium or large gap size, subjected to small or large

ARTICLE IN PRESS 106

M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

interfragmentary motions. Compared to fractures with a small gap, those with a medium gap produced a larger callus but had lower fracture stiffness, regardless the applied movement. Fractures with a large gap and large motion produced little callus and had low stiffness (Claes et al., 1997). The authors concluded that fracture gap size influenced the healing process more than interfragmentary motion. The above animal study was limited to a simple fracture geometry and few strain values. Applying the above findings to more realistic cases is difficult. Moreover, the proper understanding of the interaction between fracture gap and interfragmentary movement would help those attempts to grow new bone tissue in bone defects (Kelly and Prendergast, 2004). Mathematical models of the fracture healing process could be therefore an aid in these attempts. Several such models have been proposed (Prendergast et al., 1997; Lacroix and Prendergast, 2002; Kuiper et al., 2000; Carter et al., 1998; Claes and Heigele, 1999; Gardner et al., 2000; Bailo´n-Plaza and van der Meulen, 2003; Ito et al., 2004), all restricted to a fixed geometry of the healing callus. Since the experimental results clearly show that different gap sizes and mechanical conditions lead to different geometries of the callus, this restriction limits the value of these models. The callus geometry must be treated as other variable of the model to analyse the combined effect of fracture geometry and applied movement. Recently, we have formulated a mathematical model that uses callus geometry as a dependent variable (Garcı´ a-Aznar et al., 2002). This model uses mechanobiological rules to describe the influence of mechanical stimuli and time on cell proliferation and differentiation. Based on proliferation and differentiation, it determines callus growth and callus geometry. Whether such a model can correctly predict the influence of fracture gap and strain on the fracture healing process is not known. Some of the rules in the model were based on poorly tested assumptions. The most important of those was the rule controlling differentiation of mesenchymal stem cells into fibroblasts. Contrary to differentiation of stem cells into bone cells or cartilage cells, which is easily verified histologically, differentiation of stem cells into fibroblasts is more difficult to verify. The two cell types have a similar appearance: often stem cells in fractures are described as ‘‘fibroblast-like’’ (Bland et al., 1999; Postacchini et al., 1995). In the absence of direct histological comparison, the sensitivity of callus geometry and fracture stiffness to the mechano-biological rules can be used to test their importance. The first aim of the present study was therefore to test the hypothesis that this model can simulate the findings from the above experimental study. A second aim of the present study was to use the results from the above experiments to test its initial assumptions and find the most appropriate formulations.

2. Material and methods 2.1. Mechano-biological model description The main variables assumed were the cell concentrations of the four basic skeletal cells: MSCs ðcs Þ; fibroblasts ðcf Þ; cartilage cells ðcc Þ and bone cells ðcb Þ (Fig. 1). These four types produce the various skeletal tissues, namely granulation tissue, fibrous tissue, cartilage and bone. The percentage of these basic types was assumed to determine the mechanical properties of the local tissue. We assumed that the fracture site is invaded by proliferating and migrating MSCs in the first stage of the bone healing process. These cells may differentiate into cartilage cells, bone cells or fibroblasts, depending on the value of a mechanical stimulus. In addition, we assumed that a very high mechanical stimulus will cause MSCs to die. A further assumption was that MSC concentration ðcs Þ can increase until a saturation level, csmax ; which cannot be exceeded. Bone cells ðcb Þ were assumed to appear through either of two different ossification processes: intramembranous or endochondral. In the first process, bone cells differentiate directly from MSCs. In the second, mesenchymal cells first differentiate into chondrocytes and produce cartilage (chondrogenesis). Later, this cartilage is replaced by bone (endochondral ossification). Both processes were assumed to be guided by a mechanical stimulus. During healing, callus size and shape are constantly modified under the influence of several factors. Different interfragmentary movements have been observed to produce different callus sizes, suggesting mechanical stimulation is a particularly determinant factor. The model assumed two callus growth mechanisms: MSC proliferation and chondrogenesis. The callus stiffness was characterized through the volume of extracellular matrix for each tissue type and through its geometry (see Fig. 1). We assumed that the continuity Eq. (1), as stated by Garcı´ a-Aznar et al. (2002), governs the evolution of skeletal cells:
 DN Dcðx; tÞ ¼ þ cðx; tÞ divðvÞ  dV , Dt Dt

(1)

where N is the number of cells in a control volume dV ðN ¼ cdV Þ; DN=Dt the rate of change in the number of cells, t time (day), cðx; tÞ cell density (number of cells=mm3 ) and v growth rate (mm/day). This general expression quantifies the evolution of the number of cells in a control volume ðdV Þ: The number of cells can be modified through a change in cell concentration ðDcðx; tÞ=DtÞ or volume growth rate ðdivðvÞÞ: Following (Garcı´ a-Aznar et al., 2002) we also

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

107

Fig. 1. Scheme of the main events in the mechano-biological model of bone fracture healing, guided by the mechanical stimulus.

assumed that each cellular event during bone fracture healing can be described through Eq. (1). The process of bone healing was simulated as a process driven by a mechanical stimulus. We chose the second invariant of the deviatoric strain tensor ðJ 2 Þ as the stimulus responsible for the process. This stimulus is a function of location x and time t: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cðx; tÞ ¼ J 2 ¼ ðI  oct Þ2 þ ðII  oct Þ2 þ ðIII  oct Þ2 , (2) where I ; II and III are principal strains, and oct ¼ ðI þ II þ III Þ=3 is the octahedral strain. In the next two sections, we formulate descriptions for the evolution of cell concentration and growth for each cell type. 2.2. Cell proliferation and differentiation MSC concentration was assumed to change by proliferation, migration and differentiation (cell death was considered as a specific differentiation pathway). The evolution of stem cell concentration due to these processes was described by Dcs ðx; tÞ ¼ f proliferation ðcs ; cÞ þ f migration ðcs ; V disrupted Þ Dt  f differentiation ðc; tÞ aproliferation  cðx; tÞ  cs ¼ cðx; tÞ þ cproliferation  DðV disrupted Þ  r2 cs  f differentiation ðc; tÞ,

ð3Þ

where cs is MSC concentration (number of cells=mm3 ), V disrupted fraction of disrupted tissue, c mechanical stimulus, aproliferation ; cproliferation are constants that define stem cell proliferation (see Appendix A) and DðV disrupted Þ

the diffusion coefficient in mm2 =day: Stem cell concentration was assumed to be limited to a maximum csmax : The first term in the right-hand side of expression (3) simulates stem cell proliferation as function of mechanical stimulus. Stem cells proliferate by mitosis, so this effect was assumed proportional to the current stem cell concentration. Proliferation was also assumed to depend on the mechanical stimulus, such that it increases with stimulus until it reaches a maximum ðcdeath Þ for which cells were assumed to die. The assumption that mechanical strain increases proliferation was based on the observations by Smith and Roberts (1980), who measured a marked proliferation increase of precursor cells within 6 h following strain application. The second term of the right-hand side of (3) models random migration of mesenchymal cells through the matrix. We assumed that cells would migrate slower in disrupted tissue and modelled this behaviour by making the migration coefficient dependent on the volume fraction of disrupted tissue: DðV disrupted Þ ¼ D0 

1  V disrupted , V disrupted þ 0:01

(4)

where DðV disrupted Þ is the diffusion coefficient in mm2 =day; V disrupted the fraction of disrupted tissue, and D0 a constant defined in Appendix A. D tends to zero for completely disrupted tissue ðV disrupted ¼ 1Þ: This assumption is based on the argument that cell movement requires matrix attachment (Horwithz and Parsons, 1999), something which is very difficult in severely disrupted tissue. Directed migration, e.g. under the influence of growth factors or cytokines, was neglected. The final term that was assumed to control MSC dynamics is cell differentiation. Based on experiments, it was suggested recently that mesenchymal cell differentiation is time-dependent as well as mechano-biologically dependent (Cullinane et al., 2003). Following this

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

108

idea, we propose a set of rules that controls the process of MSC differentiation into specialized cells (such as bone cells, cartilage cells and fibroblasts) as a function of mechanical stimulus and time: f differentiation ðc; tm Þ

¼

8 hintramembranous ðc; tÞ > > > > > > > > > > > > gdifferentiation ðc; tÞ > > > > > <

if ðclim ococbone Þ

> l differentiation ðc; tÞ > > > > > > > > > > > > > cs > > > : 0

if ðccartilage ococfibrous Þ

and ðt4tbm Þ ðbone cellsÞ; if ðcbone ococcartilage Þ and ðt4tcm Þ ðcartilage cellsÞ; and ðt4tfm Þ ðfibroblastsÞ; if ðcdeath ocÞ ðdeath cellsÞ;

dent of the mechanical stimulus: tcm ¼ tcmlim

ð5Þ

with hintramembranous ðc; tÞ; gdifferentiation ðc; tÞ; l differentiation ðc; tÞ the functions that define the evolution to osteoblasts, chondrocytes and fibroblasts, respectively (number of cells=mm3 day), tim maturation time needed for each cell type i to mature, and clim ; cbone ; ccartilage ; cfibrous mechanical stimulus limits for each cell type (see Appendix A). The dependency of these rules on time is defined by the maturation time ðtim Þ; the time that stem cells need to differentiate into specialized cells ðiÞ: We assumed that this period is different for each cell type and depends on mechanical factors for differentiation to bone cells and fibroblasts. The notion of time-dependency follows Cullinane et al. (2003). Based on the time course of proliferation and gene expression in cell cultures before bone matrix formation, the sequence of osteoblast development from precursor cells has been outlined as a sequence of proliferation, maturation (i.e. differentiation) and mineralization (Aubin et al., 1995). This sequence was linked recently to a specific transcription factor (Pratap et al., 2003). We therefore assumed that mechanical stimuli that increase proliferation will delay maturation. As a first approach we have simulated this effect by a linear function: tbm ¼ tcmlim

c cbone

if cocbone ,

(6)

where tcmlim is the time limit defined in Appendix A. In case of large fracture site movements, the principal mechanism of regeneration observed is the formation of cartilaginous tissue that is replaced gradually by new bone (endochondral ossification). The cartilage cells producing cartilage differentiate from MSCs. Most of them appear simultaneously 15–20 days after fracture ðtcmlim Þ in humans (Cullinane et al., 2003) and by day 7 in rats (Iwaki et al., 1997). We neglected any difference in maturation time between cartilage cells and assumed that cartilage maturation time is constant and indepen-

(7)

There not exist enough experimental data to determine the time MSCs need to differentiate into fibroblasts. In order to study the influence of mechano-regulatory influences on fibroblast differentiation, three different functions were analysed: a linear function in which maturation time increases with mechanical stimulus (8a), a positive bilinear law in which maturation time increases linearly with mechanical stimulus above threshold stimulus cslope (8b) and a negative bilinear one, in which maturation time decreases with mechanical stimulus up to a threshold limit cslope (8c). tfm ¼

other cases;

if cbone ococcartilage .

c  ccartilage tslope þ tcmlim cfibrous  ccartilage mlim if ccartilage ococfibrous ,

8 c > < tmlim c  cslope tfm ¼ c > tslope :c mlim þ tmlim fibrous  cslope

ð8aÞ if ccartilage ococslope ; if cslope ococfibrous ;

(8b)

tfm

8 > <

c  cslope c tslope mlim þ tmlim ¼ ccartilage  cslope > : tc mlim

if ccartilage ococslope ; if cslope ococfibrous ;

(8c) where cslope is the level of stimulus which represents the change in the slope of the bilinear law proposed, tslope mlim is a time limit defined in Appendix A. Having defined the maturation time for each cell type, the remaining differentiation functions in Eq. (5) can be defined. Bone cells are known to appear during intramembranous and endochondral ossification, through differentiation from MSCs or following calcification of cartilage, respectively. Migration and proliferation of osteoblasts were considered negligible in comparison with the increase of osteoblast population due to intramembranous and endochondral ossification. We assumed that the evolution of bone cell density ðcb Þ depends on the mechanical stimulus c and maturation time ðtm Þ through the functions hintramembranous and hendochondral ; which distinguish between intramembranous and endochondral ossification: Dcb ¼ hintramembranous ðc; tÞ þ hendochondral ðc; tÞ. (9) Dt In both processes, penetration of new blood vessels is a precursor to bone formation (Glowacki, 1998). New vessels in fractures grow out from existing capillaries, which are found in the periosteum and newly formed bone (Glowacki, 1998). We chose not to model this

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

vascular ingrowth directly, but incorporated it as diffusion of osteoblasts. This recognizes the close association between ingrowth of vessels and appearance of bone cells. We further assumed that a threshold level of bone cells ðcmin b Þ would be indicative of sufficient vascularization, in which case MSCs could differentiate to bone cells directly. Both types of ossification have been modelled according to this idea. Focussing first on the intramembranous ossification, we propose: hintramembranous ðc; tÞ 8 Db  r2 cb if ðclim ococbone Þ; ðt4tbm Þ > > > > > > and ðcb ocmin > b Þ; > < if ðclim ococbone Þ; ðt4tbm Þ ¼ cs > > > > and ðcb 4cmin > b Þ; > > > : 0 other cases;

ð10Þ

with Db diffusion coefficient ðmm2 =dayÞ; cb bone cell density (number of cells=mm3 ), clim and cbone stimulus levels that control the ossification process (see Appendix A), cmin minimum bone cell concentration that controls b the rate of intramembranous ossification process, tsm maturation time that mesenchymal stem cells need to differentiate into osteoblasts. Endochondral ossification has been modelled analogously through: hendochondral ðc; tÞ 8 DðcÞ  r2 cb > > > > > > > > < ¼ cc > > > > > > > > : 0

Evolution of chondrocytes was defined by the following expression: Dcc ðx; tÞ ¼ gdifferentiation ðc; tÞ þ gendochondral ðc; tÞ, Dt

(13)

where gdifferentiation is the function that controls MSC differentiation into chondrocytes (number of chondrocytes=mm3 day) and has been defined as: gdifferentiation ðc; tÞ ( cs if ðcbone ococcartilage Þ and ¼ 0 other cases;

ðt4tcm Þ;

ð14Þ

gendochondral ¼ hendochondral is the change of chondrocyte density due to endochondral ossification (number of chondrocytes=mm3 day) and hendochondral was defined in (11). MSCs can also differentiate into fibroblasts. Migration and proliferation of fibroblasts were assumed negligible. We thus assumed: Dcf ðx; tÞ ¼ l differentiation ðc; tÞ, Dt

(15)

where l differentiation ðc; tÞ ( cs if ðccartilage ococfibrous Þ ¼ 0 other cases;

and

ðt4tfm Þ;

ð16Þ

with cs MSC concentration (number of cells=mm3 ) and ccartilage and cfibrous limits corresponding to fibroblast differentiation (see Appendix A).

if ðcocbone Þ; ðpmi 4pmin mi Þ and ðcb ocmin b Þ; if ðcocbone Þ; ðpmi 4pmin mi Þ

ð11Þ

and ðcb 4cmin b Þ; other cases

ðcbone  cÞ , 2  cbone

2.3. Callus growth model We assumed that callus growth was mainly due to mesenchymal cell proliferation and chondrocyte hypertrophy during endochondral ossification:

with cb bone cell density (number of cells=mm3 ), cc cartilage cells concentration (number of cells=mm3 ), cbone stimulus value that limits ossification, pmi percentage of mineralization of cartilage, pmin and cmin mi b constants defined in Appendix A. Compared to intramembranous ossification, we added the assumption that the diffusion coefficient decreased with increasing mechanical stimulus DðcÞ: DðcÞ ¼ Db

109

(12)

where Db =2 is the maximum diffusion coefficient. Hypertrophy of chondrocytes causes a cartilage volume increase. After cartilage calcification, the cells undergo programmed cell death (apoptosis), leaving the matrix for the invasion of blood vessels and consequently osteoclasts and osteoblasts. Migration and proliferation of chondrocytes were assumed negligible.

divðvÞ ¼ f vproliferation ðcs ; cÞ þ gvendochondral ðc; tÞ,

(17)

where f vproliferation ðcs ; cÞ defines the rate of callus growth due to proliferation, v is the growth rate (mm/day) and gvendochondral ðc; tÞ controls the rate of callus growth due to chondrocyte hypertrophy. It was assumed that the concentration of mesenchymal cells can vary between zero and a maximum or saturation cell density. When the saturation concentration of MSCs is reached the only way cells can proliferate further is by increasing the callus size at a constant level of cell concentration: f vproliferation 8 0 > < ¼ f proliferation ðcÞ > : csmax

if ðcs ocsmax Þ; if ðcs ¼ csmax Þ;

ð18Þ

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

110

Table 1 Composition used for each different tissue type (Martin et al., 1998) Volume

Debris tissue

Bone tissue

Cartilage

Calcified cartilage

Fibrous tissue

Granulation tissue

Collagen I (cI) Collagen II (cII) Collagen III (cIII) Ground substance (gs) Mineral (mi)

0 0 0.018 0.082 0

0.2848 0 0 0.0352 0.43

0 0.135 0 0.079 0

0 0.135 0 0.079 0.015

0.1861 0 0 0.07885 0

0.018 0 0 0.082 0

where f proliferation is the proliferation function of MSCs (number of MSCs=day mm3 ) defined in (3), csmax maximum MSC density (number of MSCs=mm3 ) and cðx; tÞ mechanical stimulus. During endochondral ossification, cartilage cells swell up to ten times their original volume (Breur et al., 1991). The number of cartilage cells remains almost constant during the process (Wilsman et al., 1996), which allows to define growth from cartilage hypertrophy: 1 qcc ðx; tÞ  cc ðx; tÞ qt ggrowth ðc; tÞ ¼  , cc ðx; tÞ

gvendochondral ðc; tÞ ¼ 

where ggrowth is 8 c  ccalcified > > khyper > < cbone  ccalcified ggrowth ¼ > > > : 0

ð19Þ

if ðcoccalcified Þ and ðcc occmin Þ; other cases; (20)

with cc chondrocyte density (number of chondrocytes= mm3 ), khyper ; ccalcified and ccmin constants defined in Appendix A. It was assumed that the chondrocyte density will reduce until an equilibrium value ðccmin Þ: At this point cell volume has reached a maximum and further growth ceases. 2.4. Extracellular matrix evolution At the healing site, five main tissues are found: debris tissue with a low stiffness (occupying the gap site after fracture), granulation tissue where MSCs are found, cartilage matrix, fibrous tissue and bone tissue. Properties of each of these tissues can vary. For example, properties of cartilage vary from newly formed to calcified cartilage. Two different types of bone can be found: immature and mature bone. The first one is a disorganized tissue, which is replaced by mature bone when the fracture is more stable. The matrix was characterized by its density and composition. The composition of each separate tissue was assumed fixed and defined in Table 1. The main components of the matrix were assumed to be water,

mineral, ground substance and collagens types I, II and III. The production rate of extracellular matrix volume was assumed to depend on cell type, cell density and matrix production rate per cell: qV imatrix ¼ ci  Qi , (21) qt where V imatrix is the volume fraction of tissue i, i.e. granulation tissue (s), fibrous tissue (f) or cartilage (c) in the extracellular matrix, ci is the cell density, Qi is the matrix production rate per cell (defined in Appendix A). This expression was used for each tissue type except mature bone. The production of bone matrix volume in mature bone was analysed using the internal bone remodelling formulation proposed by Beaupre´ et al. (1990a, b). With this, the expression for mature bone matrix volume evolution became: qV bmatrix ¼ krem  r_  Sv if ð0ococmature Þ, (22) qt where r_ is the formation/resorption of bone matrix volume per available bone surface per unit time, S v is the specific bone surface (Martin, 1989) and krem is the percentage of bone surface active for remodelling. All tissues were assumed biphasic, isotropic and linear elastic. The material properties of the different tissues, except for mature bone, were computed by using a mixture rule, similar to other authors (Ament and Hofer, 2000; Lacroix, 2001). The modulus of elasticity E and the Poisson’s ratio n were determined from the proportion of each component px (the subscript x denotes the specific component, see Table 1): EðMPaÞ ¼ 2000pmi þ 430pcI þ 200pcII þ 100pcIII þ 0:7pgs , n ¼ 0:33pmi þ 0:48ðpcI þ pcII þ pcIII Þ þ 0:49pgs .

(23)

Elastic properties of mature bone were determined as a function of apparent density (Jacobs, 1994): E ¼ 2014  r2:5 ; 3:2

E ¼ 1763  r ;

n ¼ 0:2 n ¼ 0:32

if ðrp1:2g=cc:Þ, if ðrX1:2g=cc:Þ.

ð24Þ

Permeability depends on type of tissue. In soft tissues, the permeability was computed from composition and porosity, following the model proposed by Levick (1987). However, for hard tissues we assumed a constant

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

permeability, distinguishing between cortical bone, trabecular bone, and calcified cartilage, with respective values of 7  1015 ; 7  1011 and 3:5  1012 mm2 : 2.5. Disrupted soft tissue A bone fracture causes local haemorrhage and tissue destruction. This damaged tissue is known as debris tissue and its spatial distribution is fully dependent on how the fracture was produced. We assumed that the only mechanism through which debris tissue appears in the fracture gap is the one described in the next subsection. Within days after fracture the inflammatory process begins to remove debris from the fracture site and form a hematoma. At the same time, MSCs proliferate and capillary sprouts grow into the blood clot in the injured area, forming granulation tissue, in response to cytokines released by tissue damage. This newly formed tissue allows the invasion and proliferation of MSCs to and at the fracture site. If the fracture is not sufficiently stabilized and excessive movement continues, additional disruption of this granulation tissue can occur (Gardner et al., 2000). This can alter the normal reparative process. We incorporated this effect in the proposed model, through the variable V disrupted : This variable quantifies the fraction of matrix disrupted with respect to the reference volume. It varies between [0,1[ where 0 represents an intact matrix and 1 a completely disrupted matrix. We assumed that granulation tissue is locally disrupted when the mechanical stimulus in the granulation tissue reaches a limit stimulus cdisrupted : Following such an event, the disrupted volume fraction V disrupted ; was assumed to be 0.95. Matrix disruption was assumed to reduce elastic modulus of the granulation tissue: E new ¼ E 0 ð1  V disrupted Þ,

(25)

111

where E 0 is the elastic modulus of the intact granulation tissue. This disrupted tissue was supposed to be uniformly distributed in the matrix and to be repaired with new matrix once the stimulus level had reduced sufficiently. We assume that repair would be downregulated when it is near the completion (see also Olsen et al., 1995). This was modelled using a linear relation between disrupted tissue fraction and repair rate. MSCs were therefore assumed to produce new matrix according to: s V_ disrupted ¼ V disrupted  V_ matrix

(26)

with V_ disrupted the rate of change of disrupted matrix s fraction and V_ matrix the matrix volume production rate of granulation tissue. 2.6. Numerical implementation The Finite Element Method (ABAQUS v.6.3 Hibbit, Karlsson and Sorensen, Inc., 1999) was used to solve all the previous equations. State variables (cell concentration cx for every cell type and volume fraction V x of each matrix component) were stored as nodal values. The callus geometry was defined by nodal positions (Fig. 2). The geometry of a simplified human tibia was considered. A cloud of nodes was generated to define the tibia. A mesh was produced from this cloud of nodes using the automatic mesh generator Qhull (Barber et al., 1996), yielding a two-dimensional mesh of triangles. Then, three different analysis steps were performed in order to determine the evolution of the state variables and the callus shape and size. In the first analysis stage, the mechanical stimulus in the matrix of different tissues was determined using a poroelastic analysis. Next, a diffusion analysis step was performed to simulate the migration of MSCs and the advance of the ossification front. Finally, a thermoelastic analysis stage was

Fig. 2. Flow chart of the iterative computational implementation.

ARTICLE IN PRESS 112

M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

performed to determine the new position of the nodes and the new callus geometry. This step uses the conceptual similarity between growth and thermal expansion of solids. In this context, the reciprocal of cell density (MSCs or chondrocytes) was used as thermal expansion coefficient, and the values of the proliferation function f proliferation ðcÞ or the chondrocyte growth function ggrowth ðc; tÞ were used as increments of temperature. The new callus geometry was defined through the displacements obtained in this thermo elastic analysis. After the three analysis steps, all nodal state variables were updated as a function of the mechanical stimulus determined in the poroelastic analysis, taking into account cell differentiation, proliferation and death. Once the updated concentration of every cell type was determined, the change of matrix volume was computed. With this, the loop was completed and a new time increment consisting of the three analysis steps was performed (Figs. 1 and 2). Before the start of the increment the automatic mesh generator was used to make a new mesh based on the new position of the nodes. 2.7. FE model description A simplified axisymmetric model of a human tibia was developed, with a width for periosteum, cortex, endosteum and bone marrow of 1, 7, 0.5 and 6.5 mm, respectively. Three values for the fracture gap size were analysed: 1, 2 and 6 mm. Symmetry about the fracture line was assumed (Fig. 3a). The definition of all the functions and parameters used are shown in Appendix A. Identical initial conditions were imposed to the three gap sizes

Fig. 4. Interfragmentary movement as function of the healing time (Claes et al., 1997).

(Fig. 3b): MSCs were assumed to exist in the periosteum, and a line of preosteoblastic cells was assumed on the layer between periosteum and cortical bone. It was also assumed that periosteum was removed at the fracture gap, and surrounding soft tissue had a poor vascularization. The gap was initially filled with debris tissue. Clearly, the interfragmentary movement at the fracture site depends on the stiffness of the fixator and the gap size. In their animal study, Claes et al. (1997) chose fixator stiffnesses such that an initial interfragmentary strain of 31% was achieved for gap sizes of 1, 2 and 6 mm. For all the animals they measured the interfragmentary movements along time (Fig. 4). These movements were used as applied displacement at the bone end in our model. 2.8. Analyses performed

Fig. 3. Boundary conditions of the poroelastic analysis (a) and initial distribution of mesenchymal stem cells (MSCs) (b).

Firstly, we assumed that fibroblast differentiation was controlled by Eq. (8c). Under this condition, the influence of the gap size was analysed. Secondly, we studied the influence of the fibroblast differentiation rule ((8a) and (8b)), the influence of diffusion vs. direct transformation during intramembranous and endochondral ossification ((10) and (11)) and the influence of the mechanical stimulus on endochondral ossification (12). To validate the predictions, we compared them with the experimental work in sheeps accomplished by Claes and coworkers (1997). Although these animal experiments were not exactly identical to our simulations, they are indicative of the influence of gap size on the fracture healing course.

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

3. Results 3.1. Influence of the gap size During the first two weeks of simulated healing, MSCs proliferated in the external periosteum causing the callus to grow (compare Fig. 5 and Fig. 3b). In the same time, they also migrated, filling the periosteal callus. Some stem cells differentiated into chondrocytes (Fig. 6). Intramembranous ossification occurred by differentiation of MSCs into osteoblasts, visible as an ossification front advancing from the periosteum (Fig. 7). Callus size and shape at day 14 strongly depended on the gap size, with the largest callus appearing in the case of a 2 mm gap, a smaller in the

113

case of 1 mm gap and the smallest callus for a gap of 6 mm (Fig. 5). Inside the gap, no MSCs were present due to the large interfragmentary strains that provoke MSCs to die in this region (Fig. 5). During week 3 and 4, MSCs differentiated mainly into cartilage cells (Fig. 6) but also into fibroblasts (Fig. 8). Fibroblasts appeared in the gap near the periosteum in the area of highest stimulus. The ossification front advanced further, due to the combination of intramembranous and endochondral ossification. At day 28, the only remaining location of MSCs was the periosteum. They had disappeared from the gap (Fig. 5). Cartilage cells started to appear 14 days after fracture (Fig. 6) through direct differentiation from MSCs (Fig. 5) for the three gap sizes. Comparing cartilage cell

Fig. 5. Evolution of MSCs.

Fig. 6. Evolution of cartilage cells.

ARTICLE IN PRESS 114

M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

Fig. 7. Evolution of bone cells.

Fig. 8. Evolution of fibroblasts.

concentration between 14 and 28 days, these cells occupied the complete external callus after 4 weeks in the 1 and 2 mm gaps, but did not form in the 6 mm gap, where most of the callus was filled with fibroblasts (Fig. 8). From day 14 onwards, endochondral ossification took place. The cartilage region transformed into bone, in the two cases with smallest gap size (Figs. 6 and 7). Swelling of cartilage cells during this process gradually modified the size and external shape of the callus until osseous bridging occurred between the two fracture fragments (Fig. 9). This feature was less prominent in the 1 mm gap, where movement and callus were smaller. Finally, at day 56, the callus was completely occupied by bone cells (Fig. 7), although some small islands of cartilage remained in the gap of 2 mm (see Fig. 6).

Fig. 9. Axial stiffness throughout healing for the three gap sizes.

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

Endochondral ossification did not occur in the 6 mm gap, where the cartilage cell concentration in the gap was very small and almost constant during the whole process (Fig. 6). Since the cartilage cells did not swell the callus size remained constant. Intramembranous ossification occurred in the periosteum, the sole bone formation mechanism in this case. Little bone was formed in the callus at the end of the simulation (Fig. 7). Fibroblasts occupied most of the callus (Fig. 8). During the first days of the simulated healing, the axial fracture stiffness was very similar in all cases (Fig. 9). However, later, fracture stiffness showed a large influence of gap size. After 25 days, a high rise of stiffness occurred in the fractures with 1 and 2 mm gap (Fig. 9). The fracture with the 2 mm gap remained less stiff than that with the 1 mm gap up to day 56, when its stiffness was still 40% smaller. The fracture with largest gap (6 mm) retained a low stiffness throughout the simulation (Fig. 9).

115

3.3. Influence of the intramembranous ossification rule Intramembranous ossification was simulated by a combination of diffusion and differentiation (10). With this expression, the evolution of bone cell concentration four days after fracture was slower than that obtained by simulating intramembranous ossification as a direct differentiation from MSCs (Fig. 11) (Garcı´ a-Aznar et al., 2002). At day 56, the end of the simulation, the fracture healing course, callus shape, size and stiffness predicted by the two ways of simulating intramembranous ossification were indistinguishable for the three gap sizes. 3.4. Influence of the endochondral ossification rule Endochondral ossification was simulated by a combination of diffusion and replacement of cartilage by bone (11). With this rule, the ossification front was observed to advance continuously (Fig. 12a). Modelling

3.2. Influence of the fibroblast differentiation rule To decide which rule best describes fibroblast differentiation (8), its influence on the healing pattern was analysed for the three different gap sizes (1, 2 and 6 mm). Use of the linear law (8a) predicted a low stiffness for the 6 mm gap and a larger stiffness for the 1 and 2 mm gaps (Fig. 10a). This law did not allow fibroblasts to appear and the three cases healed successfully. In the case of a 6 mm gap the final shape of the callus was very irregular. Use of the positive bilinear law (Fig. 10a) predicted the fracture with 2 mm gap never to heal and fibrous tissue to appear in nearly the whole callus area. A callus with low stiffness was formed. Use of the negative bilinear law (8c) predicted the final stiffness of the callus to decrease with increasing gap size.

Fig. 10. Axial stiffness for the three gap sizes and the three fibroblast differentiation laws proposed.

Fig. 11. Bone cell concentration during intramembranous ossification (a) when assuming bone cell migration precedes MSC differentiation and (b) when assuming bone cell migration does not play any role.

Fig. 12. Bone cell concentration during endochondral ossification (a) when assuming bone cell migration precedes cartilage differentiation and (b) when assuming bone cell migration does not play any role.

ARTICLE IN PRESS 116

M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

Fig. 13. Influence of the mechanical stimulus in the diffusion of osteoblasts (a) diffusion coefficient as a decreasing function of the mechanical stimulus and (b) constant diffusion coefficient.

endochondral ossification as a single process of direct replacement of cartilage cells by bone cells caused a modification of the fracture healing course (Fig. 12b) (Garcı´ a-Aznar et al., 2002). In this case, bone cells appeared in the callus far from the ossification front. 3.5. Influence of the mechanical stimulus on endochondral ossification When proposing the diffusion coefficient for osteoblasts during endochondral ossification to be a function of the mechanical stimulus (12), an osseous bridge was achieved before osteoblasts entered the gap. With a constant diffusion coefficient, osteoblasts entered the gap before completion of the osseous bridge (Fig. 13).

4. Discussion The main purpose of this work was to test the hypothesis that a cell-based mechano-biologically regulated mathematical model of skeletal tissue growth and differentiation (Garcı´ a-Aznar et al., 2002) could model the influence of variations of the fracture gap size on the pattern of healing fractures. The model was used to predict the healing process of fractures with three different gap sizes, and the predicted results were compared with those of a well-described sheep experiment (Claes et al., 1997, 1998; Augat et al., 1997). We found a close match between predicted and experimental results, in terms of callus geometry, tissue differentiation pattern and ranking of callus size and fracture stiffness. Differences between simulation and experiment in terms of absolute size and stiffness could be well explained by differences in bone sizes. The simulations used the approximate dimensions of a human tibia as a starting template. Because a human tibia is larger than that of a

sheep, the difference may cause differences in callus size and fracture stiffness. Rankings of sizes may however be affected less. This work expanded on an earlier model (Garcı´ aAznar et al., 2002). This model had a number of poorly tested assumptions. A second aim of the current study was to address these assumptions through comparison with the above experimental results. The first poorly tested assumption concerned the differentiation of fibroblasts from stem cells. Fibroblasts and stem cells are difficult to distinguish in healing fractures, so an assumption is needed. The earlier model assumed that maturation time for fibroblasts to form increased linearly with mechanical stimulus. This assumption was based on the general observation of an inverse relation between cell proliferation and cell differentiation (Malaval et al., 1999). Because the earlier and the current model both assumed that cell proliferation rates increase with mechanical stimulus, this general observation would suggest that the same increased mechanical stimulation would delay differentiation. In the current study, simulations with this positive linear rule prohibited the appearance of fibroblasts when the fracture gap was large, leaving a callus mainly consisting of granulation tissue. This was not consistent with experimental evidence. The assumption that maturation time decreased with increasing mechanical stimulus gave correct predictions of callus geometry, composition and stiffness. Such an assumption would contradict the general observation of delayed differentiation with increasing proliferation rates. One possible explanation could be that stem cell proliferation does not increase continuously with increasing stimulus, as assumed in this study, but peaks at an intermediate level. To our knowledge, the effect of mechanical strain on stem cell proliferation has not been investigated over a sufficiently wide strain regime to resolve this issue. However, a relatively quick maturation of stem cells to fibroblasts at high strain would suggest that early in the healing process populations of stem cells and fibroblasts exist together, contributing to the apparent confusion between them (Cooper and Melton, 1992; Postacchini et al., 1995). The second poorly tested assumption in the earlier model concerned bone formation. The presence of bone cells was assumed to be the only prerequisite for bone to form. Bone cells were assumed to differentiate either from stem cells, leading to intramembranous bone formation, or from cartilage cells, leading to endochondral bone formation. Differentiation to bone cells would occur wherever the mechanical stimulus level on stem cells or cartilage cells was below a threshold level. In other words, bone would form wherever the stimulus level were low. First of all, this assumption disregards the role of blood supply, which is required for bone formation to

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

occur (Winet et al., 1990). Blood supply relies on new blood vessels, which grow out from existing vessels. Such vessels are normally found in mature and newly formed bone and along the outer surface of the callus but are not observed in cartilaginous regions. Second, the assumption that bone develops wherever the mechanical stimulus is low enough for bone cells to differentiate leaves open the possibility for bone to form as isolated pockets of bone completely surrounded by granulation tissue or cartilage. To our knowledge, such isolated pockets of bone have not been observed histologically. One explanation for the absence of isolated pockets could be the role of blood supply, which tends to promote bone formation in areas adjacent to vascularized regions, i.e. bone and the callus surface. In the present study, we have made an effort to take into account the role of blood supply during bone formation. We did this by using random migration of osteoblasts from regions of bone as an indicator of vessel formation. The migrating osteoblasts play the part of migrating endothelial cells, simultaneously ensuring bone forms from existing bone only. This approach allowed us to limit the number of cell types in the model, but is likely to be over-simplistic. For one thing, it neglects vessel ingrowth from the callus surface, an area where no bone exists but which is wellvascularized. It also makes it more difficult to take into account known potent stimulators of vessel formation, such as hypoxia following the increased oxygen consumption of proliferating cells (Giordano and Johnson, 2001). When using the above approach to investigate the influence of vessel formation on the pattern of bone formation in our simulations, we found virtually no effect during intramembranous bone formation but a clear effect during endochondral bone formation. In the latter case, bone formed in cartilage away from the main ossification front when the influence of vessel formation was omitted. This effect is only possible when there is enough vascularization in the surrounding soft tissue of the callus. Fact that has not been considered in this first work. The limited influence of vessel formation on the pattern of intramembranous bone formation may be related to the small strains in granulation tissue that we assumed a prerequisite for this type of bone formation. Such small strains in a flexible tissue are most readily found adjacent to stiff regions, in other words adjacent to bone regions. Cartilage is stiffer than granulation tissue and compared to intramembranous bone formation we assumed a larger threshold strain value for cartilage to calcify. Hence, isolated bone formation is more likely in cartilage and hence the role of blood vessels more prominent. We investigated a possible role of mechanical stimuli on the rate of vessel formation by assuming a negative

117

correlation between mechanical stimulus and osteoblast migration coefficient. Such a relation follows the notion that excessive strain damages the small growing vessels such as found in healing fractures (Glowacki, 1998), thus impeding their growth. The largest effect in our simulation was seen in the fracture gap, where the influence of mechanical stimulus on vessel growth delayed bone formation. The gap area of a fracture tends to be the last to fill, which may well be due to the detrimental effect of strain on vessel growth. The final poorly tested assumption in the early model was the use of a diffusion coefficient increasing with time for simulating random migration of mesenchymal stem cells. This time-dependency was introduced as an implicit method to take into account slower migration of cells through fibrin clot and debris tissue, which is the tissue formed initially. In the current study we have introduced an explicit dependency of cell migration coefficient and the presence of disrupted tissue (or debris tissue). The assumption that MSCs migrate slower through disrupted tissue was based on the argument that the fibres in this tissue are probably spaced widely. Cell migration rates decrease with increasing fibre spacing (Kuntz and Saltzman, 1997), perhaps because in this situation not enough attachment points exists for cells. Cells rely on matrix attachments to migrate (Horwithz and Parsons, 1999). Assuming only decreased migration rates for less dense tissue may be too simplistic. In order to migrate through denser tissue, mesenchymal cells need to break down matrix because the fibre spacing is much narrower than the cell size (Lutolf et al., 2003). The need for break down makes cell migration slower through dense matrix and faster through less dense matrix (Lutolf et al., 2003). Combined, the two effects are likely to result in a biphasic relation between matrix density and migration speed. Taking into account this biphasic relation would refine the cell migration model. Like all models, this mechano-biologically regulated model of skeletal differentiation was based on several simplifications, which may limit its validity. The most important of these limitations relate to the general aspects of the model were: (1) local mechanical stimulus and time were the only factors controlling cellular events, (2) the mechanism for callus growth was limited to MSC proliferation and cartilage cell hypertrophy, (3) migration and proliferation of cartilage cells and fibroblast were assumed negligible, (4) the role of blood supply in bone formation was coupled to migration of osteoblasts from existing bone and was therefore dependent on the initial location of bone tissue. Next to these general limitations of the model, additional limitations applied to this specific study were: (5) the periosteum was the only source of MSCs, although in

ARTICLE IN PRESS 118

M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119

reality the endosteum might form an additional source, (6) the interfragmentary movement was regarded as an independent variable, but in reality the fracture healing course may affect interfragmentary movement, and (7) an axisymmetric finite element model was used, leading to an idealized cylindrical bone with axial loading only, thus neglecting finer details of geometry and load such as flexion, torsion and shear of the bone. Despite these simplifications, general agreement between simulation and experiment was found in terms of callus geometry, callus size, fracture stiffness and tissue differentiation pattern. In conclusion, the model was able to simulate the increasing callus size and delay in healing when gap size increase from 1 to 2 mm, and the very small callus volume and non-union when the gap is increased further to 6 mm, as observed by Claes et al. (1998). At the same time the model was able to predict the temporal evolution of the callus stiffness for the three different gaps. This correspondence between results and experiments suggests that the observed effects can be explained largely from the mechano-biological model presented here.





aproliferation ¼ 0:85 day1 and cproliferation ¼ 0:2 are the parameters that define stem cell proliferation. csmax ¼ 100 000 cells=mm3 is the saturation value of stem cell density (Wilsman et al., 1996). khyper ¼ 1000 chondrocytes=ðmm3 dayÞ is the parameter that controls the volume increase produced during cartilage calcification (Wilsman et al., 1996). Qb ¼ 3  106 mm3 ðosteoblast dayÞ; Qs ¼ 3  106 6 3 3 mm =ðstem cell dayÞ; Qc ¼5  10 mm =ðchondrocyte dayÞ; Qf ¼ 5  106 mm3 =ðfibroblast dayÞ is the amount of matrix per bone cell, MSC, cartilage cell or fibroblast and unit time (Martin et al., 1998). ccmin ¼ 60 000 chondrocytes=mm3 (Wilsman et al., 1996) is the equilibrium value of cartilage cell concentration. pmin mi ¼ 0:015 is the minimum mineralization level needed for endochondral ossification to take place. D0 ¼ 1 mm2 =day is the diffusion coefficient that controls stem cell migration. Db ¼ 4 mm2 =day constant that defines the diffusion coefficient in the ossification process. cmin ¼ 20 000 cells=mm2 is the minimum concentrab tion of bone cells that indicates when blood supply is completed at the ossification front and osteoblasts can differentiate directly from stem cells.

Appendix A The parameters used in the simulation were:





cbone is the stimulus level below which bone tissue may be formed; this value has been estimated as 0.03. clamellar ¼ 0:008; is the stimulus level below which woven bone begins to remodel into lamellar bone. clim ¼ 0:006 is the minimum stimulus level required for mesenchymal stem cells to differentiate into bone cells during intramembranous ossification (Turner et al., 1994). ccartilage ¼ 0:1 is the minimum stimulus value necessary for stem cells to differentiate into cartilage cells. ccalcified ¼ 0:06 is the stimulus level below which cartilage begins to calcify before endochondral ossification occurs. cfibrous ¼ 0:8 is the stimulus level below which stem cells may differentiate into fibroblasts. cslope ¼ 0:28 is the stimulus level which defines the change in slope of the bilinear differentiation rule of fibroblasts. cdeath ¼ 1 is the stimulus level above which MSCs will die (Claes et al., 1997). cdisruption ¼ 0:6 is the stimulus level above which disruption of granulation and debris tissue matrix takes place. tcmlim ¼ 15 days is the timepoint at which cartilage cells start to appear in the callus (Cullinane et al., 2003), tslope mlim ¼ 200 days is the time limit that define the curve of fibroblast differentiation rule.

References Ament, C., Hofer, E., 2000. A fuzzy logic model of fracture healing. J. Biomech. 33, 961–968. Aubin, J.E., Liu, F., Malaval, L., Gupta, A., 1995. Osteoblast and chondroblast differentiation. Bone 17, 77S–83S. Augat, P., Merk, J., Genant, H.K., Claes, L., 1997. Quantitative assessment of the experimental fracture repair by peripheral computed tomography. Calcified Tissue Int. 60, 194–199. Augat, P., Margevicious, K., Simon, J., Wolf, S., Suger, G., Claes, L., 1998. Local tissue properties in bone healing: influence of size and stability of the osteotomy gap. J. Orthop. Res. 16, 475–481. Bailo´n-Plaza, A., van der Meulen, M.C., 2003. Beneficial effects of moderate, early loading and adverse effects of delayed or excessive loading on bone healing. J. Biomech. 36, 1069–1077. Barber, C.B., Dobkin, D.P., Huhdanpaa, H., 1996. The quickhull algorithm for convex hulls. ACM Trans. Math. Software 22 (4), 469–483. Beaupre´, G.S., Orr, T.E., Carter, D.R., 1990a. An approach for timedependent bone modeling and remodeling-theoretical development. J. Orthop. Res. 8, 551–651. Beaupre´, G.S., Orr, T.E., Carter, D.R., 1990b. An approach for timedependent bone modeling and remodeling-application: a preliminary remodeling simulation. J. Orthop. Res. 8, 662–670. Bland, Y.S., Critchlow, M.A., Ashhurst, D.E., 1999. The expression of the fibrillar collagen genes during fracture healing: heterogeneity of the matrices and differentiation of the osteoprogenitor cells. Histochem. J. 31, 797–809. Breur, G.J., VanEnkevort, Farnum, C.E., Wilsman, N.J., 1991. Linear relationship between the volume of hypertrophic chondrocytes and the rate of longitudinal bone growth in growth plates. J. Orthop. Res. 9, 348–359.

ARTICLE IN PRESS M.J. Go´mez-Benito et al. / Journal of Theoretical Biology 235 (2005) 105–119 Carter, D.R., Beaupre´, G.S., Giori, N.J., Helms, J.A., 1998. Mechanobiology of skeletal regeneration. Clin. Orthop. 355S, 41–55. Claes, L., Heigele, C.A., 1999. Magnitudes of local stress and strain along bony surfaces predict the course and type of fracture healing. J. Biomech. 32 (3), 255–266. Claes, L., Wilke, H.J., Augat, P., Suger, G., Fleischman, W., 1994. The influence of fracture gap size and stability on bone healing. Trans. Orthop. Res. 19, 203. Claes, L., Wilke, H.J., Kiefer, H., 1995. Osteonal structure better predicts tensile strength of healing bone than volume fraction. J. Biomech. 28 (11), 1377–1390. Claes, L., Augat, P., Suger, G., Wilke, H.J., 1997. Influence of size and ans stability of the osteotomy gap on the success of fracture healing. J. Orthop. Res. 15 (4), 577–584. Claes, L., Heigele, C.A., Neidlinger-Wilke, C., Kaspar, D., Seidl, W., Margevicius, K., Augat, P., 1998. Effects of mechanical factors on the fracture healing process. Clin. Orthop. 355S, 132–147. Connolly, J.F., 2001. Nonoperative fracture treatment. In: Bucholz, R.W., Heckman, J.D. (Eds.), Rockwood and Green’s Fractures in Adults, fifth ed. Lippincot, PA, pp. 113–180. Cooper, C., Melton, L.J., 1992. Hip fractures in the elderly: a worldwide projection. Osteoporosis Int. 2, 285–289. Cullinane, D.M., Salisbury, K.T., Alkhiary, Y., Eisenberg, S., 2003. Effects of the local mechanical environment on vertebrate tissue differentiation during repair: does repair recapitulate development. J. Exp. Biol. 206, 2459–2471. Einhorn, T.A., 1995. Current concepts review: enhancement of fracture-healing. J. Bone Joint Surg. 77-A, 940–955. Garcı´ a-Aznar, J.M., Kuiper, J.H., Go´mez-Benito, M.J., Doblare´, M., Richardson, J.B. Computational simulation of fracture healing. In: Proceedings of the 13th Conference of European Society of Biomechanics, Wroclaw, Poland, pp. 394–395. Gardner, T.N., Mishra, S., 2000. The mechanical environment of a bone fracture and its influence upon the morphology of healing. Med. Eng. Phys. 33, 415–425. Gardner, T.N., Stoll, T., Marks, L., Mishra, S., Tate, M.K., 2000. The influence of mechanical stimulus on the pattern of tissue differentiation in a long bone fracture—an FEM study. J. Biomech. 33, 415–425. Giordano, F.J., Johnson, R.S., 2001. Angiogenesis: the role of the microenvironment in flipping the switch. Curr. Opin. Genet. Dev. 11, 35–40. Glowacki, J., 1998. Angiogenesis in fracture repair. Clin. Orthop. Related Res. 355, S82–S89. Goodship, A.E., Watkins, P.E., Rigby, H.S., Kenwright, J., 1993. The role of fixator frame stiffness in the control of fracture healing. An experimental study. J. Biomech. 26 (9), 1027–1035. Hadjiargyrou, M., McLeod, K., Ryaby, J.P., Rubin, C., 1998. Enhancement of fracture healing by low intensity ultrasound. Clin. Orthop. 355S, 216–229. Hibbit, Karlsson and Sorensen, Inc., 1999. Theory Manual, v. 5.8. HKS Inc. Pawtucket, RI, USA. Horwithz, A.R., Parsons, J.T., 1999. Cell migration—movin’ on. Science 286, 1102–1103. Ito, K., Wilson, W., van Donkelaar, C.C., Huiskes, R., 2004. Comparison of mechano-regulation parameters for tissue differentiation during fracture healing. In: 50th Annual Meeting of the Orthopaedic Research Society, Poster No. 480. Iwaki, A., Jingushi, S., Oda, Y., Izumi, T., Shida, J.I., Tsuneyoshi, M., Sugioka, Y., 1997. Localization and quantification of proliferating cells during rat fracture repair: detection of proliferating cell nuclear antigen by immunochemistry. J. Bone Miner. Res. 12, 96–102.

119

Jacobs, C.R., 1994. Numerical simulation of bone adaptation to mechanical loading. Ph.D. Thesis, Stanford University. Kelly, D.J., Prendergast, P.J., 2004. A mechanobiological analysis of osteochondral defect repair. In: 50th Annual Meeting of the Orthopaedic Research Society. Poster No. 666. Kenwright, J., Gardner, T.N., 1998. Mechanical influences on tibial fracture healing. Clin. Orthop. 355, 179–190. Kuiper, J.H., Richarson, J.B., Asthton, B.A., 2000. Computer simulation to study the effect of fracture site movement on tissue formation and fracture stiffness restoration. In: European Congress on Computational Methods in Applied Sciences and Engineering, Barcelona, Spain. Kuntz, R.M., Saltzman, W.M., 1997. Neutrophil motility in extracellular matrix gels: mesh size and adhesion affect speed of migration. Biophys. J. 72 (3), 1472–1480. Lacroix, D., 2001. Simulation of tissue differentiation during fracture healing. Ph.D. Thesis. Lacroix, D., Prendergast, P.J., 2002. A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading. J. Biomech. 35, 1163–1171. Levick, J.R., 1987. Flow through interstitium and other fibrous matrices. Q. J. Exp. Physiol. 72, 409–438. Lutolf, M.P., Lauer-Fields, J.L., Schmoekel, H.G., Metters, A.T., Weber, F.E., Fields, G.B., Hubbell, J.A., 2003. Synthetic matrix metalloproteinase-sensitive hydrogels for the conduction of tissue regeneration: engineering cell-invasion characteristics. Proc. Natl Acad. Sci. 100 (9), 5413–5418. Malaval, L., Liu, F., Roche, P., Aubin, J.E., 1999. Kinetics of osteoprogenitor proliferation and osteoblast differentiation in vitro. J. Cell Biochem. 74 (4), 616–627. Martin, R.B., 1989. Porosity and specific surface of bone. In: CRC Critical Review of Biomedical Engineering, pp. 179–222. Martin, R.B., Burr, D.B., Sharkey, N.A., 1998. Skeletal Tissue Mechanics. Springer, Berlin. Olsen, L., Sherratt, J.A., Maini, P.K., 1995. A mechanochemical model for adult dermal wound contraction: on the permanence of the contracted tissue displacement profile. J. Theor. Biol. 177, 113–128. Postacchini, F., Gumina, S., Perugia, D., De Martino, C., 1995. Early fracture callus in the diaphysis of human long bones, histologic and ultrastructural study. Clin. Orthop. 310, 218–228. Pratap, J., Galindo, M., Zaidi, S., Vradii, D., Bhat, B.M., Robinson, J.A., Y., C.J., Komori, T., Stein, J.L., Lian, J.B., Stein, G.S., van Wijnen, A.J., 2003. Cell growth regulatory role of runx2 during proliferative expansion of preosteoblasts. Cancer Res. 63, 5357–5362. Prendergast, P.J., Huiskes, R., Soballe, K., 1997. Biophysical stimuli on cells during tissue differentiation at implant interfaces. J. Biomech. 6, 539–548. Sarmiento, A., Latta, L.L., 1995. Functional Fracture Bracing: Tibia, Humerus and Ulna. Springer, Berlin. Smith, R.K., Roberts, W., 1980. Cell kinetics of the initial response to orthodontically induced osteogenesis in rat molar periodontal ligament. Calcified Tissue Int. 30, 51–56. Turner, C.H., Forwood, M., Rho, J.Y., Yoshikawa, T., 1994. Mechanical loading thresholds for lamellar and woven bone formation. J. Bone Miner. Res. 9, 87–97. Wilsman, N.J., Farnum, C.E., Leiferman, E.M., Fry, M., Barreto, C., 1996. Differential growth by growth plates as a function of multiple parameters of chondrocytic kinetics. J. Orthop. Res. 14, 927–936. Winet, H., Bao, J.Y., Moffat, R., 1990. A control chamber for tibial cortex neovascularization in the bone chamber. J. Bone Miner. Res. 5, 19–30.