Simultaneous Presence of Growth and Remodeling in the Bone Adaptation Theory

Mechanical forces acting on bone during growth affect their final shape and strength. Mechanoregulation of bone growth is maybe recognized in embryogenesis, and also in the adaptation of the adult skeleton to changes in mechanical loading. By combining equations describing bone remodeling and growth with an iterative finite element analysis, a computational model to simulate the simultaneous effects of bone remodeling and bone growth was proposed in this study. Strain-energy density was assumed as mechanical stimulus of bone adaptation process. Negative exponential decay function over time was considered as metabolic growth rate. Based upon numeric results, the model shows an acceptable behavior under various modes of loading, e.g. altering in trabecula’s orientation or its thickness. This model also shows that by neglecting growth part in the adaptation model, a considerable error would result in both final density distribution and microstructural pattern of spongy bone.


INTRODUCTION
The development, growth, and remodeling of skeletal structures is a highly regulated process beginning with mesenchymal stem cells condensations in the early embryo and finishing with the homeostatic skeleton of the adult. It is widely accepted that both genetic and epigenetic factors determine the final shape and strength of the skeleton, and many authors have specifically proposed an epigenetic role for mechanical forces [1,2,3] . Equations have been proposed to describe how mechanical forces modulate growth where a mechanobiological remodeling rate is superimposed on a baseline biological growth rate [4] .
In 1976 a rigorous mathematical bone remodeling model was proposed by Cowin and Hegedus which is so-called adaptive elasticity theory. In the adaptive elasticity theory, it is assumed that the rate of change in bone mass is correlated with the history of mechanical strain [5] . Skalak et al. later formulated a continuum model of growth [6] . Carter [7] proposed a bone remodeling model which was targeted to produce a homeostatic level of an effective stimulus and at the same time, Huiskes et al. proposed that the process of bone remodeling is aimed at producing a homeostatic value of the strain energy density in the tissue [8] .
Bone growth models incorporating both biological and mechanobiological influences have been proposed by Van der Meulen and colleagues [9] for modeling the cross-sectional growth of long bones and by Stevens et al. for modeling endochondral growth using a finite element model [4] . An alternative approach is that remodeling of bone is in response to microdamage, either to regulate microdamage to a homeostatic level or as the stimulus for the activation of the coupled responses of osteoclast and osteoblast cells [10,11,12] .
Recently, Vahdati and co-workers incorporated a cellular accommodation effect into the Huiskes et al. 's semi-mechanistic model [13] . They showed that the model is sensitive to temporal sequence in which loading is applied. In a very recent effort, Vahdati and Rouhi proposed another modification on the Huiskes et al's model [13] . They included the effects of both microcracks and disuse on activation of resorption in one unifying formulation based on latest experimental findings besides considering cellular accommodation effect. These findings have not been published, yet. In 2004, Rouhi and co-workers proposed a modification on the adaptive elasticity theory by replacing volume fraction with free surface density in the constitutive equations [14] . In another attempt, Rouhi et al. incorporated a microcrack factor in their first model and showed that not only mechanical stimuli, but also their rate and history are effective and at play in the bone remodeling process [15] . Considering the great importance of bone resorption process in the osteoporotic cases, Rouhi et al. proposed a separate model of bone resoprtion by using mixture theory with chemical reactions. In their bi-phasic mixture model, they found that not only mechanical factors are at play in the resorption process, but also chemical and biological factors have crucial effects [16] . Considering three different constituents of bone, i.e. bone matrix, bone fluid, and bone resorbing cells; Rouhi proposed a tri-phasic model of bone resorption using mixture theory with chemical reactions [17] . It is concluded that rate of bone resorption is a function of apparent density of bone matrix and bone fluid, fluid velocity, momentum supply to the fluid phase, and internal energy densities of different constituents, in the former model.
Frost demonstrated that not only mechanoregulation diagram is not plateau in the lazy zone for growing bone, but also it has a considerable time dependent positive slope. This positive slope represents modeling process ( Fig. 1) [18] .

