Farid Amirouche; and Aimee Bobko

Review Article

Austin J Biotechnol Bioeng. 2015;2(2): 1041.

Bone Remodeling and Biomechanical Processes- A Multiphysics Approach

Farid Amirouche¹* and Aimee Bobko²

¹Department of Mechanical Engineering, University of Illinois, USA

²College of Medicine, University of Illinois, USA

*Corresponding author: Farid Amirouche, Department of Mechanical Engineering, University of Illinois at Chicago, Chicago, IL, USA.

Received: March 04, 2015; Accepted: May 12, 2015; Published: May 14, 2015

Abstract

Bone remodeling is bone’s mechanism to optimize its structure in response to the body’s external and internal stimuli. From the mineralized matrix of whole bone to the lacuna canalicular network of osteocytes, mechanical and biological signaling induce changes in properties like bone density. Previous studies have made strides in simulating aspects of these processes with mathematical models. The focus is to now create a more encompassing assessment by incorporating the interdependence between bone’s macrostructure and microstructure. This work is emphasizing a multiphysics and multiscale approach in creating a model to predict the age-related changes in bone remodeling.

Introduction

From a macroscopic perspective, bone is a biological entity which is designed to obtain optimal mechanical properties to enhance the structure’s ability to support and locomote the body. Bone remodeling, in addition to providing everyday maintenance to the structure, is the phenomenon responsible for the bone’s adaptation to stimuli. In order to achieve these modifications, there is a highly active and organized system present at the cellular, microscopic level. Bone remodeling is a location specific cyclical process generally described as a resorption phase completed by cells called osteoclasts, followed by a formation phase consisting of laying down new bone matrix performed by cells called osteoblasts. Locations are designated as a site in need of bone remodeling by cells called osteocytes which become activated under increased mechanical stress. This stress can appear randomly due to the formation of microcracks or can be in direct response to a change in external mechanical loads [1]. Without the optimizing ability of the bone remodeling process, bone would be at greater risk for points of structural weakness and the potential to fracture. This is a phenomenon that due to its overlapping processes is still not well understood and has been the subject of multilevel of research.

Bone Remodeling Review

After bone is originally formed or “modeled,” the process of bone remodeling begins to take place in order to ensure the upkeep of the bone’s structural integrity. A general description of bone remodeling is the resorption of old bone, followed by the formation of new bone at designated locations, resulting in ‘targeted remodeling’ [2]. More in detail, Figure 1 shows a cyclical representation of bone remodeling divided into 5 stages: activation, resorption, reversal, formation, and quiescence [3]. Reversal is referring to the time, typically 1-2 weeks, between the resorption and formation processes at one location [4]. This cycle emphasizes the sequence of events that occur to maintain the stable architecture of bone.

Figure 1: Bone remodeling process depicted as a five stage cycle [3].

    
    
    Figure 1:  Bone remodeling process depicted as a five stage cycle [3].

Bone multicellular units

Locations undergoing bone remodeling, also known as bone remodeling units or Bone Multicellular Units (BMUs), are numerous within one piece of bone [5]. At any time, it has been estimated that the entire body has 105 to 106 BMUs present [4]. The BMU travels approximately 20-40 μm/day and has an average lifespan of 6-12 months [6,7,8]. Osteoclasts lead the cavity at the front, resorbing bone in order to keep progressing forward [9]. In order to do so, these cells create what are known as resorption cavities. As osteoclasts keep enlarging the BMU, a growing capillary follows to not only ensure nutrients are delivered to the surrounding cells but also to provide a constant influx of osteoblast and osteoclast precursor cells from the bone marrow [9,10]. The migration of angiogenesis is a key component of bone homeostasis [11,12,13]. Under direction from a variety of signaling mechanisms, the precursor cells will undergo a specific sequence of events in order to differentiate into mature osteoblasts and osteoclasts [14]. Through a signaling mechanism, these cells position themselves into a gradient with osteoclasts proceeding to the front and osteoblasts settling to the back [9]. The osteoblasts will secrete osteoid to reform whatever bone has been resorbed. Additionally, the osteoblasts will begin the differentiation process into osteocytes. While within the unmineralized osteoid, the cells may be referred to as pre-osteocytes [15]. Once the cells are completely surrounded by bone matrix and are no longer motile, they can then be referred to as osteocytes [16-19]. The differentiation process from osteoblast to osteocyte takes approximately 3 days [20]. The osteocytes are cells with slender extensions known as the osteocytic processes. Some of these processes are in contact with the processes of another osteocyte or an osteoblast. These connections take place through the existence of gap junctions. A gap junction is a cell structure made of proteins which allows the exchange of signals through the transport of ions and second messengers between adjacent cells as well as the delivery of small metabolites [21]. This creates a full-length communicating cell network within a bone.

Figure 2 provides a systemic diagram of the processes of bone remodeling phenomenon. There are inputs and outputs of the system including mechanical load, nutrients, hormones, precursor cell supply, and waste removal. These items can help to stimulate or inhibit the remodeling process. The internal components carrying out this process include the various cell types and the secreted product of osteoblasts, osteoid, which will eventually mineralize. All of these activities are highly balanced in relation to one another. If a mechanical or biological element is changed, the entire system depicted in Figure 2 must respond appropriately to ensure that the overall appearance of homeostasis is maintained. Specifically, if the balance of the system is disturbed, pathological conditions such as osteoporosis, Paget’s disease, and cancer can result [22].

Figure 2: Bone remodeling system with labeled components and processes that are continuously adapting to changing internal and external stimuli.

    
    
    Figure 2:  Bone remodeling system with labeled components and processes
that are continuously adapting to changing internal and external stimuli.

Interscale synergies

An increased mechanical load is seen as a stimulant to the bone remodeling process as it initiates the need for greater strength at the points of highest stress. In fact, it has been found that both the osteons of cortical bone and the trabeculae of cancellous bone are aligned in the same direction as the largest mechanical loads being applied to a bone [23,24]. As a result, the bone remodeling process must take action to lay down new bone in these areas. When a load is being placed upon the macroscale mineralized bone matrix, it is also resonating into strains at the nanoscale of the lacuna canalicular network of osteocytes. Cells function to maintain or increase the strength and stiffness properties back at the macroscale level. These top-to-bottom and bottom-to-top interactions otherwise known as interscale synergies, as depicted in Figure 3, have proven to be the difficulty of creating a comprehensive mathematical model of the bone remodeling phenomenon.

Figure 3: The interscale synergies occurring within a cortical bone’s structure. The thermodynamic processes occurring at each scale are interconnected to deliver the resulting whole bone properties.

    
    
    Figure 3:  The interscale synergies occurring within a cortical bone’s structure.
The thermodynamic processes occurring at each scale are interconnected to
deliver the resulting whole bone properties.

More specifically at the microscale level, when an external mechanical load is applied to a piece of bone, stress throughout the whole tissue increases proportionally, creating pressure gradients within the interstitial fluid of the nanometer scale canaliculi which surround the processes of osteocytes [25]. As a result of these pressure gradients, interstitial fluid flow occurs, creating a shear stress on the osteocytic processes. This shear stress creates a deformation in the osteocytic process membrane [26,27]. It has been shown that the necessary threshold for deformation to elicit a response within the osteocyte cell body is 4.3± 0.8μm and in the osteocytic cell process is 2.0±0.5μm [26]. The difference in these values indicate that it takes a smaller deformation in the osteocytic cell process to elicit a response, thus the osteocytic cell process is more mechanosensitive than the osteocyte cell body [26]. In addition to the mechanical stimulus, the initiated interstitial fluid flow will carry nutrients and hormones which can act on the osteocytes and their signaling response [25]. It is hypothesized that the osteocyte takes in all of the various stimuli, combines their inputs, and produces an appropriate signal in response [28-32]. This signal is estimated to have an effect within a 100μm diameter [33]. The major signaling metabolic factors secreted by the osteocyte include those such as Nitric Oxide (NO) and prostaglandin E2 (PGE2) [34,35]. It is important to note that NO has an inhibitory effect on osteoclast cells, while PGE2 has a stimulatory effect on osteoblasts [4,36]. Hence, when a mechanical load is applied to bone, the cellular reaction is to enhance bone formation and limit bone resorption. However, in the case of disuse (unloading) where there is little interstitial flow to deliver NO from osteocytes to osteoclasts or in microdamage where canaliculi are damaged along with associated osteocytic cell processes, the inhibitory effect of osteocytes on osteoclasts is diminished [1,33,37,38]. Additionally, when osteocytes undergo apoptosis of cell death, which is stimulated by the lack of mechanical load such as with immobilization or weightlessness, there is less inhibitory signaling directed to surrounding osteoclasts, resulting in a weakened structure [39-44]. NO has a half-life of a few minutes which emphasizes that its paracrine actions are only effective in local regions under the current mechanical force or lack thereof [45]. Of note, it has been shown that depriving interstitial fluid of nutrients decreases the osteocyte’s responsiveness to secrete NO, even in the presence of a constant shear stress [46]. This demonstrates the importance of the biological cell requirements to a cell’s function and the need to account for this as a factor responsible for initiation and inhibition of bone remodeling [47].

The mechanosensitive osteocyte