MATERIALS AND METHODS
It is widely accepted that immature bone structure can be altered via two major mechanisms: (1) Metabolic base growth, which is dependent on genetic patterns, nutrition, etc., (2) Mechanoregulatory mechanisms. The former is called Genetic Factor while the latter is called Epigenetic Factor [19] . Generally speaking, it can be applicable superimposing Genetic to Epigenetic effects [20] .
Computational simulations are used to investigate the mechanoregulation of bone growth. The bone growth model employed here was proposed by Van der Meulen et al. [9] . It simulates the growth of the crosssection of a long bone from an embryonic bone collar to maturity, where the rate of bone apposition, t ρ , is equal to the sum of the baseline biological rate, b ρ , and the rate due to mechanobiological effects, m ρ , as defined in Fig. 1: Mechanoregulation diagram differs in mature and immature bones. Positive slope in the lazy zone represents modeling process in growing bone [18] .
The baseline biological growth rate is a decaying exponential function of time that reduces to approximately zero by the age of 18 years, defined as follows [9,20] : where = 3.6 years, b0 • ρ is the initial biological growth rate in gr.cc -1 .day -1 , arisen from Eq. (3) in which initial r is taken to be 10 g.day -1 in this model [9] , where r is the rate of change in trabecular thickness, S v is the specific surface area and t is the density of bone matrix which is assumed to be constant and equal to 2 g.cc -1 . Figure 2 shows how r changes as an individual ages [9] . Figure 3 shown the specific surface area (S v ), as a function of porosity (p) with coefficients in mm 2 .mm -3 as follows [3]  . SSA is near zero at lower and higher extremes of p [21] The finite element stress analysis is used to calculate stress and strain within bone. The strain energy density, SED, can be expressed as follows: where σ ij and ij are the stress and strain tensors, respectively and e represents the number of each element in the model. The resulted SED is assumed to act as the mechanical stimulus of the remodeling process. The bone remodeling set of equations was introduced by Ruimerman et al. in the following form [22] : As stated earlier, since there are two mechanisms contributing in the bone adaptation process before maturation, i.e. Genetic and Epigenetic; we hypothesized here to superimpose their effects. Thus, adaptation equation will consist both growth rate, b ρ , and the remodeling rate, m ρ . It is performed by inserting a linear alternative, instead of plateau, into the lazy zone of common remodeling models, appeared in Eq.
where u R is the lower or resorption SED threshold, u D is the upper or deposition SED threshold, ρ e represents elemental density and c is a growth rate constant. The values of u R and u D are age-dependent and will approach the mature values with linear approximate variation within 18 years. The value of u R starts from zero at birth and increases until reaches its mature value, i.e. 0.002 J mm 3 . Based on Ruberg's study [23] , u R is considered as 80 percent of midpoint SED in lazy zone, while u D is assumed to be 20 percent greater than midpoint SED in the lazy zone. By this assumption, u D is also a linear function of time during first 18 years of life. The superposition process also alters the resorption and deposition zones. This is carried out via replacing common models' constant, i.e. c, with its alternative for growing bone [24,25] . As it is quite evident in Fig. 1, the proportionality constant c in immature bone is not the same as in mature bone. Thus, we have considered a linear approximate temporal variation, like the ones for lazy-zone thresholds, which starts from its initial value, i.e. c = 0.01 g.day.mm -3 ., and approaches its mature value, i.e. 0.02 g.day.mm -3 , at the age of 18 years [23] . In other words, the constant c is, in fact, a time dependent coefficient which affects the rate of bone adaptation under the effect of external load placed on the bone. A higher value of c will lead to more expansion and a thicker bone; it also acts as a handle for controlling the convergence [20] .
Concentrating on Eq. 7, it is obvious that the second term imposes a positive, age-decaying slope to mechanoregulation diagram over the lazy zone. Assuming that the maturity happens at the age of 18 years, this model converges to common remodeling models proposed by others (e.g. Jacobs [28] ,Carter et al. [19] , Huiskes et al. [12] ), for mature bone.
As a physical run-time constraint, the elemental densities were assumed not to outrage the following range (0.001, 1.74 g/cc). The lower bound is taken 0.001 instead of zero due to some computational considerations. Zero density will lead to zero stiffness, internal stress shielding, checkerboard configuration, and even singularity in global stiffness matrix of structure in finite element solution [26] . Thus, the following constraint was set: e g g 0.001 1.74 cc cc < ρ < The Van der Muelen et al.'s model uses agedependent loading cases to simulate the loading history for an individual bone (Fig. 5). This force was obtained from body weight during growth, and it is assumed that the force and age are linearly proportional. The load starts from approximately zero values and gradually increases with a mean slope of 1.11 Newton per year, until reaches 20 Newtons at the age of 18 years [27] .
The relation between the elastic modulus and the bone density is assumed to be in the following form as proposed by Jacobs in 1994 [ where e ρ refers to the density of each element. In addition to elastic modulus regulation, the model proposed here uses the following equation, also proposed by Jacobs [28] where, F is target function and e is the summation index number over elements. First and second degrees of freedom of first node at left-bottom corner and first degree of freedom of 25 th node at right-bottom corner were restrained as minimum required boundary conditions. We also prevented overconstrained condition in applying boundary conditions. Applied boundary conditions are shown in Fig. 7. These two factors also converge to zero at 25 nodes per side and higher.  [3,29] Post-processing was performed to calculate stresses and strains at gauss integration points. Then, the state of remodeling, i.e.
where e (t t) ρ + δ is the updated density of n th element which will lead to updated young's moduli and Poisson's ratios. The iterative solution procedure includes two loops which one fulfils temporal simulation and the other stands for convergence testing.
The convergence loop is the inner and the temporal loop is the outer one. Number of temporal iterations was considered 75 which approximately reveals 20 th year of age, assuming 365 days each year, and 100 days each iteration.  (13) where, V represents geometry volume and v b shows bone volume fraction. Stoppage threshold of e was set between 10 -4 and 10 -5 . Considering , in which V B is volume of bone tissue and V T represents total bone specimen volume, this formulation leads to Eq. (14). Algebraic procedure can be found in [23] . where e represents an arbitrary element within the geometry. Fig.8 shows a general flowchart of our solution.

RESULTS AND DISCUSSIONS
Four different simulations were performed to determine whether our model shows reasonable behavior under various loading conditions. The conceptual initial architecture was considered as a uniform density distribution with a density of 0.7 g/cc [29] . In the figures 9-12, final configurations of  In the first simulation, we have applied a 1 Newton force on the upper side of square which was oriented 45 o clockwise with respect to horizon (Fig. 9a). Then, we changed the orientation of force to 45 o counterclockwise, with same magnitude. Figures (9b-e) show the adaptation process due to change in external force orientation. It is obvious that hard tissue mass is gradually transferred from initial major strut to second major strut. It is strongly in agreement with Huiskes et al. findings in [12] . Both final architectures show similar configurations and mass distributions (Fig. 9). In the second test, we applied a uniformly distributed shear force on the right side of the square. As can be shown in Fig. 10, model has regulated itself by producing two struts in response to loading; first just beneath loaded nodes and second is developed with 45 o angle oblique to the first. The first strut is developed in the direction of maximum shearing stress and the second is produced in the direction of principal direction of normal stresses. Based on Wolff's law, some researchers have interpreted that bone reinforces itself in orientation which is maximally stressed [7,8,28] .
Our 2 nd simulation is in agreement with the former interpretation of Wolff's law.
In the third simulation, the bending condition was simulated by applying a bi-ramp on the upper side of the model. As can be shown in Fig. 11, the final configuration of this case has a strong qualitative correlation with distal femur macroscopic architecture. Asymmetricity found in final configuration of this loading case is interpreted in relation with asymmetric boundary conditions. Trabecular thickness, which is one of the most important geometric parameters relevant to osteoporosis and also osteopenia [30] , is the next issue we are going to discuss. In the forth and last qualifying round of tests, initially we have loaded the model with an oblique 1 Newton force. Then, we have over-and under-loaded the model by 25% of the original force (1N). Graphical presentation of such test can be seen in Fig. 12. As it is obvious in the Fig. 12, over-loading condition induces an elevation and under-loading or disuse, induces a reduction in trabeculae thicknesses. It is widely expressed by authors that the former (overloaded one) causes a hypertrophy, and the latter (the under-loaded and/or disused case) will lead to an atrophy and osteoporotic condition in a healthy bone [31,32] .
In this model, nothing was said about growth yet. We have introduced our growth included model (GIM) by Eq. 7. We also mentioned growth excluded model (GEM) in Eq. 6. After applying these two models on loading cases, specimens' masses were drawn as a function of age. All the diagrams show qualitative resemblance, but not numeric. We have included the diagrams of first test results for both GIM and GEM in Fig. 13, for instance. As the figure shows, both models start from initial mass in the order of 0.61 grams. As individual ages, the GIM accelerates more than GEM due to the growth term. This situation continues until the growth is stopped at the age of 18. The difference between GIM and GEM at the age of 18 is about 80mg. This difference represents approximately 9.3% of weight for homeostatic mature GIM in first loading case, which is 0.856 grams. The aforementioned difference can be considered as the error included in the model due to growth effect cancellation assumption.  Figure 14 contains a diagram which is calculated by subtracting mass of the growth excluded model (GEM) from the mass of growth included model (GIM) which is called as "impact of growth assumption" by the authors. Based on this figure, convergence takes place at the age of 18. This shows that bone mass calculated using the GIM will approach the one of derived employing the GEM after maturation. In addition to that, this diagram resembles a step response of a 2 nd -order closed-loop control system which can be studied considering the magnitude, nature and physiologic interpretation of its particular System parameters, i.e. damping ratio and natural frequency. Some authors have interpreted system parameters of such system as a gauge of bone health and youth and its ability to response to remodeling stimulus [33] . The rate of change in bone mass diagram is also studied in this research (Fig. 15). Figure 15 shows the rate of change of bone mass as individual ages. The values of this parameter are 0.0436g/day and 0.0674g/day at birth for GEM and GIM, respectively. This difference is relevant to the initial rate of growth. Figure 16, also shows temporal variation in the difference between GEM and GIM. The parameter decays exponentially with time and the numeric analysis shows that the time constant is equal to that of the growth rate equation, introduced in Eq. 2.

CONCLUSIONS
In this research, the authors are trying to put an emphasis on the point that growth and remodeling can be considered simultaneously to predict homeostatic structures of growing bone during the first two decades of life. It is well accepted that in a growing bone, or immature bone, both genetic and epigenetic factors are at play. In other words, both modeling and remodeling are active. Thus, in this study, the authors superimposed bone modeling equation on bone remodeling equation, in the hope of modeling immature bone adaptation.
From the simulations performed in this study, it is concluded that this model can produce strut-and platelike structures found in substructures of spongy bone; as was resulted in other studies [1,23] . All the homeostatic configurations show a fair agreement with Wolff's law [3] ; topologic optimization methods [34] ; bone microstructure; and also other researchers' works [22,35] . This study also shows that neglecting growth effect in immature bone models can cause a noticeable error into the final results. In this study, neglecting growth effect in bone adaptation process could cause an error of about 9.3% in final bone mass.
It is evident that the biomechanics society is suffering from the lack of a mechanistic model of bone adaptation. One reason of this shortage comes from the fact that bone adaptation is a multidisciplinary problem which is related to different fields (e.g. cell biology; mechanics; and chemistry). One of the most important challenges in the field of bone adaptation is mechanotransduction of bone cells. In other words, understanding how mechanosensors, i.e. osteocytes, react to different mode of loadings. It is hoped that by closer collaboration among researchers from medicine and engineering more insight will be gained in the near future, and more lights will be shed on this very essential and complex problem of biomechanics.