Despite the limitation of knowing whether the signal for bone remodeling is more mechanically or biologically based, the presence or absence of a signal from an osteocyte is of great importance. In terms of bone maintenance, as stated previously, if micro damage has occurred, preventing an osteocyte’s signal from reaching osteoclasts, resorption of the surrounding damaged mineralized matrix will begin to occur [48]. The function of osteoclasts contributes to the stimulation of the osteocyte signal by first migrating to the site of localization with an average rate of 105±10μm/hour, then resorbing bone at an average rate of 165μm3/hour (as demonstrated in vitro) for approximately 2 weeks [34,49,50]. During this time, the osteoclasts create what is known a resorption cavity [1,33]. This newly developed architectural structure changes the stress distribution at this particular site and is also recognized as a mechanical stimulus by the osteocytes [1,33]. As a result of the stresses and mechanical load changes inflicted by the newly formed resorption cavity, the mechanosensitive osteocyte will secrete appropriate signals to recruit osteoblasts to the site so that new bone may be laid down within the newly formed cavity [1,33]. Osteoblasts will only be recruited if the osteocyte’s recruitment signals are above a threshold [1]. These cells will then secrete osteoid, which eventually mineralizes to bone [34]. The formation process can take place over a 3-6 month time period [34].

Intercellular regulation

Additional components of bone remodeling modulation include intercellular regulation directly between osteoblasts and osteoclasts. Cells secreting osteoid, known specifically as Active Osteoblasts (AOBs), produce an essential regulatory factor, Receptor Activator of Nuclear Factor Kappa-B Ligand (RANKL). RANKL is the factor needed to activate the RANK receptor present on the osteoclast cells. The RANK/RANKL interaction will enhance osteoclast formation and activity, to drive the cycle of resorption/formation of bone remodeling [51-53]. Resting Osteoblasts (ROBs) secrete Osteoprotegerin (OPG) which is an antagonist to the RANK/RANKL interaction by binding to RANKL to prevent its interaction with RANK and thus protect the bone from over resorption [52]. In addition to this intercellular regulation, hormones provide an additional layer of checks and balances to the bone remodeling process in response to the needs of the whole body. We will discuss two of the more prominent and studied hormones, estrogen and Parathyroid Hormone (PTH). In addition to these two, there are many other hormones as well as growth factors that have been identified as having a role in the bone remodeling process such as growth hormone, glucocorticoids, and thyroid hormones [54]. Based on the percentage of the receptors for these hormones and factors which are occupied, the cells will elicit an appropriately graded response [55]. Estrogen has been found to stimulate osteoclast apoptosis through the action of TGF-ß [55]. Thus, the overall effect is to decrease bone resorption. Hence, in a disease such as osteoporosis, which is noted to have a high incidence in post-menopausal females, it correlates that an estrogen deficiency would result in a greater amount of bone resorption [56]. Another influential hormone is PTH which acts as an informant to osteocytes and osteoblasts regarding the state of calcium concentration levels in the body’s serum [55,57]. As a result of decreased calcium serum concentration, the parathyroid glands will up regulate the secretion of PTH. PTH will then act on its receptors among the bone cells to stimulate bone resorption in order to release calcium from mineralized bone and allow its entrance into the body’s serum [58,59]. The major mechanism of action for PTH is to bind to receptors on both active and resting osteoblasts in order to increase the production of RANKL and slow the production of OPG [55]. As a result, the number of active osteoclasts will increase, resulting in greater resorption. The dynamics of biological bone remodeling activity are under many levels of regulation from systemic control, such as hormones, as well as local control, such as RANKL and OPG. The influence of both biological and mechanical forces produces changes in bone as a summation of both types of stimuli.

Mathematical Modeling of Bone Remodeling Phenomena

A typical bone has demonstrates high strength and fracture toughness [60]. Macroscopic mechanical properties such as these are defined by bone mass, architecture, and cellular activities [61]. Thus, as previously discussed, mechanical and biological processes which are taking place throughout the hierarchy of bone structure at both the macroscale and microscale during the bone remodeling process all have effects on the overall architecture [62,63]. This reinforces the fact that bone is considered an anisotropic entity in which its apparent density and stiffness vary from point to point throughout the structure [64]. In addition to spatial dependency, Doblare, et al. make the distinction that these properties also vary with time [64]. With aging, the mechanical properties of bone change. Less mobilization leads to less external loads on whole bone, decreased hormones such as estrogen in women lower the stimulus for bone formation, and there may be an age-dependent decrease in osteocyte cells [65-68]. To study these spatiotemporal properties, mathematical models of the bone remodeling phenomena have provided great insight [6,69]. Traditionally, the results of these mathematical models have then been compared to data obtained from in vivo studies of bone properties [70]. However, thus far, the current models still possess many limitations which inhibit the models from truly replicating the integrative events taking place in living bone. In order to achieve a more realistic representation, a model must account for both the mechanical and biological processes taking place as well as the interscale synergies discussed prior. The following is a review of the current methods being utilized to construct mathematical models of the bone remodeling process. Some of the key bone properties able to be modeled as of yet include external mechanical load, bone density, stress, strain, entropy, chemical kinetics, and interstitial shear stress. It is important to note that trabecular bone undergoes a higher rate of bone remodeling than cortical bone; thus, making cortical bone more stable and thus responsible for the maintenance of bone properties over time [71]. The objective of this paper is to provide a sampling of the key progress made in the field thus far and commentary on how interscale synergies may be bridged so that the possibility and need for a mathematical model which encompasses the multiphysics, multiscale complexities is depicted.

External mechanical loads

We first look at a major stimulus of bone remodeling, the presence of mechanical loads on the external surface of the mineralized matrix of bone. An example of such loads is the forces applied during movements such as gait. Figure 4 shows an example of some of the mechanical forces applied to a femur by surrounding muscles [72]. In order to accurately model the mechanical loads being transferred from the muscles, the force, direction, angle, time of duration, and distance of effectiveness should be characterized and accounted for. As these mechanical loads are present at the same time, their individual effects must then be summated. Ideally, to fully characterize the external mechanical loads being applied to a femur, the properties of all muscles involved could be included.

Figure 4: A femur labeled with the forces of some of the muscle involved with gait [72].

    
    
    Figure 4:  A femur labeled with the forces of some of the muscle involved
with gait [72].

Probability of bone surface location being subject to remodeling

Huiskes, et al. introduces a concept to look at why and where bone remodeling takes place at specific sites. The basis of the idea states that bone remodeling is either occurring for one of two reasons. One is that it is randomly acting to repair microcracks as a part of regular bone maintenance. The other is that the process is being activated at specific locations which are under an increased amount of mechanical load. The probability of a bone surface location being subject to remodeling, is defined as

hypothesis I: p(x, t) = 10% (spatially random) (1) [1]

hypothesis II: p(x, t) = c[a- P(x, t)], if P < a (strain dependent)

This theory allows for the incorporation of both the major callsto- action of bone remodeling including bone maintenance and bone adaption. However, as previously stated, the incorporation of biological stimuli is again not being considered in this concept.

Bone density

Being able to predict changes in bone density may be able to provide insight into how to treat clinical ailments such as osteoporosis. A variety of factors previously discussed influence the formation or resorption of bone which have a direct impact on bone density. The change in total bone mass has been defined as the difference of the amount of bone mass accumulated and the amount of bone mass degraded by the formation and resorption processes, seen as

$\frac{d m_{t o t a l}}{d t} = A c c u m u l a t i o n - D e g r a d a t i o n = Δ F o r m a t i o n - Δ R e s o r p t i o n (2)$

Bone density is a unique property for mathematical models of bone remodeling to consider as it can be defined at the macroscale for total bone density (g/cm²) as well as for a defined smaller volume or local volume of bone at the microscale (g/μm²). A visual representation of a local volume of bone is shown in Figure 5 which is represented by a cube which can change density over time, represented by $\frac{d m}{d t},$ based on the effects of stimuli such as the load, FL, being applied as a result from the total external force, F_T (as seen in Figure 4), acting upon the femur.

Figure 5: a) A femur under the influence of total force, F_T (as seen in Figure 4). b) A cube representing a small, local, volume of bone from the femur which is under the influence of force of mechanical load, F_L, derived from F_T. The cube changes density over time, dm/dt as a result of the formation and resorption processes of bone remodeling based on the effects of stimuli such as the load F_L

    
    
    Figure 5:  a) A femur under the influence of total force, F_T (as seen in Figure
4). b) A cube representing a small, local, volume of bone from the femur
which is under the influence of force of mechanical load, F_L, derived from F_T.
The cube changes density over time, dm/dt as a result of the formation and
resorption processes of bone remodeling based on the effects of stimuli such
as the load F_L.

In current models such as Ruimerman, et al. [33], the total change in bone mass has been defined mathematically in a biological manner by incorporating the actions of two of the major bone cell types, osteoblasts and osteoclasts, defined as

$\frac{d m_{t o t} (x, t)}{d t} = \frac{d m_{b l} (x, t)}{d t} - \frac{d m_{c l} (x, t)}{d t} (3)$

This equation follows the work by Huiskes, et al. [1] who defined the relationship for the local change in bone mass, like for the cube seen in Figure 5, as

$\frac{d m}{d t} = τ {P (x, t) - k_{t r}} - r_{o c} f o r P (x, t) > k_{t r} (4)$

$\frac{d m}{d t} = - r_{o c} f o r P (x, t) \leq k_{t r}$

Huiskes, et al. also incorporates the use of cellular bone remodeling signaling by utilizing a decay function as recruitment stimulus for osteoblasts to become active and participate in bone formation, seen as

$P (x, t) = \sum_{i = 1}^{n} f_{i} (x) μ_{i} R_{t i} (t) (5)$

The decaying effect indicates that the signal has a specific effective distance range, exhibited by osteocytes in vivo, as mentioned previously; thus, the role of osteocytes and their signaling mechanisms in bone remodeling regulation is recognized.

Both Huiskes, et al. and Ruimerman, et al. make the assumption that osteoblast recruitment, ρ(x ,t), can be equated with bone formation and if the amount of bone mineral resorbed by osteoclasts, r_oc, is subtracted, the change in local bone mass over time can be measured [1,33]. It is not known if this assumption regarding bone formation is indeed true as it is unclear if every recruited osteoblast is activated to secrete osteoid. In addition, neither model considers the effects of biological stimuli such as the effects of hormones. Furthermore, in both models, the mechanosensitivity of the osteocyte, μ, is being considered as one value rather than accounting for the differences in mechanosensitivities exhibited among the osteocyte cell body and the osteocyte cell processes. In addition, there is no representation of the shear stress taking place at the lacuna canalicular network which is initiating the mechanical stimulus for osteocyte cell signaling to take place. By accounting for these elements, a more comprehensive multiphysics, multiscale model may be able to be derived.

A model which attempts to bridge from the microscale properties defined by Huiskes, et al. to the macroscale by relating their densities to one another is presented by Coelho, et al. [73]. An insightful representation of what is called the bone remodeling law, where it is assumed that structural stiffness will be maximized and metabolic cost, k, will be minimized is applied in order to achieve such a relationship [73]. In this model, it is recognized that the time scale of the bone remodeling process is much larger than the scale for an applied force to whole bone, accounting for the temporal dependency defined earlier. As a result, the forces are evaluated as a summation of P loads over time [73]. This design allows for the model to account for both processes, applied force and bone remodeling, in a more accurate manner. In order to carry out this two-scale approach, Coelho, et al. defined variables to macroscale properties of whole bone as well as microscale properties among the trabeculae of cancellous bone [73]. The macro-density, ρ, can be related to the micro-density, μ, as seen in the following

$ρ (x) = \frac{1}{| Y |} \int_{Y} μ (x, t) d Y, \forall x \in Ω (6)$

The bone remodeling law for this model, which accounts for the different time scales of mechanical loading and bone remodeling, is defined as

$\min \underset{ρ m i n \leq ρ (x) \leq ρ m a x}{ρ (x)} F (ρ) w i t h F (ρ) =$ $\int_{Ω} Φ (ρ, u^{1}, ..., u^{p}) d Ω + k (\int_{Ω} ρ (x) d Ω) (7)$

$w h e r e Φ (ρ, u^{1}, ..., u^{p}) = \max_{μ (y)} f (μ)$

$\int_{Y}^{μ_{m i n} \leq μ (y) \leq 1} μ (y) d Y = ρ (x)$

$w h e r e f (μ) = \sum_{r = 1}^{P} α^{r} [\frac{1}{2} E_{i j k l}^{H} (μ) ε_{i j} (u^{r}) ε_{k l} (u^{r})]$

This model provides a foundation for interconnecting the macroand micro- scales. However, it would still be of great interest to define total density in terms of both space and time, ρ(x ,t). It may be possible to add to this model the changes of density as a function of time through probability calculations from data collected by current technology using DEXA scans. In addition, the limitations of the model provided by Coelho, et al. lack the inclusion of bone cell populations and activities including biological stimuli and signaling as well as the incorporation of the lacuna canalicular system.

Mechanosensitivity and intercellular signaling

Huiskes, et al. only goes so far as to define the osteocyte’s mechanosensitivity as a general variable, μ. Adachi, et al. provide a model which specify an actual relationship for mechanosensitivity which incorporates the physical parameters of canaliculi as well as osteocyte processes, seen as

$S_{o c} (x) = \int_{0}^{2 π} d φ \int_{0}^{π / 2} \frac{2 r_{p}}{r_{c}^{2}} p_{c} (n) | \bar{τ_{p} (n)} | d θ (8)$

For simplicity, the cell processes are assumed to extend isotropically in the radial direction, thus the canalicular volume fraction would not depend on orientation n. Adachi, et al. go on to mathematically define the intercellular signaling communication between osteocytes and surface cells (such as osteoblasts and osteoclasts), written as

$S_{s f} (x_{s f}) = \int_{Ω} w (l) S_{o c} (x) d Ω (9)$

It is important to note that in addition to the other mentioned, the Adachi, et al. model does not account for the incorporation of biological stimuli such as hormones or the biological regulation through factors such as the RANK-RANKL system.

A model which does emphasize the additional role of biological element contribution (for example, from calcium or PTH) in the outgoing signal from an osteocyte, along with the thoroughly discussed mechanical stimulus is provided by Rieger, et al [57]. This model accounts for both the mechanical and biological components in the signal traveling to an osteocyte by adding a weight factor for each part, shown as

Ψ = WMecha < fMecha > +WBio < fBio > (10) [57]

This method offers a way to quantitatively measure the mechanobiological stimulus reaching the osteocytes. However, while this is a valid conceptual conclusion, a limitation of this model is that it is still not known what the appropriate physiologic weight factors are.

Interstitial shear stress

Interstitial shear stress along the osteocyte process is the cause of osteocyte cell activation. This stress is modified by properties of the canaliculi structure such as the size and permeability of these channels. As discussed previously, the osteocyte cell process has been proven to more sensitive to shear stress deformation than the cell body. Figure 6 shows the effect of an external mechanical force, F_T (from Figure 4), which is transduced into an internal mechanical load force, F_L (from Figure 5), in creating an interstitial shear force, τ, within the canaliculi of one of the osteocyte processes shown. The osteocyte will then take in these mechanical stimuli in addition to biological stimuli, and after a summation, the cell will produce an appropriate signal which will be transmitted via gap junctions to the laboring cells of bone remodeling, the osteoblasts and osteoclasts.

Figure 6: This is an expansion of the sequence displayed in Figure 5. This image now shows an in-depth look of an osteocyte present within the small volume of bone shown in b). c) depicts an osteocyte cell with several cell processes lying in the lacunocanalicular network (surrounding white space). It is in this lacunocanalicular network which is filled with interstitial fluid. The interstitial fluid changes its flow rate and thus shear stress on the basis of the mechanical forces present. Thus, in this case, the external mechanical force, F_T, seen in a) was transmitted as the internal mechanical load force, F_L, to the local volume of bone seen in b) and finally transduced as a shear stress in c).

    
    
    Figure 6:  This is an expansion of the sequence displayed in Figure 5. This
image now shows an in-depth look of an osteocyte present within the small
volume of bone shown in b). c) depicts an osteocyte cell with several cell
processes lying in the lacunocanalicular network (surrounding white space).
It is in this lacunocanalicular network which is filled with interstitial fluid. The
interstitial fluid changes its flow rate and thus shear stress on the basis of the
mechanical forces present. Thus, in this case, the external mechanical force,
F_T, seen in a) was transmitted as the internal mechanical load force, F_L, to the
local volume of bone seen in b) and finally transduced as a shear stress in c).

Adachi, et al. derive the interstitial shear stress being applied to osteocyte processes through the relationship stated as

$τ_{p} (n) = \frac{r_{c}}{γ} [A_{1} I_{1} (\frac{γ}{q}) - B_{1} K_{1} (\frac{γ}{q})] \nabla p (x) . (11)$

n where constants A₁ and B₁ are given by

$A_{1} = \frac{K_{0} (γ) - K_{0} (γ / q)}{I_{0} (γ / q) K_{0} (γ) - I_{0} (γ) K_{0} (γ / q)}$

$B_{1} = \frac{I_{0} (γ / q) - I_{0} (γ)}{I_{0} (γ / q) K_{0} (γ) - I_{0} (γ) K_{0} (γ / q)}$

Hambli and Korta also developed a model to simulate internal bone remodeling and the effects of interstitial fluid in the lacunocanalicular network. Their model is more progressive as it includes the velocity of the fluid as an additional contributing factor utilizing Bernoulli’s and Darcy’s law:

$V = - \frac{1}{2} \frac{K}{μ} \frac{ρ V_{1}^{2}}{L_{c a n}} [{(\frac{\frac{B}{A - (ε_{r r}) w 1} - R_{m}^{2}}{\frac{B}{A - (ε_{r r}) w 2} - R_{m}^{2}})}^{2} - 1] (12)$

Rate of change in interstitial bone fluid calcium concentration

Pennline recognizes the importance of calcium homeostasis in the regulation of bone formation and resorption, so he provides a definition of the rate of change of calcium in interstitial bone fluid, expressed as

$\frac{d C a}{d t} = - r_{F B} (B / B_{0}) + r_{B F} (C / C_{0}) + r_{F} + r_{C F} - r_{F C} - r_{d C a} C a (13)$

Developing more equations such as these for other important metabolic factors in the interstitial fluid, for instance, NO and PGE2 secreted from the osteocyte, would provide more information regarding the local environment concentration changes the osteocyte withstands and responds to.

Rate of change of cortical bone area

Pennline presents another property of bone that could be modeled and has the potential to reveal insight into the bone remodeling process [4]. Looking at a cross-section sample of cortical bone, Pennline first references Martinto define what will be considered as the intracortical bone balance, expressed as

Q = Q_BN_F - Q_CN_R (14) [75]

The rate of bone formation is then related as

$Q_{B} = \frac{π (R_{c}^{2} - R_{h}^{2})}{T_{F}} (15)$

and bone resorption is given by $Q_{c} = \frac{π R_{c}^{2}}{T_{R}} (16)$

Pennline then goes on to reference Martin’s definition of porosity, seen as

$p_{v} = \frac{A_{v}}{A} (17)$

All of these relationships are then used to establish an equation for the rate of change of cortical bone area in a cross-section sample, expressed as

$\frac{d A_{m}}{d t} = Q_{B} N_{F} (\frac{B}{B_{0}}) - Q_{C} N_{R} (1 - p_{v}) (\frac{C}{C_{0}}) (18)$

While the equation above can provide a suitable estimate of the change in bone area over time, its major limitation is that bone is being defined in a simplified form as a perfect cylinder. However, this is not the case. One way to better account for the real structure of bone may be to use Finite Element Analysis (FEA) for instance in order to maintain the actual shape and 3D elements of bone. Also, in addition to accounting for the change in concentration of osteoblasts and osteoclast cell types, osteocyte concentration could also be incorporated. With these changes in cell type densities, it would be important to consider the rate of cell apoptosis as in the model developed by Buenzli [76].

Stress and strain

From the induction of mechanical loads on the bone, stress and strains arise within the mineralized matrix and serve as stimuli to the bone remodeling process [77]. For this last section, we will present a look into how these bone properties can be approached from a mechanobiological perspective within the multiscale system. For example, consider a small area, ΔA, around P of the cutting plane of a deformed configuration shown in Figure 7. Let the internal resultant forces be represented by ΔF acting on ΔA.

Figure 7: a) Simplified image representing the diaphysis of a bone. Different forces being placed upon the bone are represented by F₁, F₂, F₃, and F₄. The blue line represents the cutting plane with a center Point P. b) Same representation now cut along the plane with the stress tensor, t_n, normal stress, σ_n, and the shear stress, σ_s.

    
    
    Figure 7:  a) Simplified image representing the diaphysis of a bone. Different
forces being placed upon the bone are represented by F₁, F₂, F₃, and F₄.
The blue line represents the cutting plane with a center Point P. b) Same
representation now cut along the plane with the stress tensor, t_n, normal
stress, σ_n, and the shear stress, σ_s.

The stress is a measure of the intensity of internal forces generated within the body of the cross sectional area and will vary from point to point. We define the stress tensor at Point P as:

$t_{n} = {lim}_{Δ A \to 0} \frac{Δ F}{Δ A} (19)$

With the normal stress defined as

$σ_{n} = {(t_{n})}_{i} {\bar{n}}_{i} \infty (20)$

and the shear stress defined as

$σ_{s} = {[{| t_{n} |}^{2} - {(σ_{n})}^{2}]}^{1 / 2} (21)$

In addition, due to the law of action and reaction, on the opposite side of the cut plane, we will have -tn, -sn, and -ss. This can be further generalized to 3 planes with

tn = t_xn_x + t_yn_y + t_zn_z(22)

and can be simplified with the use of Cauchy’s stress tensor relation ${(t_{n})}_{i} = σ_{i j}^{T} n_{j} (23)$

Now, if we proceed to a smaller scale and examine a point on the trabecular bone at an arbitrary cross sectional area and assume that the accumulation vector of the deformed configuration is given by

$\bar{a} = a_{x \vec{i}} + a_{y \vec{j}} + a_{z \vec{k}} (24)$

and let the external force per unit mass of the bone be given by

$\bar{f} = f_{x i}^{- 1} + f_{y j}^{- 1} + f_{z k}^{- 1} (25)$

Then, by using the balance of linear momentum, we can now write the following to define bone density in terms of stress:

$\sum_{i = 1}^{3} ρ a_{i} = ρ f_{i} + \nabla \cdot σ (26)$

${(\nabla \cdot σ)}_{x} = \frac{δ σ_{x x}}{δ x} + \frac{δ σ_{y x}}{δ x} + \frac{δ σ_{z x}}{δ x},$ ${(\nabla \cdot σ)}_{y} = \frac{δ σ_{x y}}{δ y} + \frac{δ σ_{y y}}{δ y} + \frac{δ σ_{z y}}{δ y},$ ${(\nabla \cdot σ)}_{z} = \frac{δ σ_{x z}}{δ z} + \frac{δ σ_{y z}}{δ z} + \frac{δ σ_{z z}}{δ z} (27)$

which can be simplified to

$ρ = (\frac{\nabla \cdot σ}{a - f}) (28)$

Next, we consider the strain which can be defined by

$[ε] = \frac{1}{2} {([\nabla_{o} u] + [\nabla_{o} u])}^{T} (29)$

$\in_{i j} = \frac{1}{2} (u_{i, j} + u_{j, i}) (30)$

Note, that in most mechanical systems ∇ºis of the order of 10_-6, but in this case, ∇º is of the order 10^-9 or less. There are 3 sets of equations that govern the displacements, strains, and stresses in an element. Assuming small deformations and linear elastic material, we have in an index notation

$ε_{i j} = \frac{1}{E} [- v σ_{k k} δ_{i j} + (1 + v) σ_{i j}] (31)$

With

$E = \frac{μ (3 μ λ + 2 μ)}{λ + μ}$

and Poisson’s ratio

$v = \frac{λ}{2 (λ + μ)} (32)$

Linear Rotational Formulation: If we assume that the constitutive equation above in terms of stress-strain relation are applied to a deformable body, assuming an isotropic material to simplify the calculations, we obtain

$\in_{i j} = \frac{1 + v}{E} σ_{i j} - \frac{v}{E} σ_{k k} δ_{i j} + α (T - T_{o}) δ_{i j} (33)$

Viscoelastic Behavior Formulation: There are several scenarios that explain the bone formation, resorption, and its dynamic response to load. With time, as in stress relaxation, the energy to increase bone formation also decreases. So at the initial state of bone formation, a highforceis imposed and maintained until the bone stress begins to decrease. Furthermore, if the stress maintained is constant and the bone strain increases, a phenomenon close to creep, then one has to assume that the bone formation can continue for a period of time, under the same stress or loading conditions. The rate of bone formation is directly associated with the velocity as in the case of a Maxwell model,

$\dot{u} = \frac{\dot{F}}{μ} + \frac{F}{η} (34)$

The constitutive equations for numerical simulation for at least 2D models result in

${\int_{v} (\nabla δ u)}^{T} σ d v - \int_{s} δ ϱ^{T} t d S = \int_{d v} δ u^{T} f d S (35)$

Chemical Virtual Work Traction acting at the (stress-strain) Boundary of the cell

There are two sources of elastic response to deformation: change of internal energy and change of entropy. Assuming we have conservation of energy, we can examine how the laws of thermodynamics connect to internal energy: Ρ (energy per unit mass), S (entropy per unit mass), T (temperature), p (pressure), V (volume), ρ (density), σ_ij (stress), ∈^ij (strain), and σ_ij ‘∈_ij’ (stress-strain deviations). If we let dQ be the energy transferred to the body, then

$dε = dQ + \frac{1}{ρ} σ_{ij} d \in_{ij} (36)$

dQ = T dS (37)

$d ε= T d S + \frac{1}{ρ} σ_{ij} d \in_{ij} (38)$

With $σ_{i j} = - p δ_{i j} + σ_{i j}^{'}$ and $\in_{i j} = \frac{1}{3} \in_{x x} δ_{i j} + \in_{i j}^{'}$

$d ε = T d S - p d V + \frac{1}{ρ} σ_{i j}^{'} d \in_{i j}^{'}$

If the material is perfectly elastic, then the existence of a strain-energy function can often be justified on the basis of thermodynamics. However, living tissues are not perfectly elastic, and therefore, they cannot have an energy function in terms of strain-energy function in the thermodynamics sense. If we can ensure that the rate of bounding and unbounding is such that the strain rate is small and can be ignored, as is the case for bone remodeling, than the stress-strain can be uniquely defined.

Mapping Functions and Bone Remodeling

The cells live under continuous stress which causes them to participate and engage in a metabolic internal reaction as well as protein resorption and mass change due to the energy being diminished and produced. While preserving the constitutive laws of physics, the mechanical change in terms of mass, geometry, play a role in the local and global space and time frame (aging process).

The constitutive laws are directly related through variables:

$σ_{i j} = σ (\in_{i j}, ρ_{i}, f_{i}, μ_{i} ...) (40)$

$ρ_{i} (t) = ρ_{i} (t = 0) + \int_{0}^{t} ρ_{\dot{i}} (\hat{ι}) d \hat{ι} (41)$

$\overset{\cdot}{ρ} = G (ρ_{i}, f_{i}, \in_{i j}, σ_{i j}, o t h e r f a c t o r s) (42)$

$\overset{\cdot}{f} = H (ρ_{i}, f_{i}, \in_{i j}, σ_{i j}, ...) (43)$

The cross-bridge concept of creating mapping functions is to help us understand the inter-relationship between the biological, mechanical functions can be defined using probabilistic finite elements, where the cross bridge is at a finite point in space and time (n,t). The Lagrangian stress can be appreciated with the instantaneous energy and mass as it relates to stiffness. The bounding and unbounding of matters and energy are probabilistically determined based on a set of variables that are determined and defined by

$n (ε, t) = Q^{- 1} \int_{0}^{t} f [(ε - δ) (τ)] Q (ε, t) d Z + n_{0} (ε) Q^{- 1} ε (t) \frac{n!}{r! (n - r)!} (44)$

where f and Q are functions describing the bounding and unbounding of energy processes.

What remains to be defined is the extraction of these functions for different mechanisms including chemical reactions, transport phenomenon, energy dissipation, and the creation of matter.

Developing a Comprehensive Model

The goal of this research was to identify the applicability, needs, and limitations in building a comprehensive multiphysics, multiscale model which would ultimately be able to predict age-related changes in bone remodeling. New developments with current technology such as μ CT allow for imaging and quantitative 3D visualization of smaller scale structures of bone such as osteons and the lacuno canalicular network [78-82]. By being able to obtain more defined and realistic properties of these structures, mathematical models will greatly improve [83]. Animal and clinical studies such as Webster, et al. which was able to calculate strain energy density gradients of mice vertebrae will offer methods for obtaining these values [84]. Ideally, it is best to not only incorporate bone’s properties and interactions but also those properties and interactions of surrounding supporting structures, such as muscles and ligaments [85]. In addition to accounting for all of these essential factors, the model needs to recognize the spatialtemporal dependency of the bone remodeling process [86]. In Figure 8, we propose a demonstration of the co-dependency of the time scale of the bone remodeling process on the spatial structure of the bone matrix.